2.1 Limitations of the current fluid theory

The fluid theory was developed for the ideal fluid (i.e. related to a MDF) by Euler in the 18th century, and more significantly for non-ideal flows in the 19th century with the well-known works of Navier, Cauchy, Poisson, Stokes, Reynolds and others who found fluid equations describing transport and turbulence. In parallel, the kinetic theory of gas was developed and allowed for a deeper description and understandings of transport and turbulence. However, the fluid theory has not been able to recover some results of the kinetic theory yet. Indeed, the first limitation which has been recurrent in the last century is the use of non-adapted basis functions for the representation of the distribution function. One of the first noticeable works linking non-Maxwellian deviations of the distribution function and the fluid moments was published by Grad (Reference Grad1949, Reference Grad1963). In these references, Grad used the Hermite polynomials and obtained the well-known 13-moments fluid model. Since this result was obtained, all analytic representations of the distribution function have used one of the following orthogonal basis sets: Hermite, Laguerre or Legendre polynomials, the Bessel functions or the Fourier series (see figure 1).

Figure 1. Orthogonal basis functions commonly used to represent NMDFs. The limit at infinity of Hermite (a), Laguerre (b) and Fourier (c,d) basis functions are not zero. The Bessel basis functions (e) are quasi-periodic and the Legendre basis functions (f) are relevant only for discretized distribution functions. None of them are efficient (e.g. a MDF can be described with at least a few hundreds or thousands of terms depending on the required level of accuracy).

These representations are not efficient for the description of non-Maxwellian steady-state distribution functions or even for the fluctuation part. Indeed, for the Hermite and Laguerre polynomials, as well as for the Fourier series, their limit at infinity is not 0 so these representations are not adapted for localized modifications of the distribution function. Moreover, due to the fact that these basis functions are not compact, a very large number of terms are required to describe a compact deviation of the distribution function. The multiplication of these basis functions by a background Maxwellian (e.g. Grad Reference Grad1949) is still not efficient enough even if the distribution function becomes compact because the departure of the distribution function from the Maxwellian (i.e. $f/f_{0}$ where $f$ is the NMDF and $f_{0}$ the MDF) is also not well described by all of these orthogonal basis functions (Izacard Reference Izacard2016a ). For the Bessel functions and the Legendre polynomials, their limit goes to zero but the quasi-periodicity of the Bessel functions and the non-efficiency of the Legendre polynomials are additional reasons for the failure toward an efficient description of NMDFs. Indeed, the Legendre polynomials are adapted for a discretized distribution function which is not efficient for the fluid reduction due to the large number of terms (hence, the large number of fluid equations). Furthermore, even for the description of a MDF, one needs a few hundreds or thousands of terms with a discretized distribution function instead of only 3 parameters with the continuous exponential function. In the end, the fluid reduction can be efficient only by using continuous, non-quasi-periodic and compact basis functions which is the case for the MDF. An important related limitation found in Hirvijoki et al. (Reference Hirvijoki, Candy, Belli and Embréus2015) is the use of a sum of shifted MDFs (called in the literature the radial Gaussian basis functions, RGBF)Footnote ¹ . This use is a limitation because for an arbitrary collisionality, a MDF is the steady-state solution of an isolated system (i.e. when the self-collisional time is much smaller than the characteristic time of exchange of energy with external sources). Then, it is rarely self-consistent to superpose at least two MDFs corresponding respectively to the bulk plasma and to the super-thermal population. Both reach their thermodynamic equilibrium by neglecting all interactions (i.e. collisions) with other species or injected particle sources. However, the presence of the super-thermal population is due to non-negligible collisions with a source of energy (e.g. neutral beam injection (NBI) or $\unicode[STIX]{x1D6FC}$ -particles). The last sentence can enter into conflict with the fact that both species are at their thermodynamic equilibrium, which is reached only when we neglect the collisions with other species. In summary, the sum of MDFs and more generally the RGBF are valid only in the specific regime of collisions

where $a\in \{\text{th},\text{f}\}$ , $b\in \{\text{th},\text{f},\text{s}\}$ and $a\neq b$ ( $\text{th}$ , $\text{f}$ and $\text{s}$ represent respectively the thermal, fast and source populations). In other words, $\unicode[STIX]{x1D708}_{a-a}$ represents the self-collisionality of the thermalized or the fast population and $\unicode[STIX]{x1D708}_{a-b}$ represents the collisionality of the interaction between thermalized population, fast particles or sources. This regime is too constraining for many experiments where a finite ratio needs to be taken into account.

Figure 2. Example of basis functions relevant to INMDFs. These basis functions efficiently describe localized modifications of NMDFs and parameters can be added ‘à la Maxwellian’ to allow robust representations.

Figure 3. (a) represents the comparison of the fitting of artificial data (black crosses) with a Maxwellian (blue curve), a Hermite polynomials of order $N=24$ (red curve) and an INMDF (green curve). The artificial data are constructed from an INMDF and 5 % of noise. (b) shows the relative error in per cent between the fitted distribution function and the artificial data as a function of the number of terms, $N$ , used by the representation of the distribution function.

Finally, one of the arguments found in the literature for the use of the representations mentioned above is that their basis functions are orthogonal. However, this argument is only needed for the conveniences of analytic developments for a very large number of terms and is not physically motivated. We use here the generalized formula of the interpreted NMDF (called INMDF) introduced by Izacard (Reference Izacard2016a ,Reference Izacard b )

with $m_{k}=2\text{E}[(n_{k}+1)/2]+1$ , $n_{k}\in \mathbb{N}$ where $\text{E}[x]$ is the floor function and $a_{k}$ , $b_{k}$ , $c_{k}$ , $d_{k}$ , $e_{k}$ are fluid coefficients for all integer $k$ ( $N$ is the maximum number of terms) and $\text{v}$ is the velocity phase-space coordinateFootnote ² . The set of coefficients $\{a_{k},b_{k},c_{k},d_{k},e_{k}\}$ are the fluid hidden variables. The solution shown in figure 2 solves the problem of the non-adapted orthogonal basis functions used in the literature. This figure represents the first terms of possible basis functions that can efficiently describe NMDFs, especially with a super-thermal tail. We remark that these basis functions are not orthogonal in comparison to the previously mentioned ones. Moreover, this set of basis functions is so general that it can reproduce a sum of MDFs as well as the Hermite polynomials (O. Izacard, Private communications with J. Candy, 2013–2014). The use of the first basis function $f_{0}(\text{v})$ allows us to describe MDFs with 3 hidden variables: the density $n$ , the fluid velocity $v$ and the temperature $T$ . Each of these 3 parameters has a well-known physical interpretation. Following this methodology, the INMDF (see (Izacard Reference Izacard2013, Reference Izacard2016a ,Reference Izacard b ) and (3.3)) is the first example of a NMDF representing a super-thermal tail that uses the second basis function $f_{1}(\text{v})$ associated with $3$ additional hidden variables: the kinetic flux $\unicode[STIX]{x1D6E4}$ , the central flow $c$ and the width of the heat spread $W$ . These new parameters have been physically interpreted in Izacard (Reference Izacard2016b ) as a displacement of a population of particles in the velocity phase space due to an external source of energy. The efficiency of the INMDF to describe a localized super-thermal tail is shown in figure 3 which compares different orders of the Hermite polynomials and the INMDF to fit artificial data. The artificial data are created from an INMDF which includes 4 % of super-thermal particles and 5 % of noise. The crosses in figure 3(a) are a plot of the artificial data. The blue, red and green curves correspond respectively to the least-squares fitting of the artificial data by a MDF, an INMDF or a Hermite polynomial of order $N=24$ . We remark that due to the limited number of velocity grid points (i.e. 101 points), the best fit from a Hermite polynomial is obtained with 24 terms (see figure 3 b) and if we want to reduce the least-squares errors we have to increase the velocity phase-space resolution. The fitted MDF (blue curve) is defined by the three first moments of the artificial data. Figure 3(b) represents the relative error in per cent between the artificial data and the fitted curves as a function of the number of used parameters. The Maxwellian uses $N=3$ parameters, the INMDF uses $N=6$ and the Hermite polynomial gives a minimal error (for the artificial data shown in figure 3 a) when 24 coefficients are kept. The error of the Hermite polynomial is due to the fluctuations at velocities far from the fluid velocity $v$ (i.e. far away from $\text{v}/\text{v}_{max}\sim 0.5$ ) as well as the mismatch of the distribution function close to the central flow $c$ ( $\text{v}/\text{v}_{max}\sim 0.6$ ). In the Hermite representation, a large number of terms are required in order to reduce the divergence at infinity induced by the first added term. In other words, the Hermite polynomial is not efficient for the description of local departures of a MDF, in the same way that a polynomial will not be efficient in general for the representation of a cosine. The INMDF is clearly much more accurate and efficient in reproducing super-thermal particles with a minimum number of parameters (i.e. hidden variables) even when the super-thermal tail represents more than 4 % of particles as shown in figure 3(a).

Summarizing the first limitation given here, the choice of the basis function used for the description of NMDFs is crucial since the number of hidden variables (i.e. the fluid parameters) is completely related to this choice. The INMDFs generalized by (2.2) seem to be one of the best ways to describe localized deviations from the MDF, at least in the presence of super-thermal particles. The reader may develop other non-orthogonal basis sets in order to describe other phenomena (see other non-orthogonal basis sets used for neoclassical particles in Izacard Reference Izacard2016a ).

A second limitation found in the literature is the acceptance of misunderstandings about the fluid theory and the fluid closure. Here is a list of some of these misunderstandings accepted by the majority of the community.

1. (1) F-(i) Many researchersFootnote ³ agree about the fact that the fluid codes are useless with respect to (gyro)kinetic codes for the description of many phenomena in plasmas since they cannot capture kinetic effects. We must agree that, until now, no rigorous and efficient developments of fluid equations catching kinetic effects, even in the collisionless limit, have been published. Indeed, this work is the proof that the fluid theory can efficiently describe kinetic effects. The efficiency of the fluid reduction is obviously dependent on the choice of the basis functions used to represent the distribution function in the velocity phase space.

2. (2) F-(ii) Another misunderstanding about the fluid theory is that the fluid closure appears only due to the fluid reduction of the kinetic Boltzmann equation. This statement is incorrect since the fluid closure appears due to the representation of the distribution function by a finite number of terms. It can be either by the discretization of the distribution function in the velocity phase space, or by the truncation of the number of basis functions used to approximate a continuous distribution function. This means that all existing (gyro)kinetic codes are related to a fluid closureFootnote ⁴ .

3. (3) F-(iii) By performing the fluid reduction of the kinetic equation without using a truncation of the representation of the distribution function in the velocity phase space, we obtain the hierarchical fluid equations for an infinite number of fluid moments. By assuming the truncation of the dynamical evolution of a finite number of fluid moments, we obtain an additional unknown fluid moment in the last equation (i.e. $N+1$ unknowns for $N$ equations). This is well known as the problem of the fluid closure and many studies have used a relation of this additional fluid moment as a function of the previous fluid moments. We can categorize the different fluid closures found in the literature. The first category is the local collisional fluid closures such as found in Braginskii (Reference Braginskii1965). The second is the non-local collisional fluid closure such as found in Hammett & Perkins (Reference Hammett and Perkins1990), Dimits, Joseph & Umansky (Reference Dimits, Joseph and Umansky2014). The third one is the local collisionless fluid closure such as found in Grad (Reference Grad1949, Reference Grad1963) where the distribution function is represented by a truncation of the Hermite polynomials. The latter gives the well-known 13-moments fluid model but as explained above, this representation is not efficient for the description of localized deviations of NMDFs. Moreover, another collisionless fluid closure which is not physically correct is found in De Guillebon & Chandre (Reference De Guillebon and Chandre2012), where the authors assume that $M_{k}=0$ for $k>N$ . Their choice is too constraining for non-Maxwellians because even the MDF cannot be described by this method. Indeed, for a MDF, only the even moments $\{P_{2k+1},k\in \mathbb{N}\}$ are cancelled and the odd moments $\{P_{2k},k\in \mathbb{N}\}$ are functions of the density, fluid velocity and temperature, and none of the moments $M_{k}$ are equal to $0$ (the moments $M_{k}$ and $P_{k}$ are defined in appendix A). In the case of the INMDFs given by (2.2), the collisionless and collisional local fluid closures are obtained by writing the additional moment as function of the hidden variables ( $M_{k}=A(u)$ , where $A$ is an operator that act on the fluid hidden variables listed in $u$ ) instead of a function of the previous moments. In fact, there is no physically motivated reason to have reversible relations between the hidden variables ( $u=\{(a_{k},b_{k},c_{k},d_{k},e_{k}),\text{for }k\in \mathbb{N}\}$ ) and the fluid moments (i.e. such as for the MDF). This non-reversibility is one of the reasons why the collisionless fluid closure of NMDFs has never been developed until now. However, it is important to remark that even if this relation is not reversible, the next generation of fluid code (based on the generalized equations given in § 3) can evolve the dynamical equations of the hidden variables $u$ by inverting the matrix $\unicode[STIX]{x1D63D}$ where $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}M_{k}=\unicode[STIX]{x1D63D}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}u$ for $\unicode[STIX]{x1D709}$ a space or time coordinate (i.e. $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}u=\unicode[STIX]{x1D63D}^{-1}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}M_{k}$ ).

4. (4) F-(iv) It has been argued that the fluid theory is valid only for negligible mean-free-path lengths in comparison to all characteristic scale lengths. This means that the fluid theory is only valid when the collision frequency is higher than all characteristic frequencies, as explained by the quoted citation above. Hence, one may think that the fluid theory is more valid when the collisionality goes to infinity. However, many experiments are not consistent with this infinite collisionality assumption and the fluid theory can still produce a good representation if we use a consistent set of fluid equations (i.e. at least with our collisionless fluid closure shown in § 3). Indeed, for an infinite collisionality the distribution function tends to be closer to a MDF but for finite collisionality NMDFs may be observed. The fluid theory is not invalid for a finite collisionality (i.e. departures of Maxwellian distributions when the mean free path is not negligible) but the validity of the fluid theory is assured by a large number of particles (in order to have a good statistical kinetic representation), and is particularly related to the representation of the distribution function with continuous basis functions in the velocity phase space. In other words, the generalized fluid theory is valid when the kinetic theory is.

As a conclusion of the limitations of the fluid theory found in the literature, we argue that the description of fluid models including non-Maxwellian kinetic effects even in the collisionless limit has never been obtained self-consistently since 1872 and 1877 when Maxwell (Reference Maxwell1872) and Boltzmann (Reference Boltzmann1877) proved that only the MDF is the solution of the Boltzmann equation without external sources. This makes sense because all attempts to represent variations in the velocity phase space of NMDFs were (i) not efficient, reducing to more than a few hundreds or thousands of fluid moments and (ii) not self-consistent with an intermediate collisionality. The second reason of this failure is the misconception of the fluid reduction from the kinetic Boltzmann description, and more generally of the fluid theory: it has been difficult to expand this fluid theory here for non-equilibria based on the properties of MDFs without falling into the usual limitations, but a summarized point of view of the state of the art of the fluid reduction since the 1870s clarifies our novel approach.

The introduced INMDFs represent an extension of the usual fluid theory since (i) it has been shown (Izacard Reference Izacard2016b ) to fit a numerical NMDF which resolves the TS-ECE discrepancy of the electron temperature measurement in JET (De La Luna et al. Reference De La Luna, Krivenski, Giruzzi, Gowers, Prentice, Travere and Zerbini2003; Beausang et al. Reference Beausang, Prunty, Scannell, Beurskens, Walsh and De La Luna2011), by using only 6 additional parameters, (ii) the collisionless and collisional fluid closure shown in § 3 prove that the fluid closure is entirely related to the choice of the description of the distribution function in the velocity phase space and (iii) the INMDF can be consistent with a new solution of the Boltzmann equation for intermediate collisionality regimes in the presence of external sources.

The INMDF seems to prove that all accepted misunderstandings about the fluid theory mentioned previously can be bypassed.

2.2 Inconsistency between all fluid models and non-Maxwellians

We have previously argued that all existing fluid models are related to a MDF. This part explains the reasons of this strong assumption. Starting from the Boltzmann equation and an assumed non-Maxwellian steady state such as the one measured in experiments or computed from PIC or Fokker–Planck codes, we know that if the NMDF is smooth enough (i.e. continuous and infinitely derivable) then it is a solution of the collisionless limit of the Boltzmann equation (i.e. the Vlasov (Reference Vlasov1938) equation $\unicode[STIX]{x2202}_{t}\,f+D[\cdots \,]=0$ where $D[\cdots \,]$ is the advection operator described in appendix A for the fluid equations). Hence, the collisionless steady-state limit of all fluid models should give us some information about the link between these fluid models and their associated NMDFs.

In a general point of view, the fluid models can be written in the form

where $D_{\unicode[STIX]{x1D709}}$ (respectively $S_{\unicode[STIX]{x1D709}}$ ) are the dissipative (respectively source) ad hoc terms for $\unicode[STIX]{x1D709}\in \{n,v,T\}$ . The fluid hierarchy is such that the equation of the temperature involves the next fluid moment, i.e. the heat flux $q$ . Usually, the fluid closure assumes a heat conduction relation

where $\unicode[STIX]{x1D712}$ is the constant heat conductivity coefficient (other fluid closures are discussed in § 2.3 where the local spatial gradient is replaced by a Hilbert or a Fourier transform). The ad hoc dissipative and source terms approximate kinetic effects related to the fluctuations of the NMDF around its Maxwellian steady state (i.e. it is a perturbation theory). This technique using ad hoc transport coefficients is not self-consistent because the kinetic effects (approximated by the transport coefficients, see § 3) are static at the same time that standard fluid quantities (i.e. $n$ , $v$ and $T$ ) are dynamical. Our generalized fluid theory shown in § 3 is self-consistent because the kinetic effects and the associated transport are also dynamically evolving. However, the steady-state limit of these fluid models given by (2.3)–(2.5) becomes

with the fluid closure given by

This fluid closure is the well-known adiabatic (or double adiabatic, when pressure anisotropy is retained) approximation (Chew, Goldberger & Low Reference Chew, Goldberger and Low1956), and its validity conditions are well known as well, especially at the level of the linear kinetic theory. All ad hoc dissipative and source terms are dropped in steady state because the distribution function is assumed to be Maxwellian. Indeed, it turns out that the steady-state limits of almost all existing fluid models found in the literature are only related to MDFs and are not consistent with the description of non-Maxwellian effects. This means that all non-Maxwellian steady-state distribution functions measured in experiments or computed from PIC or Fokker–Planck codes cannot be described self-consistently and rigorously by any of the existing fluid models. We remark that the most advanced technique helping to better mimic kinetic effects in the state-of-the-art existing fluid codes consist of using radial and poloidal profiles of the transport coefficients $D$ and $\unicode[STIX]{x1D712}$ such as in Canik et al. (Reference Canik, Maingi, Soukhanovskii, Bell, Kugel, Leblanc and Osborne2011), Groth et al. (Reference Groth, Brezinsek, Belo, Beurskens, Brix, Clever, Coenen, Corrigan, Eich and Flanagan2013), Meier et al. (Reference Meier, Soukhanovskii, Bell, Diallo, Kaita, Leblanc, Mclean, Podesta, Rognlien and Scotti2014). This technique has been developed to match midplane profiles of density and temperature obtained from experimental measurements in expectation of better describing the dynamics of the scrape-off layer (SOL) and the divertor plasma. It has been found that this profile fitting technique reduces the radiation shortfall (Groth et al. Reference Groth, Brezinsek, Belo, Beurskens, Brix, Clever, Coenen, Corrigan, Eich and Flanagan2013) and can help to simulate the transport enhancement between the two X-points of a snowflake magnetic configuration (Rognlien Reference Rognlien2014; Izacard et al. Reference Izacard, Scotti, Soukhanovskii, Rensink, Rognlien and Umansky2016c ). However, even if this technique allows us to describe different transport levels at different positions, these transport coefficients are assumed ad hoc and there is no direct relation with the distribution function. As described below, another extension to these fluid models is the use of non-local transport coefficients.

2.3 Non-self-consistent solution: non-local heat transport and others

One of the most advanced attempt to approximate some kinetic effects in the fluid equations is to use the non-local heat transport theory. In some cases, instead of using the usual fluid closure $q=-\unicode[STIX]{x1D712}\unicode[STIX]{x1D735}T$ where $\unicode[STIX]{x1D712}$ is an ad hoc coefficient constant in time and fixed in order to match some experimental data, only the use of a non-local definition of the dissipative coefficients can enhance the validation of simulations against experiments (Dimits et al. Reference Dimits, Joseph and Umansky2014) and can be consistent with some kinetic effects (Hammett & Perkins Reference Hammett and Perkins1990; Chang & Callen Reference Chang and Callen1992; Snyder, Hammett & Dorland Reference Snyder, Hammett and Dorland1997). It turns out that this choice uses the linearization of the Landau (Landau Reference Landau1946; Van Kampen Reference Van Kampen1955; Case Reference Case1959; Bacus Reference Bacus1960; Ong Reference Ong1963; Bernstein, Trehan & Weenink Reference Bernstein, Trehan and Weenink1964) or Fokker–Planck equations (similar to Ji & Held Reference Ji and Held2013), or a Hilbert or Fourier transform (Sedlàček Reference Sedlàček1967), and is always associated with fluctuations of the distribution function. A non-local heat flux is defined by $q=\int F[n,v,T]\,\text{d}^{3}x$ where $F$ is an ad hoc functional of the fluid quantities (usually motivated by the description of the phenomenon under investigation). This non-local transport theory can be robust because the non-local definition efficiently adds a time evolution of the ad hoc dissipative transport coefficients as a function of the fluid quantities, but there is no robust and simple link between a specific shape of the NMDF and the ad hoc functional $F$ . Similarly, more advanced fluid closures at higher-order ranks have been proposed (Snyder et al. Reference Snyder, Hammett and Dorland1997; Scott Reference Scott2010; Sulem & Passot Reference Sulem and Passot2014) but they always involve perturbations of the distribution function. The fluid closure of Grasso & Tassi (Reference Grasso and Tassi2015) is obtained by adding strong constraints from Scott (Reference Scott2010) due to the interest of the authors in finding the Hamiltonian structureFootnote ⁵ . Even though the non-local transport is a way to describe modifications of the anomalous transport due to turbulence and possibly non-Maxwellian kinetic effects, this non-local theory is still a perturbation theory and is not consistent with any steady-state NMDF which should naturally include a local closure as shown in § 3.

3.1 Local collisionless fluid closures

The fluid equations derived from the moments of the Vlasov equation (i.e. the collisionless limit of the Boltzmann equation), are written in a highly compact form using the moments $M_{k}$ rather than the moments $P_{k}$ (see definitions in appendix A). The fluid equations (previously noted $\unicode[STIX]{x2202}_{t}M_{k}=D[\cdots \,]$ ) read

for $k\in \mathbb{N}^{\star }$ . By assuming that the distribution function is described by the INMDF given by (3.3), all fluid moments $M_{k}$ are functions of the $6$ hidden variables $(n,v,T,\unicode[STIX]{x1D6E4},c,W)$ (see appendix B), especially the moment $M_{6}$ which appears in the last dynamical moment equation $\unicode[STIX]{x2202}_{t}M_{5}$ . The collisionless fluid closure associated with the INMDF (Izacard Reference Izacard2016a ,Reference Izacard b ) that is given by

where $\unicode[STIX]{x1D6E5}=n(2\unicode[STIX]{x03C0}T)^{-1/2}\exp (-v^{2}/(2T))$ and $\unicode[STIX]{x1D6E5}_{I}=\unicode[STIX]{x1D6E4}(2\unicode[STIX]{x03C0}W^{3})^{-1/2}\exp (-c^{2}/(2W))$ , reads

The physical interpretation of the hidden variables $(n,v,T,\unicode[STIX]{x1D6E4},c,W)$ is given by Izacard (Reference Izacard2016b ). From the point of view of numerical codes, we have to remark that an algorithm for each time step can be constructed as the following:

1. (i) A-(1) Each time step of the generalized fluid code needs inputs. For each time step, except the initial condition, the fluid hidden variables $u=(n,v,T,\unicode[STIX]{x1D6E4},c,W)$ associated with the INMDF are given. Only the initial time step may be different: we can use a numerical distribution function $f_{init}$ or the numerical values of its $6$ first fluid moments $M=(M_{0},M_{1},M_{2},M_{3},M_{4},M_{5})$ . In both cases, a fitting algorithm is required in order to match $f_{init}$ or $M$ with the analytic computation assuming the INMDF and the fitted initial conditions of the hidden variables contained in $u$ .

2. (ii) A-(2) The second step of this algorithm is to numerically compute the right-hand side of the dynamical equation of the fluid moments $M$ (see appendix A). For that, we need to compute the fluid moments $M$ as function of the hidden variables $u$ (see appendix B) $M=A(u)$ where $A$ is a polynomial operator that acts on the hidden variables  $u$ , as well as the spatial derivatives of the fluid moments $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}M_{k}=\unicode[STIX]{x1D63D}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}u$ where $\unicode[STIX]{x1D709}\in \{x,y,z\}$ and $\unicode[STIX]{x1D63D}$ is a simple matrix that is a function of the hidden variables $u$ . In other words, we compute the right-hand side of the dynamical equations of the fluid moments $\unicode[STIX]{x2202}_{t}M$ as functions of the hidden variables $u$ and their spatial derivatives.

3. (iii) A-(3) The third step of this algorithm is to invert the matrix $\unicode[STIX]{x1D63D}$ (a function of $u$ ) which also links the time derivatives $\unicode[STIX]{x2202}_{t}M=\unicode[STIX]{x1D63D}\unicode[STIX]{x2202}_{t}u$ . Then, $\unicode[STIX]{x2202}_{t}u=\unicode[STIX]{x1D63D}^{-1}\unicode[STIX]{x2202}_{t}M_{k}$ where $\unicode[STIX]{x1D63D}^{-1}$ is a simple inverse matrix that is a function of the hidden variables $u$ .

4. (iv) A-(4) The last step is the computation of the time evolution of the hidden variables $u$ using one of the available numerical schemes. Then, the $6$ fluid hidden variables $u$ are updated and the loop of the algorithm is closed.

An important remark is the fact that the closure on the moment $M_{6}$ is part of the step A-(2) in this algorithm because the fluid moment $M_{6}$ is known as a function of the hidden variables $u$ .

At this point, we have introduced here $6$ fluid equations associated with the INMDF whereas the current fluid theory is mainly based on the common $3$ fluid equations associated with MDFs (see § 2.3 for discussion on exceptions). However, it is possible to reduce the degrees of freedom of the INMDFs by assuming $2$ constraints relevant to physical motivations of the system under consideration. As an example, it is possible to use $2$ constraints if one wants to describe INMDFs where the central flow $c(n,v,T,\unicode[STIX]{x1D6E4})$ and the width of the heat spread $W(n,v,T,\unicode[STIX]{x1D6E4})$ are given as functions of the $4$ dynamically independent hidden variables $u=(n,v,T,\unicode[STIX]{x1D6E4})$ . As an example from the construction of the coefficients $r$ and $s$ in Izacard (Reference Izacard2016b ) where $c=v+r\sqrt{T}$ and $W=s^{2}T$ , we can constrain a specific value for the coefficients $r$ and $s$ in order to remove the independence of $c$ and $W$ from the $4$ other fluid hidden variables. In this case, the collisionless fluid closure that appears in the dynamical equation $\unicode[STIX]{x2202}_{t}M_{3}$ is given by

where the constraints on $c(n,v,T,\unicode[STIX]{x1D6E4})$ and $W(n,v,T,\unicode[STIX]{x1D6E4})$ are required here.

These two 4-moments and 6-moments fluid models are the first ones efficiently related to a NMDF that can describe a tail or a dissymmetry. Both models are directly reduced from the INMDF. A large number of reduced models can be obtained from this set of fluid equations and the perturbative theory may help in adding more terms related to turbulence, in the same way it has been developed in the previous decades for fluid models close to MDFs. Indeed, we can construct fluid models for turbulence (i.e. fluctuations of the distribution function) close to an INMDF.

Figure 4. Validation studies of (3.6) and (3.7) proving the kinetic origin of the particle diffusion via the competition between the growth of the INMDF given by (3.3) and the thermalization relaxation simulated by forcing $\unicode[STIX]{x1D6E4}=0$ at specific times. Based on this minimal INMDF fluid model, the collisionless simulation (a) shows propagation of blob structures and highly collisional simulation (b) shows particle diffusion.

As shown for example for the characteristic curve of the Langmuir probe by Izacard (Reference Izacard2016a ,Reference Izacard b ) where the INMDF was able to replace an ad hoc diffusion term, a population of super-thermal particles can naturally enhance the transport. It is then trivial to model fluctuations by a sum of INMDF deformations of the distribution function. We might be able to recover the anomalous transport due to the presence of turbulence by only considering NMDFs. Moreover, the dissipative terms present in the turbulence models can be recovered by specific constraints of the hidden variables. For example, it turns out that if we impose the kinetic flux $\unicode[STIX]{x1D6E4}=-D\unicode[STIX]{x1D735}n$ , the divergence of the first moment $M_{1}$ which appears in the continuity equation $\unicode[STIX]{x2202}_{t}n$ becomes $-\unicode[STIX]{x1D735}\boldsymbol{\cdot }M_{1}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(nv)+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(D\unicode[STIX]{x1D735}n)$ . Then, if the diffusion coefficient is homogeneous, we recover the usual particle diffusion $-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6E4}=D\unicode[STIX]{x1D735}^{2}n$ . Another interesting observation is found when we assume the constraints $v=c=0$ and constants $T$ and $W$ . We then obtain from the dynamical moments equations $\unicode[STIX]{x2202}_{t}M_{k}$ for $k\in [0,1]$ (see the moments in appendix B) of an INMDF the two following equations

Then, with the linearization of the term $\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D6E4}$ , we find the relation $D=T/\unicode[STIX]{x1D6FE}$ where $\unicode[STIX]{x1D6FE}$ is the linear growth rate of the kinetic heat flux $\unicode[STIX]{x1D6E4}$ . In other words, the diffusion coefficient $D$ has its origin in a static linear growth rate of the kinetic heat flux. As a validation test, figure 4 shows two numerical simulations of the reduced INMDF fluid model given by (3.6) and (3.7). The initial value of the density $n(x,it=0)$ is a Gaussian function (solid black curve) and that of the kinetic heat flux $\unicode[STIX]{x1D6E4}(x,it=0)$ is $0$ (dashed black curve). Moreover, the non-periodic spatial gradients are evaluated with the numpy.gradient function in Python and the time integration with a first-order difference scheme. If we let the system $(n,\unicode[STIX]{x1D6E4})$ evolve in time, we observe in figure 4(a) the splitting of the Gaussian of the density in two propagating Gaussians in opposite directions (blue, red and green curves) and the propagating velocity is directly a function of the constant temperature $T=0.5$ (we also fix $W=T/2$ for the INMDF). However, after $600$ iterations in time (at $t=it\times \text{d}t$ , where here $it=600$ and the time step $\text{d}t=20$ in arbitrary units) we chose to manually reset $\unicode[STIX]{x1D6E4}(x,it=600)$ to $0$ because in this reduced INMDF fluid model there are no collisions that assure a positive distribution function. This manual reset simulates the thermalization of the INMDF toward the MDF (when $\unicode[STIX]{x1D6E4}=0$ ). We note the effect of resetting $\unicode[STIX]{x1D6E4}(x,it=600)=0$ (green curves) in figure 4(a): each Gaussian splits in 2 new Gaussians propagating again in opposite directions (yellow, cyan and magenta curves). However, it turns out that due to the high value of the time step and spatial gradients, as well as our arbitrary choice of the temperature $T$ and width of heat spread $W$ , the distribution function becomes negative after a few time steps and its minimum value saturates when the dynamics becomes purely advective. We remark that other choices of $T$ and $W$ would have resulted in a saturated advection with a positive INMDF but the dynamics would have been slower. Figure 4(b) is the same simulation but the kinetic heat flux $\unicode[STIX]{x1D6E4}$ is reset to $0$ at every time step before getting a negative distribution function (see criteria in Izacard Reference Izacard2016b ). The finite collisionality is simulated by the duration $\unicode[STIX]{x0394}it$ of the thermalization which varies from $6$ (at $it=0$ ) to $181$ (at $it=12\,000$ ) time steps due to the reduction of the spatial density gradients by the diffusionFootnote ⁶ . As a summary of these reduced INMDF fluid numerical simulations, we prove here the origin of the diffusion via kinetic effects (i.e. via the existence of the INMDF due to spatial gradients) on shorter time scales than those of the thermalization. We obtained similar numerical results for the heat conduction with $\unicode[STIX]{x2202}_{t}(nT)=-3\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6E4}W)$ and $\unicode[STIX]{x2202}_{t}(\unicode[STIX]{x1D6E4}W)=-\unicode[STIX]{x1D735}(nT^{2})$ where we can choose constant and homogeneous $n$ and $W$ . As a first improvement of these models limited by the artificial thermalization by forcing $\unicode[STIX]{x1D6E4}=0$ , the implementation of the Bhatnagar–Gross–Krook (also called BGK or Krook) collision operator (Gross & Krook Reference Gross and Krook1956) in these minimal reduced INMDF fluid models is under investigation and will be reported in future publications, as well as the numerical and analytic evaluation of the Fokker–Planck collision operator (see § 3.2). These results are unique thanks to the use of the INMDF since by using the momentum $nv$ of a MDF instead of $\unicode[STIX]{x1D6E4}$ there is no physical motivation to reset the momentum to $0$ on the time scales of the thermalization. These observations prove that (i) the kinetic heat flux $\unicode[STIX]{x1D6E4}$ exists at least in part due to the existence of spatial density gradients, (ii) the usual static particle diffusion and heat conduction are not consistent with dynamical density and temperature gradients (and with our dynamical growth rate of the kinetic heat flux) and (iii) the INMDF has a physical existence via our understanding of the origin of particle and heat diffusionFootnote ⁷ . These significant observations based on very simple INMDF fluid models are the keys that allow us to be optimistic for a potential understanding of the turbulence transition. These results are further discussed in the perspectives in § 4 and will be the purpose of future publications.

Until now we have explained only the collisionless fluid closure. The next step is to introduce the collisional fluid closure that is responsible of the relaxation toward the MDF equilibrium.

3.2 Local collisional fluid closures

There are two types of collisional fluid closures: the self-collisions of the particle species with themselves and the collisions with energetic particles coming from a source of energy (e.g. NBI, radio frequency waves, alpha particles, other species). The former tends to reduce the deviation from the MDF which is the thermodynamic equilibrium without sources, but the presence of the latter constantly generates deviations from the MDF and a steady-state NMDF is possible at finite collisionality. Here we argue that this steady state can be described by the kinetic theory as well as by the generalized fluid theory introduced by our work.

For the first type of collisional fluid closure, the nonlinear Fokker–Planck collisional operator is required for the description of low angle scattering by two-body collisions. By assuming that the steady-state distribution function is an INMDF $f_{I}$ instead of a MDF, we compute the moments of the self-collision Fokker–Planck operator

and the fluid equations read

where $C_{k}$ are functions of the hidden variables of the INMDF. We found (Izacard Reference Izacard2016a ) that the terms $C_{k}$ introduce a new function $I_{k}(a,b)$ , called the plasma collision function, that is defined by

In fact, the nonlinear Fokker–Planck collision operator can be written in a form involving only one Rosenbluth potential (Rosenbluth et al. Reference Rosenbluth, Macdonald and Judd1957) as shown by Gaffey (Reference Gaffey1976). The correction of the Rosenbluth potential from $F(x)=1/(nv)\int _{-\infty }^{\infty }$ $\times |\boldsymbol{v}_{i}-\boldsymbol{v}|f(\boldsymbol{v})\,\text{d}^{3}\boldsymbol{v}$ where $x=\text{v}_{i}/v$ for a MDF $f=f_{0}$ (i.e. the correction from $F_{0}(x)=(x+1/(2x))\text{Erf}(x)+\exp (-x)/\sqrt{\unicode[STIX]{x03C0}}$ ) introduces the plasma collision function $I_{k}(a,b)$

where $A=\text{v}_{i\bot }^{2}/(2W)$ and $B=\text{v}_{i\bot }(c-\text{v}_{iz})/(2W)$ . Moreover, even an anisotropic MDF (e.g. with different perpendicular $T_{\bot }$ and parallel $T_{\Vert }$ temperatures) introduces this plasma collision function. In order to obtain the correction $F_{I}(x)$ , an anisotropic three-dimensional version (Izacard Reference Izacard2016a ) of the INMDF is required. We remark that the tensor version of the fluid equations are given in appendix A. More details of the development of the nonlinear Fokker–Planck collision operator from a three-dimensional version (Izacard Reference Izacard2016a ) of the INMDF will appear in a future publication since the focus of this work is to understand the limitations of the current fluid theory and to explain the key points allowing the self-consistent generalization of the fluid theory including kinetic effects.

For the second type of collisional fluid closure we need to add as a source or sink of energy the collision with other species. The distribution function of these other species can be assumed or evolved by the Fokker–Planck equation. In the latter case, if there are no external sources in both Fokker–Planck equations, both steady states can only be Maxwellians and are not relevant to experiments. However, for example, the addition of a NBI in one or the other Fokker–Planck equations is the key to obtaining NMDFs steady states for both species. The source term, by assuming an INMDF $f_{I}$ as a first species and a source $f_{s}$ as a second species, is defined by

and the fluid equations become

Following the techniques published by Gaffey (Reference Gaffey1976), it could be possible to derive the steady-state distribution function $f_{s}$ of a NBI as a response to the INMDF background plasma, and then to compute the Fokker–Planck collision operator and the source term given by (3.13). From this computation, we should be able to check if the INMDF is the first analytic self-consistent solution of the Boltzmann equation in the presence of sources such as NBI.

In summary, we have demonstrated the method to construct generalized fluid models including kinetic effects. We have made the distinction between the collisionless and collisional fluid closures that are both local and defined by the velocity phase-space representation of the distribution function. For the next generation of fluid codes, we have shown an algorithm that evolves the dynamics of the fluid hidden variables, defining both the distribution function and the fluid moments. The unification of kinetic and fluid theories is obtained. The highlights of some results and a list of perspectives are detailed in § 4.