2.1 Free energy equations

In order to understand the effects and limitations of various closures, it is useful to conceptualize the turbulent dynamics in the context of an energy equation. The free energy (Hatch et al. Reference Hatch, Jenko, Bratanov and Navarro2014) is given by

(2.6)\begin{equation} \varepsilon_{\boldsymbol{k}, n}=\varepsilon_{\boldsymbol{k}}^{(\phi)} \delta_{n, 0}+\varepsilon_{\boldsymbol{k}, n}^{(f)}, \end{equation}

with field component

(2.7)\begin{equation} \varepsilon_{\boldsymbol{k}}^{(\phi)}=\frac{1}{2}\left(\tau+1-\varGamma_{0}\left(k_{{\perp}}^{2}\right)\right)^{{-}1}\left|\phi_{\boldsymbol{k}}\right|^{2}, \end{equation}

and entropy component

(2.8)\begin{equation} \varepsilon_{\boldsymbol{k}, n}^{(f)}=\frac{1}{2} {\rm \pi}^{1 / 2}\left|{f}_{\boldsymbol{k}, n}\right|^{2}. \end{equation}

The free energy evolution equation can be obtained from (2.4) and (2.5)

(2.9)\begin{equation} \frac{\partial \varepsilon_{\boldsymbol{k}, n}^{(\phi)}}{\partial t}={\rm J}_{\boldsymbol{k}}^{(\phi)} \delta_{n, 0}+N_{\boldsymbol{k}, n}^{(\phi)}, \end{equation}

and

(2.10)\begin{equation} \frac{\partial \varepsilon_{\boldsymbol{k}, n}^{(f)}}{\partial t}=\omega_T Q_{\boldsymbol{k}} \delta_{n, 2}-C_{\boldsymbol{k}, n}-{\rm J}_{\boldsymbol{k}}^{(\phi)} \delta_{n, 1} +{\rm J}_{\boldsymbol{k}, n-1 / 2}-{\rm J}_{\boldsymbol{k}, n+1 / 2}+N_{\boldsymbol{k}, n}^{(f)}. \end{equation}

The terms on the right-hand sides of (2.9) and (2.10) represent various energy injection, dissipation and transfer channels. The only energy sink – collisional dissipation $C_{{k}, n}=2 \nu n \varepsilon _{{k}, n}$ – is directly proportional to the Hermite number n multiplied by the free energy. The energy source $\omega _T Q_{{k}}=\omega _T \mathrm {Re}[-({{\rm \pi} ^{1 / 4}}/{2^{1 / 2}}) {\rm i} k_{y} {f}_{2}^{*} \bar {\phi }]$is proportional to the perpendicular heat flux $Q_{\boldsymbol {k}}$.

There are also two conservative energy transfer channels. The nonlinear energy transfer $N_{\boldsymbol {k},n}^{(f)}$ redistributes energy in $k$ space but does not transfer energy between different $n$ and is not a net source or sink (it vanishes under summation in $k$-space).

Here, ${\rm J}_{\boldsymbol {k}}^{(\phi )}=\mathrm {Re}[-{\rm i} k_{z} \phi ^{1 / 4} \bar {\phi }^{*} {f}_{\boldsymbol {k}, 1}]$ is the energy transferred between the field component at $n = 0$ and the entropy component (i.e. Landau damping).

For our purposes of studying closures, the most important terms are the linear phase-mixing terms ${\rm J}_{{k}, n-1 / 2}=\mathrm {Re}[-{\rm \pi} ^{1 / 2} {\rm i} k_{z} \sqrt {n} {f}_{{k}, n}^{*} {f}_{{k}, n-1}]$ and ${\rm J}_{\boldsymbol {k}, n+1 / 2}=\mathrm {Re}[{\rm \pi} ^{1 / 2} {\rm i} k_{z} \sqrt {n+1} {f}_{\boldsymbol {k}, n}^{*} {f}_{\boldsymbol {k}, n+1}]$. These terms also represent a conservative energy transfer channel, albeit in velocity space. They conservatively transfer energy between $n$ and $n - 1$, $n + 1$ respectively but do not transfer energy in k-space. One way to characterize the closure problem is determining the proper value of ${f}_{n+1}$ so that ${\rm J}_{\boldsymbol {k}, n+1 / 2}$ sends the proper amount of energy to higher-order moments – or, as the case may be, receives the proper amount of energy from higher-order moments. Below, in § 4, we will analyse several closures in terms of their capacity to recover the proper (turbulent, kinetic) rates of energy transfer in phase space.

3.1 HP-style closures

Here, we provide a brief description of HP-style closures. Derivations and verification of these closures can be found in appendix A. The HP closure (Hammett & Perkins Reference Hammett and Perkins1990; Smith Reference Smith1997) is designed so that the dispersion relation, also referred to as the kinetic response function, arising from the hierarchy of closed moment equations matches the linear kinetic dispersion relation arising from the Vlasov–Poisson kinetic system. The exact kinetic response function, which involves the plasma dispersion function, $Z(\omega )$, is

(3.1)\begin{equation} R_{00}(\omega) ={-}{\rm i}Z(\omega). \end{equation}

The HP closure for the $N$th moment takes the form $f_N = \sum _{i=0}^{N-1} A_i f_i(\omega )$. Combining this closure ansatz with the hierarchy of moment equations results in an approximate response function $R_{00}^a(\omega )$, a polynomial in $\omega$ involving the closure coefficients, $A_i$.

The HP closure enforces a match in the low frequency limit, $\omega \rightarrow 0$, so the Taylor expansion for plasma dispersion function can be used, which turns (3.1) into a polynomial in $\omega$. The closure coefficients $A_i$ can then be chosen so that $R_{00}^a(\omega ) = R_{00}(\omega )$. A detailed derivation of the HP closure can be found in appendix A.

The HP closure for the 4th moment in our system is

(3.2)\begin{equation} f_{\boldsymbol{k},4} = {\rm sgn}(k_z) A_3 f_{\boldsymbol{k},3} + A_2f_{\boldsymbol{k}, 2}, \end{equation}

where the coefficients are $A_3 = -1.759{\rm i}$ and $A_2 = 0.755$. These coefficients are the same for all $\boldsymbol {k}$, so the only $\boldsymbol {k}$-dependence for this closure comes from the ${\rm sgn}(k_z)$.

In order to test HP-style closures in collisional regimes, we consider a generalization of the HP closure, developed by Snyder in Snyder, Hammett & Dorland (Reference Snyder, Hammett and Dorland1997), which also includes the effects of collisionality. We have developed a collisional extension of the HP closure, the Hammet–Perkins–collisional (HPC) closure, which is inspired by Snyder's method but modified to match the low frequency limit through second order.

The procedure for arriving at the HPC closure for the $N$th moment, $f_{\boldsymbol {k}, N}$, is as follows. First, use the HP method to determine the closure for the $N+1$th moment, then substitute this expression for $f_{\boldsymbol {k}, N+1}$ into the linearized time evolution equation for $f_{\boldsymbol {k}, N}$ and take the low frequency limit of this equation ($\partial f_{\boldsymbol {k},N}/\partial t = 0$). Differentiating this equation with respect to time and then using the low frequency limit of the time evolution equation for the $N-1$th moment yields a collisional closure for the $N$th moment.

The HPC closure for the 4th moment in our system is

(3.3)\begin{equation} f_{\boldsymbol{k},4} = \frac{-3.051{\rm i}k_z\nu - 1.759{\rm i}k_z^2{\rm sgn}(k_z)}{1.838\nu^2 + 3.709 k_z\nu {\rm sgn}(k_z) + k_z^2}f_{\boldsymbol{k},3} + \frac{0.755 k_z^2}{1.838\nu^2 + 3.709 k_z\nu sgn(k_z) + k_z^2}f_{\boldsymbol{k},2}. \end{equation}

This is a second-order accurate (for small $\omega$) closure for $f_{\boldsymbol {k}, 4}$ in terms of $f_{\boldsymbol {k}, 3}$ and $f_{\boldsymbol {k},2}$ including collisional effects. A detailed derivation of this closure and its coefficients can be found in appendix B.

Note that if one takes the collisionless limit, $\nu \rightarrow 0$, of (3.3), the collisionless closure given in (3.2) is recovered.

Both the HPC and HP closures were initially designed for models based on the conventional fluid moments in which the $n$th fluid moment is calculated by integrating the kinetic distribution times velocity to the $n$th power. Subsequent work generalized the procedure for Hermite-based systems (Smith Reference Smith1997). The relationship between the Hermite moments and the fluid moments is very simple and is shown in appendix C.

3.2 The LMM closure

We now ask the question of how a closure may be generalized for a nonlinear system in which the turbulent dynamics continually perturbs the relationships between the low-order moments retained in the system.

This is motivated, in part, by several recent results showing discrepancies between linear and nonlinear phase-mixing dynamics. Plunk (Reference Plunk2013) and Kanekar et al. (Reference Kanekar, Schekochihin, Dorland and Loureiro2015) investigate the effect of a stochastic forcing term on Landau damping rates, demonstrating large deviations from the linear expectations for some parameters. Parker et al. (Reference Parker, Highcock, Schekochihin and Dellar2016), Schekochihin et al. (Reference Schekochihin, Parker, Highcock, Dellar, Dorland and Hammett2016) and Meyrand et al. (Reference Meyrand, Kanekar, Dorland and Schekochihin2019) demonstrate a ‘fluidization’ of collisionless plasma turbulence – i.e. a large reduction of Landau damping rates due to the cancellation of the forward velocity space cascade due to turbulence. Likewise, Hatch et al. (Reference Hatch, Jenko, Navarro, Bratanov, Terry and Pueschel2016b) observes Landau damping rates far smaller than the linear predictions in a turbulent system (see figure 10 of that paper). We thus posit that, in order to capture the phase mixing rates appropriate for a turbulent kinetic system, a closure should be endowed with the versatility to adapt to the nonlinear state. To this end, we propose a closure scheme that learns directly from the turbulent kinetic system.

To illustrate the closure strategy, consider the Hermite-based system described in § 2 at two different truncation levels: (i) a four-moment fluid system, and (ii) a kinetic system of N Hermite moments, where N is large enough that the system is effectively kinetic (in our simulations we opt for $N=48$). For a given wavevector, ${\boldsymbol {k}}$, an eigenvector of the linear operator is simply a vector with the complex values of each moment – i.e. a four-dimensional (4-D) vector in the fluid system and an ND vector in the kinetic system.

The following closure approach is conceptually similar to the HP approach. For a given set of physical parameters (gradient drive, collisionality), solve for the linear eigenvector of the kinetic system. Then use the relationship between $g_4$ and $g_3$ from this kinetic eigenvector to close the fluid system. If the linear eigenvector persists unmodified in the nonlinear state, this approach would be sufficient. However, as described above, important nonlinear modifications are observed in turbulent systems. Consequently, our strategy is to ‘learn’ an appropriate closure directly from the turbulent kinetic system.

We do so by extracting from a nonlinear kinetic simulation an ‘optimal’ basis for the nonlinear turbulent state at each wavevector. For a four-field fluid model, we extract this optimal basis from the first five moments of the kinetic system in order to retain the information necessary to close the system. The turbulent fluid state is then projected onto these basis vectors (with the fifth moment of each removed). Since these basis vectors are attached also to the kinetic information (i.e. the fifth moment), this projection can be used to close the system. Mathematical details are described in the next subsection.

Since an optimal basis was extracted from a kinetic simulation, we would expect this procedure to be effective at the parameter point of the kinetic ‘training’ simulation. The utility of this method, however, will depend on the closure retaining efficacy in some non-negligible parameter domain surrounding the training point. We demonstrate below that this is the case.

We call this closure strategy the LMM closure because (i) it ‘learns’ the closure coefficients from the full turbulent kinetic system, and (ii) it employs multiple modes (basis vectors) in order to better capture the dynamical variations in the turbulent state.

We end this section by noting some connections with other lines of research. First, this closure approach is related to various strategies for projection-based model reduction (Sirovich Reference Sirovich1987; Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993a; Feldmann & Freund Reference Feldmann and Freund1995; Freund Reference Freund2003; Rozza, Huynh & Patera Reference Rozza, Huynh and Patera2008; Peherstorfer & Willcox Reference Peherstorfer and Willcox2016), wherein basis vectors are extracted (often via SVD) from data describing a complex system to reduce the complexity of the underlying models.

We also note some connections with the closure proposed in Sugama, Watanabe & Horton (Reference Sugama, Watanabe and Horton2001) and Sugama, Watanabe & Horton (Reference Sugama, Watanabe and Horton2003). This closure scheme employs two modes (the ITG mode and its complex conjugate) in order to enforce a ‘non-dissipative’ closure – i.e. it eliminates any energy transfer between the fluid moments and higher-order moments. Consequently, it produces damping rates that are far below (i.e. zero) the linear values, qualitatively similar to the nonlinear results cited above. However, the true turbulent system allows energy to shift dynamically between lower- and higher-order moments. Consequently, we view this closure as a compelling idea, but one that is perhaps too restrictive.

We also note the connection between the LMM closure and the line of research exploring the role of damped eigenmodes in plasma microturbulence (Terry, Baver & Gupta Reference Terry, Baver and Gupta2006; Hatch et al. Reference Hatch, Terry, Jenko, Merz and Nevins2011a,Reference Hatch, Terry, Jenko, Merz, Pueschel, Nevins and Wangb, Reference Hatch, Jenko, Navarro, Bratanov, Terry and Pueschel2016b; Whelan, Pueschel & Terry Reference Whelan, Pueschel and Terry2018), which shows that multiple modes co-existing at a single wavevector play a crucial role in turbulent energetics. Our ‘multi-mode’ closure also acknowledges the activity of multiple eigenmodes per wavevector and defines the closure coefficients in terms of the relative amplitude of these modes in the nonlinear state.

3.3 Implementation of the LMM closure

Here, we describe the mathematical details of the approach outlined in the previous section. The closure requires a nonlinear kinetic simulation to formulate a set of basis vectors. In our case, we use 48 Hermite moments for the full kinetic simulation. Any number of subsequent fast fluid simulations can then be run requiring explicit computation of only $f_0$, $f_1$, $f_2$ and $f_3$. In this section we will use bold uppercase letters to denote matrices and bold lowercase letters to refer to vectors.

The full kinetic simulation is used as follows. Let $\boldsymbol {F}_{N\times M}$ ($M$ is the number of time points and $N$ is number of moments retained in the fluid model plus one) be the matrix created from the simulated distribution function at a single wavevector. The distribution function at a single wavevector is written $f_i(t)$, where $i = 0,1,\ldots,N-1$ denotes the Hermite number and $t$ takes on discrete values $t_j$ with $j = 0,1,\ldots,M-1$ (the wavevector is suppressed for clarity), so that element $ij$ of $\boldsymbol {F}$ is $\boldsymbol {F}_{ij}= f_{i}(t_j)$

(3.4)\begin{equation} \boldsymbol{F} = \begin{bmatrix} f_{0}(t_0) & f_{0}(t_1) & \dotsb & f_{0}(t_{M-1}) \\ f_{1}(t_0) & f_{1}(t_1) & \dotsb & f_{1}(t_{M-1}) \\ \vdots & \vdots & \ddots & \vdots \\ f_{N-1}(t_0) & f_{N-1}(t_1) & \dotsb & f_{N-1}(t_{M-1}) \\ \end{bmatrix}. \end{equation}

The SVD of $\boldsymbol {F}$ is given by

(3.5)\begin{equation} \boldsymbol{F}_{N\times M}=\boldsymbol{U}_{N\times N} \boldsymbol{\varSigma}_{N\times N} \boldsymbol{V}^H_{N\times M}, \end{equation}

where $\boldsymbol {U}$ and $\boldsymbol {V}$ are unitary and $\boldsymbol {\varSigma }$ is diagonal with real entries. General background information about this extremely useful matrix decomposition and be found in Golub & Van Loan (Reference Golub and Van Loan2013) and a review on its application to turbulence as proper orthogonal decomposition can be found in Berkooz, Holmes & Lumley (Reference Berkooz, Holmes and Lumley1993b). The columns of the matrix $\boldsymbol {U}$ are called the left singular vectors. In our application, they define $N$ basis vectors for the distribution function. The rows of $\boldsymbol {V}^H$ are the time traces of the amplitude of each of these vectors. The diagonal entries in $\boldsymbol {\varSigma }$ define the singular values, which encompass all the amplitude information. The utility of the SVD lies in its property that the outer product between the first basis vector and the first time trace (weighted by the corresponding singular value) reproduces more of the fluctuation data (as measured by the Frobenius norm) than any other possible decomposition of this form. Likewise the superposition of the first two ($n$) outer products captures more of the fluctuation data than any other rank two ($n$) decomposition and so forth. For convenience, we define a matrix $\boldsymbol {B}$, which weights the basis vectors by their corresponding singular values so that they include the amplitude information: $\boldsymbol {B}=\boldsymbol {U}\boldsymbol {\varSigma }$.

For the purposes of our desired four-moment model, we select $N=5$ (i.e. only a small subset of the 48 total Hermite moments). Since different Hermite moments are only connected to their direct neighbours, this is sufficient to fully exploit the information in the simulation defining the natural (kinetic, turbulent) relations between $f_3$ and $f_4$.

Let $\boldsymbol {f}$ represent the column vector of the first four moments at a single time step: $\boldsymbol {f} = [ f_0\; f_1\; f_2\; f_3]^{\rm T}$. In each time step of a subsequent fluid simulation, we numerically advance $\boldsymbol {f}$ explicitly via (2.4). The truncated moment, $f_4$, is calculated as follows. First, we project the state vector $\boldsymbol {f}$ onto the basis formed by the columns of $\boldsymbol {B}$. This entails finding the projection coefficients that define the amount of each SVD mode in the turbulent state at a given point in time. We will call the column vector containing these projection coefficients $\boldsymbol {c}$. We can do this by removing the row corresponding to the unknown $N$th moment (the 5th row) from $\boldsymbol {B}$ and extracting $\boldsymbol {c}$ from the following equation:

(3.6)\begin{equation} \boldsymbol{f} = \boldsymbol{M} \boldsymbol{c}, \end{equation}

where $\boldsymbol {M}$ denotes the submatrix of $\boldsymbol {B}$ consisting of the first 4 rows and all 5 columns of $\boldsymbol {B}$, i.e. $\boldsymbol {M}$ is the submatrix produced by removing the last (5th) row of $\boldsymbol {B}$. This gives

(3.7)\begin{equation} \boldsymbol{c} = \boldsymbol{M}^{\dagger} \boldsymbol{f}, \end{equation}

where $^{\dagger}$ denotes the pseudo-inverse.

Now that we have $\boldsymbol {c}$, a length $N$ vector of the inferred mode amplitudes, we can predict $f_{4}$ by applying these mode strengths to the previously removed row of $\boldsymbol {B}$, $\boldsymbol {b}_5$. This gives

(3.8)\begin{equation} f_{4} = \boldsymbol{b}_5 \boldsymbol{c} = \boldsymbol{b}_5\boldsymbol{M}^{\dagger} \boldsymbol{f} = \boldsymbol{c}_{{\rm LMM}}\boldsymbol{f}, \end{equation}

where $\boldsymbol {b}_5$ is the 5th row of $\boldsymbol {B}$ and $\boldsymbol {c}_{{\rm LMM}} = \boldsymbol {b}_5 M^{\dagger}$ is the vector containing the 4 LMM closure coefficients. This procedure is repeated at each wavevector $\boldsymbol {k}$ to obtain a full set of coefficients that can be used to conduct an LMM-closed simulation.

An important technical aspect of this procedure is the time grid used for extracting the basis vectors. When generating the matrix $\boldsymbol {F}$, we limit the matrix to the last 70 % of the time domain of the kinetic simulation. This is done to ensure that we extract basis vectors that reflect the turbulent state of the system and not the linear growth phase that occurs at the beginning of the simulation. Each time trace is $\sim 170\,000$ time steps long and when generating $\boldsymbol {F}$, we sample the distribution function every 20 time steps. This means that $\boldsymbol {F}$ contains $\sim 6000$ time points. There is some variation depending on the nonlinear time step, which is adapted to satisfy a Courant–Friedrichs–Lewy (CFL) criterion.

This procedure results in a closure that reflects the natural relations between moments in the turbulent kinetic system and adapts to the relative amplitude of each basis vector in the nonlinear state.

Regarding computational cost, the LMM closure comes down to the dot product between two length 4 vectors: the closure coefficients, $c_{{\rm LMM}}$, and the lower-order moments, $\boldsymbol {f}$. The closure coefficients are computed ahead of time and saved to a file, which is loaded at the beginning of the simulation. During the simulation, the computational expense of the LMM closure is very similar to that of the HP closure; the HP closure requires two complex multiplications per wavevector per time step (one for each of the two HP closure coefficients), and the LMM closure requires four complex multiplications per wavevector per time step. This is much less demanding than the pseudo-spectral computation of the nonlinearity, so the increased expense is negligible. The main additional expense is in running nonlinear kinetic simulations for training. If this can be done sparsely, then the LMM closure is viable.

5.1 Tests of energy dissipation

In order to gain insight into the nonlinear dynamics and its effect on the closure problem, we investigate the energy evolution equation, (2.10). In (2.10), the contribution from phase mixing defines the energy flux to higher-order moments Hatch et al. (Reference Hatch, Jenko, Bratanov and Navarro2014). More specifically, $\chi _{n+1/2} \equiv {\rm J}_{n+1/2}/(|k_z|\,|f_n|^2)$, the normalized rate at which energy is transferred to/from higher-order moments (the phase-mixing rate), is defined by a correlation between two neighbouring moments. The linear physics defines a fixed, dissipative, relationship between $f_n$ and $f_{n+1}$. In the presence of turbulence, however, the various moments are continually perturbed by the nonlinearity, resulting in correlations that can differ substantially from the linear expectation.

These considerations are illustrated in figure 2, which shows the distribution (accumulated over time) of the energy transfer rate between the 3rd and 4th moments, $\chi _{3+1/2}$, for kinetic, LMM-closed, HP-closed and HPC-closed simulations. The average dissipation, $\bar \chi _{3+1/2}$ resulting from the HP and HPC closures is much larger than the dissipation present in the kinetic system. We note that the HP and HPC closures would likely perform better by this metric with the inclusion of more moments, which may be explored in future work. The LMM-closure, being based on the nonlinear system, produces dissipation that matches the kinetic level quite closely. We note that the LMM closure coefficients used are extracted from a training simulation with different parameters than the simulation being examined (both parameter points are noted in the title of the figure), which indicates the effectiveness of the LMM closure even in simulations with parameters different from the training parameters.

Figure 2. Probability distribution functions (a) and box and whisker plots (b) showing the distribution of $\chi _{3+1/2}$, the energy transferred between the 3rd and 4th moments, in the Kinetic, LMM-closed, HPC-closed and HP-closed systems for $\omega _T = 7$ and $\nu = 0.05$ at the most unstable wavevector, $k_x = 0$, $k_y = 0.75$, $k_z = 0.4$. The LMM coefficients used to produce this plot were extracted from the kinetic simulation at $\omega _T, \nu = 6, 0.01$. Red dashed lines show the average value of each distribution, $\bar \chi _{3+1/2}$.

Figure 3 shows the ratio of the average value of $\chi ^{{\rm Closed}}_{3+1/2}$ to the average value of $\chi ^{{\rm Kinetic}}_{3+1/2}$ at the most unstable wavevector for the HP, HPC, and LMM closures for every point in the extended parameter scan. Three different sets of LMM closure coefficients are used to produce this plot: one set is extracted from the kinetic simulation at $\omega _T, \nu = 6, 0.01$, one from $\omega _T, \nu = 9,0.1$ and one from $\omega _T, \nu = 12,0.5$. These three different LMM closures are detailed in § 5. At each point in our parameter grid, we use the LMM closure that was trained closest to that grid point to produce the values of $\bar \chi ^{{\rm LMM}}_{3+1/2}$ shown in this plot.

Figure 3. Ratios of the average value of $\chi _{3+1/2}^{{\rm closed}}/\chi _{3+1/2}^{{\rm Kinetic}}$ at the most unstable wavevector throughout parameter space for the HP, HPC, and LMM closures. Values below 1 indicate not enough dissipation and values above 1 indicate too much dissipation. This heatmap is set up such that fractions far from 1 in either direction are penalized the same way. For example, ratios of 0.5 and 2 will be the same colour. As shown here, the LMM closure matches kinetic dissipation levels much better than the HP and HPC closures throughout most of our parameter grid.

The HP and HPC closures overestimate dissipation levels throughout most of the parameter grid. They perform best at low gradient drive ($\omega _T$) with deteriorating performance as gradient drive is increased. The LMM closure matches kinetic dissipation levels significantly better than both HP closures throughout the parameter space.

It is important to note that the substantial disagreement between kinetic and HP/HPC dissipation levels is somewhat misleading. The discrepancy between the heat flux saturation levels examined in the next section is much less than is suggested by the large difference in phase-mixing rates. While figure 3 clearly shows that the HP/HPC phase-mixing rates to higher $n$ are about a factor of 5 higher than the kinetic result, this difference might not be so important if these phase-mixing rates are still slow compared with nonlinear cascades rates to higher $\boldsymbol {k}$, where hyperdiffusion provides dissipation. We leave more detailed investigation of this possible explanation to future work.

5.2 Comparison of heat fluxes

Ultimately, we would like closed simulations to reproduce the most important macroscopic behaviour of gyrokinetic simulations, notably the saturated value of the radial heat flux, $Q = \sum _{k_x,k_y,k_z} Q_{\boldsymbol {k}}$, where $Q_{\boldsymbol {k}}$ is defined in the discussion surrounding (2.10). We view this metric – the proximity of the saturated heat flux for a given closure scheme to that of the kinetic simulation – to be the most relevant metric for the performance of the closure.

Time traces of the heat fluxes produced in the kinetic simulation and six closed simulations are shown for each combination of input parameters, $\omega _T$ and $\nu$, in figure 4. The final saturation levels of each simulation type at each set of input parameters were calculated by averaging over the last 30 % of the time trace. Each plot in figure 5 shows the per cent error in saturated heat flux, $(Q^{{\rm Closed}}-Q^{{\rm Kinetic}})/Q^{{\rm Kinetic}}$, for all parameter combinations for each closure scheme. Versions of figures 4 and 5 including $\nu =0$ can be found in appendix D. Figure 8 in appendix D contains larger versions of the panels in figure 4 for easier inspection. Comparisons are complicated somewhat by occasional shifts in transport levels that occur unpredictably in time, which introduces a level of uncertainty that cannot be eliminated within the scope of this paper. This is a manifestation of metastable states, recently elucidated in Christen et al. (Reference Christen, Barnes, Hardman and Schekochihin2021).

Figure 4. Time traces of the total radial heat flux ($Q$) for Kinetic (blue), HP (orange), HPC (green), truncated (red), LMM-Middle (purple), LMM-Right (brown) and LMM-Left (pink) simulations for temperature gradient drives ($\omega _T$) ranging from 5 to 15 (increasing downward by panel) and collision frequencies ($\nu$) ranging from 0.01 to 0.5 (increasing to the right by panel). The metric of performance is the final saturation level. The vertical blue lines show the cutoff point $- 70\,\%$ of the simulation time – after which each heat flux curve is averaged to get the final saturation level. Figure 8 in appendix D contains larger versions of the panels in this figure for easier inspection.

Figure 5. Per cent error in saturated heat flux for each closure (HP, HPC, Truncation, LMM-Middle, LMM-Right, LMM-Left, LMM-Optimal) as compared against the kinetic simulation throughout parameter space. The circles on the 4th, 5th and 6th figures indicate which simulations from which the LMM coefficients were extracted. Per cent errors are calculated as $(Q^{{\rm Closed}} - Q^{{\rm Kinetic}})/Q^{{\rm Kinetic}} \times 100$ where $Q^{{\rm closed}}$ is calculated by averaging the last 30 % of the time trace of the heat flux from the closed simulation and $Q^{{\rm Kinetic}}$ is calculated by averaging the last 30 % of the time trace of the heat flux from the kinetic simulation. The r.m.s. error for each closure is also shown above each plot. The 7th plot, LMM-Optimal (3 Points), displays error of the LMM closure trained at the nearest parameter point using three training points. The 8th plot, LMM-Optimal (2 Points), displays the error of the LMM closure trained at the nearest parameter point using only the left and right training points.

As expected, truncation performs poorly at low values of collisionality, but improves at high collisionality where the simulations become more fluid like. Truncation still performs relatively poorly at high collisionality for low gradient drive, but performs well when both collisionality and gradient drive are large – i.e. for parameters at which other effects (gradient drive or collisions) dominate phase mixing.

The HP closure works well in the low collisionality regime for which it was designed (note the small errors at $\nu =0.01$). However, its performance deteriorates as collisionality is increased.

The HPC closure is designed to simultaneously include the effects of collisions and phase mixing in the appropriate limits. As expected, it exhibits notable improvement over the HP closure in the high collisionality regime while also retaining the strong performance of the collisionless HP closure at low collisionality. This closure performs poorly only at intermediate levels of gradient drive ($\omega _T = 7-9$) and in one simulation at $\omega _T,\nu = 15,0.2$. This closure appears to be highly effective and its performance is only surpassed by the LMM closure trained at multiple parameter points, described below.

The LMM-Middle closure based on the kinetic simulation at $\omega _T, \nu = 9, 0.1$ surprisingly performs poorly in the simulation at its training parameter point, likely due to the propensity of this system toward metastable states (note the sudden jump in the LMM time trace toward the end of the simulation). The model does, however, perform well in nearby regions in the middle of our scanned parameter space. In fact, this closure extends throughout parameter space quite well, displaying low errors everywhere except in the top left corner (low collisionality and gradient drive).

The LMM-Right closure based on the kinetic simulation at $\omega _T, \nu = 12,0.5$ produces very low errors at high collisionality and gradient drive, in the middle of the parameter region and even at low gradient drive if collisionality is high. However, it also starts to deteriorate as simulations venture farther from the training point – at low collisionality and gradient drive.

The LMM-Left closure based on the kinetic simulation at $\omega _T, \nu = 6,0.01$ works well at low gradient drive regardless of collisionality. At high values of gradient drive, however, its performance is poor. Specifically at the highest value of gradient drive, $\omega _T = 15$, this closure has extremely large errors. Inspecting figure 4 or 8, one can see the reason for this. The LMM-Left heat flux initially appears to saturate at an accurate level, but part way through the simulation, it jumps to a higher level and re-saturates. This behaviour can also be seen in kinetic simulations at $\omega _T, \nu = 15, 0.2$ and $\omega _T, \nu = 15, 0.5$. When this jump occurs in the closed simulation but not the kinetic, it results in very large errors. These simulations illustrate the potential hazard of applying an LMM closure in a regime too far removed from its training point.

We note several sets of simulations with sudden jumps between saturation levels, likely indicative of metastable states (Heinonen & Diamond Reference Heinonen and Diamond2020; Christen et al. Reference Christen, Barnes, Hardman and Schekochihin2021). This is particularly clear in the lower right-hand corner of the parameter grid, but can be observed throughout the parameter space. It is likely that some fraction of the errors can be attributable to this phenomenon. One example is the LMM-Middle training point, described above. The large HPC error at $\omega _T=15,\nu =0.2$ is another candidate. Unfortunately, this error cannot be eliminated within the scope of this work. We do, however, probe this phenomenon in some more detail.

First, we ran a second set of kinetic simulations with identical physical parameters to the first set. The initial condition is also identical save for a random component, which, in some cases can result in deviations in the long-time simulations. We compared the heat flux saturation levels between the two sets of simulations to quantify the uncertainty in kinetic saturation levels. The per cent difference between the saturation level in each initial simulation was compared with the saturation level of its counterpart with slightly different initial conditions. There was a 5 % r.m.s. difference between the saturation levels in the two sets of simulations. There was one simulation at $\omega _T,\nu =7,0.5$ in which the simulation from the first set experienced a jump from one saturated state to another but the simulation from the second did not. So for this particular parameter point, there was a 20 % difference in the saturation levels.

This illustrates the metastable nature of our system, which introduces some uncertainty into the determination of saturated heat flux levels. That is, because turbulent systems are chaotic, the timing of when a system switches from one metastable state to another, and when it might switch back, can also be chaotic. A reduced closure model of a turbulent system is only expected to reproduce the long-time averages of the exact system. The existence of long-lived metastable states and the possibility that a dynamical system might switch back and forth between them chaotically over time means that, in some cases, longer simulations than usual might be needed for accurate time averages.

As another test, we probed the kinetic simulation at $\omega _T = 15,\nu =0.5$. We changed the collisionality of this simulation from $\nu =0.5$ to $\nu =0.2$ and restarted the simulation with the endpoint of the previous simulation as the initial condition. The results of this test are shown in figure 6. The system finds a saturated state considerably higher than its counterpart with identical physical parameters. The simulation which starts at $\nu =0.5$ initially saturates at $Q\approxeq 840$ but resaturates at $Q\approxeq 760$ once $\nu$ is switched to $\nu =0.2$. The simulation where $\nu =0.2$ the entire time is saturated at $Q\approxeq 460$. It appears that $Q=760$ and $Q=460$ are both stable states for the system with parameters $\omega _T,\nu = 15,0.2$. Which of these metastable states the system falls into depends on the previous state of the system. This test clearly demonstrates that our system has multiple metastable states and exhibits hysteresis. We refer the reader again to a recent, more detailed, study of this phenomenon (Heinonen & Diamond Reference Heinonen and Diamond2020; Christen et al. Reference Christen, Barnes, Hardman and Schekochihin2021).

Figure 6. Time traces of the total radial heat flux ($Q$) from 4 simulations. The blue time trace shows heat flux from the simulation at $\omega _T=15,\nu = 0.5$. At $t\approxeq 1000$, we restarted this simulation from a checkpoint keeping $\nu =0.5$ (orange) and also restarted it from a checkpoint changing the collisionality to $\nu =0.2$ (red). The green line shows the heat flux from a simulation kept at parameters $\omega _T = 15,\nu =0.2$ the entire time. The green and red lines saturate at different levels of $Q$ indicating that there is more than one stable state for the parameters $\omega _T=15,\nu =0.2$ and that which of these states the system falls into depends on the history of the system.

We note that the LMM-Right and LMM-Middle closures are similarly accurate to the HP and HPC closures throughout the parameter space, even taking into account the simulations far removed from the training point. The LMM-left closure exhibits the worst performance of all the closures due to the large late-time jumps in saturation levels at high gradient drive.

Ultimately, we envision the LMM closure applied in scenarios where kinetic training simulations can be supplied sparsely throughout parameter space. Consequently, we show in figure 5 the closure called LMM-Optimal, which is defined by selecting the LMM-closure that is trained at the nearest point in parameter space. For example, at $\omega _T, \nu = 5, 0.05$, the LMM-Left closure with coefficients extracted from the $\omega _T, \nu = 6,0.01$ kinetic simulation is used and at $\omega _T, \nu = 15,0.5$, the LMM-Right closure with coefficients extracted from the $\omega _T, \nu = 12,0.5$ kinetic simulation is used. In figure 5, we include an LMM-Optimal closure that makes use of all three training points and another that uses only the left and right training points. Using two training points instead of three only slightly reduces performance; using three training points results in an r.m.s. error of 8 % and using two yields an r.m.s. error of 12 %. Tables 1 and 2 in appendix D show which training simulation was used to run the LMM-closed simulation at each grid point for the LMM-Optimal (3 Point) and LMM-Optimal (2 Point) plots in figure 5. Choosing the coefficients this way – using the LMM coefficients extracted from the nearest kinetic simulation to each grid point – results in excellent performance. Only one parameter point exhibits an error above 20 % and the r.m.s. error is only 8 % when three training points are used.

The utility of the LMM closure would depend on efficient exploration of parameter space. A database of closure information could be maintained as a parameter space is explored. Expensive kinetic simulations could be performed sparsely, guided by a rigorous statistical framework like Bayesian optimization. Fast fluid simulations, could then select the closure information at the nearest available parameter point.