In search of a data-driven symbolic multi-fluid ten-moment model closure

Abstract

The inclusion of kinetic effects into fluid models has been a long standing problem in magnetic reconnection and plasma physics. Generally, the pressure tensor is reduced to a scalar which is an approximation used to aid in the modelling of large scale global systems such as the Earth's magnetosphere. This unfortunately omits important kinetic physics which have been shown to play a crucial role in collisionless regimes. The multi-fluid ten-moment model, however, retains the full symmetric pressure tensor. The ten-moment model is constructed by taking moments of the Vlasov equation up to second order, and includes the scalar density, the vector bulk-flow and the symmetric pressure tensor for a total of ten separate components. Use of the multi-fluid ten-moment model requires a closure which truncates the cascading system of equations. Here we look to leverage data-driven methodologies to seek a closure which may improve the physical fidelity of the ten-moment multi-fluid model in collisionless regimes. Specifically, we use the sparse identification of nonlinear dynamics (SINDy) method for symbolic equation discovery to seek the truncating closure from fully kinetic particle-in-cell simulation data, which inherently retains the relevant kinetic physics. We verify our method by reproducing the ten-moment model from the particle-in-cell (PIC) data and use the method to generate a closure truncating the ten-moment model which is analysed through the nonlinear phase of reconnection.

1 Introduction

Magnetic reconnection is the process in which a magnetic field embedded in plasma undergoes a topological restructuring. This process, which converts the energy stored in magnetic fields into particle kinetic and thermal energy, is ubiquitous throughout the universe and plays an important role in such diverse events as sawtooth crashes in fusion plasmas (Zweibel & Yamada 2009; Yamada, Kulsrud & Ji 2010), magnetic substorms in the Earth's magnetosphere (Hastie 1997) and solar coronal mass ejections (Masuda et al. 1994).

Magnetic reconnection was first studied using magnetohydrodynamics (MHD). In ideal MHD, which describes a plasma as a single fluid of infinite conductivity, Ohm's law in the lab frame is given by $\boldsymbol {E} + \boldsymbol {u}\times \boldsymbol {B} = \eta \boldsymbol {J} = 0$ and implies that magnetic flux is frozen into the plasma flows. Hence, magnetic reconnection is topologically prohibited.

In real plasmas, however, the right-hand side of Ohm's law is not zero. In resistive MHD, finite conductivity is introduced as the resistive term $\eta \boldsymbol {J}$. Parker (1957) used this model to develop the first self-consistent description of magnetic reconnection. However, many plasmas of interest are nearly collisionless and the reconnection rates predicted by Sweet and Parker do not match with observations. Petschek's model (Petschek 1964) addressed the shortcomings of the Sweet–Parker model by realizing that the geometry of the reconnection region is crucial and that reducing the aspect ratio of that region can explain the observed faster reconnection rates. This means that the reconnection region does not span the global length scale but is much shorter and connected to the global geometry by slow shocks. Petschek's model, however, is not realized in numerical simulations unless one triggers a shorter reconnection region by anomalous resistivity or other modifications of the physics.

It was later realized that the one-fluid MHD description of plasmas is too limited to incorporate all the essential reconnection physics. The next steps were Hall-MHD and extended MHD models that take into account that electrons and ions do not move together at the small scales of a reconnection region. From the electron momentum equation, one can derive a generalized Ohm's law (Vasyliunas 1975; Birn et al. 2001):

(1.1)\begin{equation} \boldsymbol{E} + \boldsymbol{u}\times\boldsymbol{B} = \eta\boldsymbol{J} + \frac{\boldsymbol{J}\times\boldsymbol{B}}{n|e|} - \frac{\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{P}_e}{n|e|} + \frac{m_e}{n|e|^2} \times \left[\partial_t \boldsymbol{J} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left(\boldsymbol{u}\boldsymbol{J} + \boldsymbol{J}\boldsymbol{u} - \frac{\boldsymbol{J}\boldsymbol{J}}{n|e|} \right) \right]. \end{equation}

The generalized Ohm's law (1.1) describes the various non-ideal contributions to the electric field in a plasma. The resistive dissipation term, sometimes denoted as the collisional term, is given by $\eta \boldsymbol {J}$. The second term on the right-hand side is called the Hall term which, in and of itself, cannot support magnetic reconnection. It is, however, known to change the geometry of the reconnection process to enable ‘fast’ reconnection, where the time scales of the process are of the order of the Alfvenic transit times. The third term is the divergence of the electron pressure tensor, and the last term is the electron inertial term.

In fluid simulations, the electron pressure tensor is usually approximated by a scalar pressure, however, it is known from fully kinetic simulations that off-diagonal terms can support reconnection electric fields. It has generally become clear that kinetic physics, being the first-principles description of a plasma, are important for a complete description of magnetic reconnection, in particular, in the collisionless regime. Advances in computing power have enabled comparatively large fully kinetic simulations using the particle-in-cell (PIC) method, even in three dimensions, that can achieve a good separation between the global size of the simulation, and ion and electron scales, even though they usually require a reduced ion/electron mass ratio to keep the computational cost manageable. These kinds of simulations are limited to local simulations of reconnection, although their extent may be 10s or 100s of ion inertial scales ($d_i$). In contrast, many physical systems of interest, e.g. Earth's magnetosphere or the Sun's corona, are many orders of magnitude larger and will remain inaccessible to fully kinetic simulations for the foreseeable future.

Fluid simulations are computationally significantly cheaper than kinetic simulations due to their reduced dimensionality, making them computationally more attractive but at the cost of a loss of kinetic information in the closure formulation.

An alternative fluid formulation employed in this context is the multi-fluid moment model. This model is derived without any approximations by taking moments of the Vlasov equation and keeps the full set of Maxwell's equations.

1.1 Multi-fluid and ten-moment model

The multi-fluid model is constructed by taking sequentially increasing velocity space moments of the Vlasov equation:

(1.2)\begin{equation} d_t f_s = \partial_t f _s+ \boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla}_r f_s + \frac{q}{m}\left( \boldsymbol{E} + \boldsymbol{v}\times \boldsymbol{B}\right) \boldsymbol{\cdot} \boldsymbol{\nabla}_v f_s = 0 \end{equation}

for each species $s$. Following Wang et al. (2015), the following moments are defined:

(1.3a)\begin{gather} n\equiv \int f \,{\rm d}\boldsymbol{v}, \end{gather}

(1.3b)\begin{gather}u_j \equiv \frac{1}{n} \int v_j f \,{\rm d}\boldsymbol{v}, \end{gather}

(1.3c)\begin{gather}\mathcal{P}_{ij}\equiv m\int v_iv_j f \,{\rm d}\boldsymbol{v}, \end{gather}

(1.3d)\begin{gather}\mathcal{Q}_{ijk}\equiv m\int v_iv_jv_k f \,{\rm d}\boldsymbol{v}. \end{gather}

The evolution equations for these moments are derived by multiplying the Vlasov equation by consecutive powers of $v$ and integrating out velocity space. Truncating the system of equations after the second order gives the so-called ten-moment equations (1.4), which describe the evolution of density, the three components of momentum and each of the six unique components of the symmetric pressure tensor. Each species is maintained as a separate set of fluid moments, but a species index has been left out above for brevity.

(1.4a)\begin{gather} \partial_t n + \partial_{j} (nu_j) = 0, \end{gather}

(1.4b)\begin{gather}m\partial_t(nu_i) + \partial_{j}\mathcal{P}_{ij} = nq\left( E_i + \epsilon_{ijk}u_jB_k\right), \end{gather}

(1.4c)\begin{gather}\partial_t \mathcal{P}_{ij} + \partial_{k} \mathcal{Q}_{ijk} = nqu_{[i}E_{j]} + \frac{q}{m}\epsilon_{[ikl}\mathcal{P}_{kj]}B_l. \end{gather}

The square bracket notation above indicates the minimal sum over free indices that yields a completely symmetric tensor. The electromagnetic fields evolve according to the Maxwell equations:

(1.5a)\begin{gather} \boldsymbol{\nabla} \times \boldsymbol{E} ={-}\frac{\partial\boldsymbol{B}}{\partial t}, \end{gather}

(1.5b)\begin{gather}\boldsymbol{\nabla} \times \boldsymbol{B} = \mu_0 \boldsymbol{J} + \frac{1}{c^2}\frac{\partial\boldsymbol{E}}{\partial t}. \end{gather}

The moments in (1.3) are used for the derivation of the ten-moment model given by (1.4). Frequently, however, working with the centeredmoments (1.6) is preferred because of their convenient physical interpretation as the pressure stress tensor and heat flux tensor.

(1.6a)\begin{gather} P_{ij} \equiv m\int (v_i - u_i)(v_j - u_j) f \,{\rm d}\boldsymbol{v}, \end{gather}

(1.6b)\begin{gather}Q_{ijk} \equiv m\int (v_i - u_i)(v_j - u_j)(v_k - u_k) f \,{\rm d}\boldsymbol{v}. \end{gather}

The ten-moment model may be re-written in terms of the centeredmoments, which are related according to

(1.7a)\begin{gather} \mathcal{P}_{ij} = P_{ij} + 2nmu_iu_j, \end{gather}

(1.7b)\begin{gather}\mathcal{Q}_{ijk} = Q_{ijk} + u_{[i}\mathcal{P}_{jk]} - 2nmu_iu_ju_k. \end{gather}

Noticeably, each moment equation contains the next higher order moment. This trend continues ad infinitum resulting in an open system of equations. While this system of equations is exact, to be of practical use, it needs to be truncated. We truncate the ten-moment model by seeking a closure that replaces $Q$ with an expression containing known quantities, such as the lower order moments. In the adiabatic case, $P_{ij}$ is taken to be isotropic and $Q$ is zeroed out, giving the five-moment model. In this work, we are specifically interested in the pressure tensor's impact on reconnection and hence choose to retain its full structure.

Previous closures have been proposed such as the CGL family of closures (Chust & Belmont 2006) and the Hammett–Perkins closure (Hammett & Perkins 1990).

In this work, we aim to use data-driven methodologies to derive a symbolic closure relation truncating the cascading system of equations at the ten-moment model.

We build off the work of Wang et al. (2015), investigating the collisionless regime using data from a fully kinetic Harris sheet PIC simulation. We then apply a data-driven technique, sparse identification of nonlinear dynamics (SINDy) (Brunton, Proctor & Kutz 2016), to distill the raw particle output into a symbolic closure which is compared to the approximate local closure:

(1.8)\begin{equation} \partial_m Q_{ijm} \approx v_t |k_0|\left(P_{ij} - p\delta_{ij} \right)\!. \end{equation}

In the above, $v_t$ is the thermal velocity, $k_0$ is a typical scale-defining wave-number, $p$ is the scalar pressure attained by averaging the diagonal of the pressure stress tensor and the moments are centred. Originally defined by Wang et al. (2015), (1.8) is referred to as the local approximate Hammett–Perkins closure because it replaces the continuous $k$ in the general non-local Hammett–Perkins closure, with a single $k_0$ value in physical space.

In the literature, various options for closures to the ten-moment model have been investigated (Ng et al. 2017, 2020) by implementing a given proposed closure into the multi-fluid moment model, choosing parameters like a typical wavelength of kinetic instability $k_0$ and running prototypical reconnection problems like the Harris sheet or island coalescence. The time evolution and snapshots of the fields are then compared to results of a PIC simulation.

In this work, however, we exploit the additional information available from PIC simulations to test and derive closures. Since we do have the fully kinetic particle information available, we can, for example, calculate the third-order heat flux moment and compare the actual heat flux observed in the simulation to the assumptions of a given closure. This process has the advantage that we can directly compare terms at any given time, whereas when comparing multi-fluid and PIC simulations, the state at any given time is an accumulation of all differences that have occurred up to this time.

2 Background

2.1 Particle-in-cell method

The particle simulation code (PSC) (Germaschewski et al. 2016) is a modern, load-balanced, GPU accelerated, fully kinetic PIC code. The PIC method is a numerical method which discretizes the electromagnetic fields on a 3-D grid and advances them using Maxwell's equations while approximating the particle distribution function through macroparticles which are advanced in continuous phase space.

As a fully kinetic algorithm, the PIC method solves the full Vlasov–Maxwell system of equations. Thus, moments of the distribution function evolved by the PIC method exactly satisfy the ten-moment model up to noise and other numerical errors. Systemic noise inherent to the PIC method is introduced primarily through the sampling of macroparticles to construct distribution functions.

This systemic noise was demonstrated to be a significant problem when applying the SINDy method to PIC data by Alves & Fiuza (2022). By spatially integrating the PIC data, they were able to reduce noise and apply SINDy to effectively recover the Vlasov equation, recover fluid equations, and discover an adiabatic closure. In § 3, we borrow this technique and apply it to each of the moments and fields output by the PSC before any further operation is performed.

2.2 Sparse identification of nonlinear dynamics

SINDy is a framework for the discovery of symbolic equations. Developed by Brunton et al. (2016), it looks to construct parsimonious solutions from a library containing candidate terms which are calculated from raw data. It may be used to take data from either experiment or simulation and generate symbolic governing equations which describe a dynamical system.

This is accomplished by solving the sparsity-promoting regression problem

(2.1)\begin{equation} \mathcal{L} = \min_{\boldsymbol{w}} \|\varTheta\boldsymbol{w} - \boldsymbol{y} \|_2 + \lambda R(\boldsymbol{w}), \end{equation}

where $\varTheta$ is the library containing candidate terms, $\boldsymbol {w}$ is the vector of coefficients associated with each term in the library, $\boldsymbol {y}$ is the regression target and $R(\boldsymbol {w})$ is some regularization function on $\boldsymbol {w}$ scaled by constant $\lambda$. Commonly, $R$ is the $L_1$-norm in which case the above becomes lasso regression. For a detailed description of the framework, see Zheng et al. (2019) and Champion et al. (2020).

Our experiments approximate (2.1) by solving the least squares problem with thresholding. This involves iteratively solving the least squares problem $y=\varTheta \boldsymbol {w}$ and applying an upper and lower bounded threshold to the discovered $\boldsymbol {w}$, thus restricting the coefficients to a pre-specified range.

The success of SINDy discovering a parsimonious equation is dependent on selecting upper and lower bounds which maximally encourages sparse solutions while minimally allowing error. Solutions to problems of this sort are known as pareto-optimal if the family of solutions cannot improve on one measure without hindering the other. Therefore, the optimal bounds for the thresholded least squares problem are associated with the solution that lives on the pareto-front where the number of included terms may not be further minimized without increasing the solution error.

Thus, for each instance in which we apply SINDy, we first solve the thresholded least squares problem with various bounds to construct the pareto curve. The optimal bounds are then selected from the pareto-front and the problem is re-solved with those bounds.

This method was selected to symbolically model our closure because it is fully interpretable, emulates nature by leveraging parsimony, and has previously been used successfully to model both kinetic and fluid plasma systems (Kaptanoglu et al. 2021; Alves & Fiuza 2022).

3 Experiments

3.1 Harris sheet simulation

This study uses PSC to simulate a collisionless Harris sheet reconnection problem. Initially in kinetic equilibrium, the magnetic field is initialized as $\boldsymbol {B} = B_0\tanh (y/\lambda _B)\widehat {\boldsymbol {x}}$ and the densities as $n_e=n_i=n_0\operatorname {sech}^2(y/\lambda _B) + n_b$. The magnetic field is then given a small sinusoidal perturbation which initiates the reconnection process.

The run uses a $1280 \times 640$ grid with a mass ratio of 25 and $16\times 10^9$ particles which resolves to $0.676\Delta x$/Debye length.

This simulation uses an abnormally high number of particles per cell for each species, approximately 10 000. This was done to lower the level of noise inherent to PIC methods, as demonstrated by Juno et al. (2020) (see table 1).

Table 1. Harris sheet simulation numerical parameters.

3.2 Method

Our method, which relies on SINDy at its core, was used for verifying the ten-moment model, analysing the existing local approximate Hammett–Perkins closure and searching for an improved closure.

We begin by loading the moment and field data generated by the PSC, and calculating relevant spatial and temporal derivatives. We then spatially integrate each term and construct a library of terms consisting of every right-hand side term from the ten-moment model.

In verifying the ten-moment model, each order's term library $\varTheta$ is constructed by selecting out only terms with consistent units. Specifically, the zeroth-order's term library $\varTheta$, contains only two terms from from the single zeroth-order equation (reducing the method to standard regression), the first-order's term library contains nine possible terms from the three first-order equations, and the second-order's term library contains 24 possible terms from the six second-order equations. Pairing terms with consistent physical units acts as an a priori constraint on which terms may be accepted into the library.

With the constructed library, we then calculate the left-hand side of each equation in the ten-moment model. For consistency, we consider the left-hand side to be the derivative containing terms of the equation, with the exception of the zeroth order, where the partial time derivative is the left-hand side and the particle divergence is the right-hand side.

The left-hand side is then used as a regression target wherein we apply the SINDy framework on the associated term library to discover the associated symbolic equations constituting the right-hand side. This method is summarized in Algorithm .

Algorithm 1 Equation Synthesis

3.3 Ten-moment model and method verification

As the multi-fluid moment equations (1.4) are exact, we can verify them directly from the Harris sheet PIC data. To show this, we individually calculate the left-hand side and right-hand side of each equation and demonstrate their equivalence.

In the following analysis, we use the antisymmetric normalized $L_2$ error (3.1) as a method for four separate comparisons. The first is the numerical error between the left-hand side and right-hand side calculated directly from the PIC simulation data. This metric is used to verify that the kinetic data satisfies the ten-moment model and gives an indicator of the residual noise in the data after integration. The second comparison is the numerical error between the right-hand side calculated directly from the data and the discovered right-hand side equation applied to the data. This is an indicator for how well our method is able to reproduce theory. The third comparison is the numerical error between the left-hand side calculated from the data and the discovered right-hand side equation applied to the data, which gives an indicator for how well our discovered equation can explain the left-hand side of each moment relation.

The final comparison we make is the normalized $L_2$ metric between the coefficients given by the known multi-fluid equations (1.4) and those discovered using the regression method. This is referred to as the coefficient error.

(3.1)\begin{equation} L_2(\boldsymbol{x}_1,\boldsymbol{x}_2)= \frac{\|\boldsymbol{x}_1 - \boldsymbol{x}_2\|_2}{\|\boldsymbol{x}_1\|_2}. \end{equation}

Using this analysis, we find that the regression method is able to reproduce the multi-fluid equations if the $L_2$ (left-hand side, discovered) error is equal to or approximately equal to the $L_2$ (left-hand side, right-hand side) error. Ultimately, we are looking for symbolic not numeric closure approximations. With this in mind, the coefficient error is the truest indicator for the success of the method in reproducing (1.4).

We present our findings in tables 2–4 and figures 1–3, where the data presented are from a representative step T-17 $\omega _{ci}^{-1}$ (ion cyclotron frequency), during the nonlinear phase of reconnection. Experiments were run at various time steps throughout the nonlinear phase and the results were shown to not differ greatly.

Figure 1. The zeroth-order moment equation or the continuity equation verified from the Harris sheet PIC data. Here, we contrast the left-hand side in panel (a) with the true right-hand side in panel (b) and the discovered right-hand side in panel (c).

Figure 2. Top to bottom are the $x$-component, $y$-component and $z$-component of the first-order (momentum) equation with the left-hand side in panels (a,d,g), the true right-hand side in panels (b,e,h), and the discovered right-hand side in panels (c,f,i).

Figure 3. Top to bottom, left to right are the xx, yy, zz, xy, xz, yz components of the second-order moment equation. Each contains the regression target in panels (a,b,g,h,m,n), the true right-hand side in panels (c,d,i,j,o,p), and the discovered relation in panels (e,f,k,l,q,r). Visually, the left-hand side, right-hand side and discovered source term of each component match well and are consistent with theory. This is true of both the on-diagonal and off-diagonal terms, which demonstrates the reliability of our data. While some off-diagonal discovered terms may be missing, this is most likely due to their magnitude being small in this particular scenario.

Table 2. The true right-hand side, discovered right-hand side and systematic errors of the zeroth-order moment equation.

Table 3. The true right-hand side, discovered right-hand side and systematic errors for each component of the first-order moment equations.

Table 4. The true right-hand side, discovered right-hand side and systematic errors for each component of the second-order moment equations.

Systematic error is introduced into the analysis when approximating the time and space derivatives with the centred finite-difference method. This error appears in the $L_2$ (left-hand side, right-hand side) and the $L_2$ (left-hand side, discovered) errors presented in the analysis and does not indicate a violation of the conservation laws.

There is potential that this source of error could pose a problem in discovery of derivative-containing terms. So while this would not be a problem for the method verification of § 3.3, where the term library $\varTheta$ does not contain derivative terms, it could raise issues for closure discovery in § 3.5. To investigate whether this error was a problem, the method verification of § 3.3 was repeated for each equation using only the partial time derivative term as the regression target. This modification required the spatial derivative terms be included in $\varTheta$. Results showed that noise induced by the spatial derivatives did not pose a problem in right-hand side discovery. The derivative terms were found with error of the order of the results as presented.

With semantic consistency in mind, we thus stick to the convention of the ten-moment model and use the derivative-containing terms as the left-hand side throughout the following sections.

3.3.1 Zeroth-order verification

The zeroth-order equation, calculated from kinetic data, is verified visually, as seen by the close matching in figure 1 and the low $L_2$ (left-hand side, right-hand side) error. One might believe that the residual noise would distort the ability of SINDy to isolate the correct dynamics. This, however, is not the case, as demonstrated by the low coefficient error, 0.065. This demonstrates the robustness of this method to residual noise and the derivative approximation systematic error. This is important for validating the approach for use in closure discovery where many library terms will contain derivatives.

3.3.2 First-order verification

The inclusion of field terms in the first-order equations stabilizes the numerics significantly. This is demonstrated by the superb left-hand side–right-hand side matching, the low $L_2$ (left-hand side, right-hand side) error and the low $L_2$ (left-hand side, discovered) error for each of the three components. Similarly, from the coefficient error, we see that the method was able to reproduce the first-order equations with 0.01–0.03 normalized $L_2$ coefficient error. These are strong results which validate both the use of kinetic data for use in the ten-moment model and the use of SINDy to discover dynamics using the data.

3.3.3 Second-order verification

The second-order equations demonstrate different behaviour on-diagonal versus off-diagonal. The $L_2$ (left-hand side, right-hand side) numerical errors are low in both the on-diagonal and off-diagonal terms, validating the second-order equations of the ten-moment model. Similarly, the low $L_2$ (right-hand side, discovered) in conjunction with the low $L_2$ (left-hand side, discovered) numerical error indicates that the terms included in the discovered solution account for the major contributions to the true right-hand side for each component. This similarly implies that the terms excluded from the discovered right-hand side must be of small magnitude. The small magnitude of these terms leads to discrepancies in the discovered right-hand side of the xy, xz and yz directions, where the coefficient error is $> 0.4$. This high coefficient error is due to the exclusion of terms from the discovered solutions. Any term with a small magnitude will be eliminated by SINDy and is a constraint of the method. Thus, this inconsistency is explainable and expected in the case of minor-contributing terms.

3.4 Local approximate Hammett–Perkins closure

The local closure (1.8) replaces the continuous wave-number $k$ in the Fourier space response with just one typical wave-number $k_0$. Thus, omitting $k_0$ when calculating the closure, we would expect the left-hand side and right-hand side to differ by a constant factor. Further, the closure is most important near the reconnection x-point where we expect important departures from ideal behaviour. Generally, these trends seem to be verifiable in the kinetic data seen in figure 4; there is, however, room for improvement in several areas. One obvious discrepancy is that each direction of the heat-flux divergence tensor scales by a different $k_0$ from the closure. Second is the clear disparity between the heat-flux divergence and the closure at regions of the domain distant from the x-point.

Figure 4. The local approximate Hammett–Perkins closure calculated during the nonlinear phase of reconnection. In each coordinate, panels (a–c,g–i) show the divergence of the heat flux, while panels (d–f,j–l) represent the prediction of the approximate closure. The unknown factor $k_0$ for each component has been calculated by averaging ${\boldsymbol {\nabla } Q_{ij}}/({v_t(P_{ij} - p\delta _{ij})})$ across the entire domain.

Each component presented in figure 4 is represented with $k_0$ calculated by taking the average ${\boldsymbol {\nabla } Q_{ij}}/({v_t(P_{ij} - p\delta _{ij})})$ across the entire domain. In all cases, this results in the pressure terms being washed out with a lower magnitude than the heat flux divergence. The on-diagonal terms appear to match structurally while the off-diagonal terms differ. Here, $\boldsymbol {\nabla } Q_{xy}$ has some quadrupolar structure while the pressure terms appear to be bipolar. The large scale structures of $\boldsymbol {\nabla } Q_{xz}$ and $\boldsymbol {\nabla } Q_{yz}$ are represented in the associated pressure terms, but the fine-scale structures are omitted.

These discrepancies are significant and can lead to issues with physical fidelity when put into practical use as part of the ten-moment model.

3.5 Closure discovery

Applying our methodology to the improvement or replacement of the above given closure, we decided to restrict ourselves to the nonlinear phase of the reconnection process. This range may be seen in figure 5. Pareto-optimal bounds are required for SINDy to discover a parsimonious solution. Constructing a separate pareto curve for each element of $\boldsymbol {\nabla } Q_{ij}$, we modified Algorithm  to that shown in Algorithm . A representative pareto curve used in discovery of the yy-component closure is shown in figure 6.

Figure 5. The L-2 norm of each component of the heat flux divergence tensor, taken across the entire simulation domain. Vertical bars represent the nonlinear phase of reconnection.

Figure 6. A representative pareto curve constructed by solving SINDy with various upper and lower bounds for the yy-component of heat flux divergence. The bounds that yield the elbow are selected for model discovery.

Algorithm 2 Equation Synthesis

In $\varTheta$, we include each component of the pressure tensor scaled by the thermal velocity, the energy transfer term $\boldsymbol {j}\boldsymbol {\cdot }\boldsymbol {E}$, the Poynting flux terms $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {S}$, each component of the tensor $(\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla })P$ and each term in the scalar product $(P\boldsymbol {\nabla })\boldsymbol {\cdot }\boldsymbol {U}$. These terms were chosen for inclusion in $\varTheta$ as each describes an energy flux and thus has the same physical units as heat flux divergence.

Table 5 summarizes our findings for a representative step (T-17 $\omega _{ci}^{-1}$) during the nonlinear phase and figure 7 demonstrates the temporal validity of the closure.

Table 5. For each component of the heat flux divergence tensor, we give the local approximate Hammett–Perkins closure, the discovered closure and the relevant $L_2$ errors. The $k_0$ for each component of Hammett–Perkins closure is calculated by averaging ${\boldsymbol {\nabla } Q_{ij}}/({v_t(P_{ij} - p\delta _{ij})})$ across the domain.

Figure 7. Snapshots of the discovered closure taken at intervals of 1 $w_{ci}^{-1}$ throughout the nonlinear phase, represented row-wise. The diagonal components of $\boldsymbol {\nabla } Q_{ii}$ are presented column-wise with the associated error in each of the headers. Of note is the consistent performance of the closure as the nonlinear phase progresses.

The analysis was repeated at various times through the nonlinear phase and the showed that the solutions did not vary greatly. Thus, we present the predictions of the closure given by table 5 at evenly spaced time intervals through the nonlinear phase with the associated errors given in the headers of figure 7.

The on-diagonal closures were able to reduce the error significantly when compared to the given closure. The off-diagonal terms demonstrated no such improvement over the given Hammett–Perkins closure. However, our proposed closure comes with some caveats to be discussed further in § 4.

4 Discussion and conclusion

Fully kinetic data may be used to verify the ten-moment multi-fluid model as demonstrated up to and including second-order conservation laws. Specifically, zeroth-, first- and second-order equations were satisfied under visual inspection and with numerical error $L_2$ (left-hand side, right-hand side) $\leq 0.28$ implying the multi-fluid ten-moment model may be validated with kinetic data. This is evidence that the SINDy methodology may be applied to kinetic data, supporting the findings of Alves & Fiuza (2022).

Applying the SINDy methodology, we were able to validate SINDy as a candidate data-driven approach up to and including the on-diagonal second-order terms, as evidenced by the low $L_2$ coefficient error of the discovered solutions. Difficulties were encountered in the discovery of the off-diagonal second-order right-hand side terms (as demonstrated by high coefficient error), where the discovered solutions narrowly deviated from the conservation laws but found more parsimonious solutions with lower numerical error. This is due to the excluded terms having a low order of magnitude leading to elimination by SINDy and not evidence of violation of the conservation laws.

Applied to heat flux divergence, this data-driven approach demonstrates strong results. The existing local approximate Hammett–Perkins closure gives us a baseline to measure improvement. The on-diagonal terms of the heat flux divergence were modelled with significantly reduced numerical $L_2$ error over the baseline. Furthermore, the discovered closures appear to hold through a significant portion of the nonlinear phase, as demonstrated in figure 7. The off-diagonal terms were unsuccessfully modelled and further work is needed to understand this discrepancy.

Other shortcomings of this approach lie in perturbation phase restrictions. While the discovered closure holds well through the nonlinear phase, there is no such guarantees during the linear phase. The local approximate Hammett–Perkins closure was derived using quasi-linear theory, as such, we conjecture that the discovered closure will apply in the linear regime as well. Unfortunately, we cannot demonstrate this at this time, due to the low signal-to-noise ratio in the linear phase.

As a data-driven method, we make no attempt to explain the physics of the discovered closures here. Our sole claim is that the discovered closures will yield improved physical fidelity during the nonlinear phase when compared to the local Hammett–Perkins closure. Model validation is paramount when working with data-driven methodologies, as such future directions in this research will involve applying the discovered closure to a multi-fluid simulation where we hope to see improvements over local Hammett–Perkins closure. Further, we will be looking at expanding the methodology with the goal of pursing a closure which accurately represents the off-diagonal terms.