2 Solar wind heating: a huge problem

As introduced above, the solar wind is a weakly collisional, strongly turbulent medium (Bruno & Carbone Reference Bruno and Carbone2013). Several observations indicate that the solar wind is incessantly heated during its travel through the heliosphere: the temperature decay along the radial distance is indeed much slower than the decay expected with adiabatic models of the wind expansion (Marsch et al. Reference Marsch, Muhlhauser, Schwenn, Rosenbauer, Pilipp and Neubauer1982; Goldstein et al. Reference Goldstein, Neugebauer, Phillips, Bame, Gosling, McComas, Wang, Sheeley and Suess1996; Marino et al. Reference Marino, Sorriso-Valvo, Carbone, Noullez, Bruno and Bavassano2008; Cranmer et al. Reference Cranmer, Matthaeus, Breech and Kasper2009). Therefore, some local heating mechanisms must play a significant role to supply the energy needed to heat the plasma. Numerous scenarios have been proposed to understand the plasma heating and a long-standing debate about which processes are preferred is still waiting for a clear and definitive answer (see Bruno & Carbone (Reference Bruno and Carbone2013) and references therein). Among these processes, it is widely known that the turbulence efficiently contributes to the local heating of the solar wind (Sorriso-Valvo et al. Reference Sorriso-Valvo, Marino, Carbone, Noullez, Lepreti, Veltri, Bruno, Bavassano and Pietropaolo2007; Marino et al. Reference Marino, Sorriso-Valvo, Carbone, Noullez, Bruno and Bavassano2008; Sahraoui et al. Reference Sahraoui, Goldstein, Robert and Khotyaintsev2009), since it can efficiently transfer a significant amount of energy towards smaller scales, where dissipative mechanisms are at work. In fact, in a turbulent flow, much more energy is transferred towards smaller scales with respect to a laminar flow: the ratio between the energy transfer flux due to turbulence at a certain scale with respect to the heating production due to dissipation at the same scale is proportional to the Reynolds number $Re$ , thus indicating that the energy transfer towards smaller scales gets more efficient as the flow becomes more turbulent.

In the simple neutral fluid scenario, the cascade is arrested once that the dissipative scale is reached (Frisch Reference Frisch1995). On the other hand, the cascade evolves in a more complex way in a plasma: the presence of other processes (for example dispersion and kinetic effects) strongly modify the cascade before reaching the dissipative scale. A relatively wide agreement has been achieved concerning the importance of turbulence for transferring energy towards smaller scales and many scenarios have been proposed to explain the transition from the inertial range towards the kinetic scales and the nature of dissipative processes. These scenarios are often based on the ‘collisionless’ assumption, that is justified by the fact that the Spitzer–Harm collisional time (Spitzer Reference Spitzer1956) is much larger than other dynamical times. We would remark that two important caveats should be considered.

First, any mechanism which does not consider collisions is not able to describe the last part of the heating process, namely the heat production due to the irreversible dissipation of phase space structures and the approach towards the thermal equilibrium. For example, several mechanisms (e.g. nonlinear waves) can indeed increase the particle temperature, evaluated as the second-order moment of the particle distribution function, by producing non-Maxwellian features as beams of trapped particles. However, this temperature growth due to the beam production does not represent a temperature growth in the thermodynamic sense, because the beam presence makes the system out of equilibrium. The particle beam can be instead interpreted as a form of free energy stored into the VDF. This energy is not in general converted into heat by means of irreversible processes but it can be also transformed in other forms of ordered energy (e.g. through micro-instabilities) (Lesur, Diamond & Kosuga Reference Lesur, Diamond and Kosuga2014). Collisions are the unique mechanism able to degrade this information into heat by approaching the thermal equilibrium, thus producing heating in the general thermodynamic and irreversible sense.

Second, the evaluation of the Spitzer–Harm collisional time strictly assumes that the VDF shape is close to the equilibrium Maxwellian. This assumption may not be true in the solar wind (Marsch Reference Marsch2006; Servidio et al. Reference Servidio, Valentini, Perrone, Greco, Califano, Matthaeus and Veltri2015), where the VDF’s shape is strongly perturbed by kinetic turbulence. In this direction, by focusing on the collisional relaxation in a homogeneous force-free plasma where collisions are modelled with the fully nonlinear Landau operator (Landau Reference Landau1936), we recently showed that fine velocity structures are dissipated much faster than global non-thermal features such as temperature anisotropy (Paper I). The entropy production due to the relaxation of the VDF towards the equilibrium occurs on several characteristic times. These characteristic times are associated with the dissipation of particular velocity space structures and can be much smaller than the Spitzer–Harm time (Spitzer Reference Spitzer1956), this indicating that collisions could effectively compete with other processes (e.g. micro-instabilities). In this perspective, high-resolution measurements of the particle VDF in the solar wind are crucial for a proper description of the heating problem (Vaivads et al. Reference Vaivads, Retinó, Soucek, Khotyaintsev, Valentini, Escoubet, Alexandrova, André, Bale and Balikhin2016).

In principle the combination of the turbulent nature of the solar wind with its weakly collisionality may constitute a new scenario to describe solar wind heating. In fact, turbulence is able to transfer energy towards smaller scales. Then, when kinetic scales are reached, since the plasma is weakly collisional, the VDF becomes strongly distorted and exhibits non-Maxwellian features, such as beams, anisotropies and ring-like structures (Belmont et al. Reference Belmont, Mottez, Chust and Hess2008; Chust et al. Reference Chust, Belmont, Mottez and Hess2009; Servidio et al. Reference Servidio, Valentini, Califano and Veltri2012, Reference Servidio, Valentini, Perrone, Greco, Califano, Matthaeus and Veltri2015). The presence of strong gradients in velocity space tends to naturally enhance the effect of collisions, which – ultimately – may become efficient for dissipating these structures and for producing heat.

Based on these last considerations, numerous studies have been recently conducted in order to take into account collisional effects in a weakly collisional plasma such as the solar wind (Filbet & Pareschi Reference Filbet and Pareschi2002; Bobylev & Potapenko Reference Bobylev and Potapenko2013; Pezzi et al. Reference Pezzi, Valentini, Perrone and Veltri2013, Reference Pezzi, Valentini, Perrone and Veltri2014; Escande, Elskens & Doveil Reference Escande, Elskens and Doveil2015; Pezzi, Valentini & Veltri Reference Pezzi, Valentini and Veltri2015b ; Pezzi et al. Reference Pezzi, Valentini and Veltri2016a ; Banón Navarro et al. Reference Banón Navarro, Teaca, Told, Groselj, Crandall and Jenko2016; Hirvijoki et al. Reference Hirvijoki, Lingam, Pfefferlé, Comisso, Candy and Bhattacharjee2016; Tigik et al. Reference Tigik, Ziebell, Yoon and Kontar2016), where collisions are usually introduced by means of a collisional operator at the right-hand side of the Vlasov equation. The choice of the proper collisional operator remains an open problem. Several derivations from first principles (e.g. the Liouville equation) indicate that the most general collisional operators for plasmas are the Lenard–Balescu operator (Balescu Reference Balescu1960; Lenard Reference Lenard1960) or the Landau operator (Landau Reference Landau1936; Akhiezer et al. Reference Akhiezer, Akhiezer, Polovin, Sitenko and Stepanov1986). Both operators are nonlinear ‘Fokker–Planck’-like operators which involve velocity space derivatives and three-dimensional integrals. The Landau operator introduces an upper cutoff of the integrals at the Debye length to avoid the divergence for large impact parameters, while the Balescu–Lenard operator solves this divergence in a more consistent way through the dispersion equation. Therefore, the Balescu–Lenard operator is more general than the Landau operator from this point of view. However, both operators are derived by assuming that the plasma is not extremely far from the thermal equilibrium. Hence, both operators could lack the description of interparticle collisions in a strongly turbulent system. The numerical approach of operators is also much more difficult for the Balescu–Lenard operator with respect to the Landau operator, because it involves the evaluation of the dispersion function. Finally, we would also point out that, as far as we know, an explicit derivation of the Boltzmann operator for plasmas starting from the Liouville equation does not exist (Villani Reference Villani2002). Despite the fact that the adoption of the Boltzmann operator for describing collisional effects in plasmas is questionable from a theoretical perspective, it still represents a valid option since the Boltzmann and Fokker–Planck-like operators such as the Landau operator are intrinsically similar (Landau Reference Landau1936; Bobylev & Potapenko Reference Bobylev and Potapenko2013).

The computational cost to evaluate both the Landau and Balescu–Lenard operators numerically is huge: for $N$ grid points along each direction of the $3D$ – $3V$ numerical phase space (three dimensions in physical space and three dimensions in velocity space), the computation of the Landau operator would require approximately $N^{9}$ operations at each time step. In fact, for each point of the six-dimensional grid, a three-dimensional integral must be computed. To avoid this numerical complexity, several simplified operators have been proposed. We may distinguish these simpler operators in two classes. The first type of operators, such as the Bathanar–Gross–Krook (Bhatnagar, Gross & Krook Reference Bhatnagar, Gross and Krook1954; Livi & Marsch Reference Livi and Marsch1986) and the Dougherty operators (Dougherty Reference Dougherty1964; Dougherty & Watson Reference Dougherty and Watson1967; Pezzi, Valentini & Veltri Reference Pezzi, Valentini and Veltri2015a ), model collisions in the realistic three-dimensional velocity space by adopting a simpler structure of the operator. The second class of collisional operators works instead in a reduced, one-dimensional velocity space assuming that the dynamics mainly occurs in one direction. Although this approach is quite ‘unphysical’ (collisions naturally act in three dimensions), these operators can satisfactorily model collisions in laboratory plasmas devices, such as the Penning–Malmberg traps, where the plasma is confined into a long and thin column and the dynamics occurs mainly along a single direction (Anderson & O’Neil Reference Anderson and O’Neil2007a ,Reference Anderson and O’Neil b ; Pezzi et al. Reference Pezzi, Valentini, Perrone and Veltri2013; Pezzi, Camporeale & Valentini Reference Pezzi, Camporeale and Valentini2016b ).

3 Numerical approach and simulation results

As described above, to highlight the importance of nonlinearities present in the collisional operator, here we compare the effects of using the full Landau operator with a model of the linearized Landau operator, obtained by simplifying the structures of the Landau operator coefficients. We restrict ourselves to the case of a force-free homogeneous plasma and we just model collisions between particles of the same species. Our interest is in fact to understand how collisional effects change when the mathematical kernel of the collisional operator is modified. Based on these assumptions, we numerically integrate the following dimensionless collisional evolution equations for the particle distribution function $f(\boldsymbol{v},t)$ :

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}f(\boldsymbol{v},t)}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x03C0}\left(\frac{3}{2}\right)^{3/2}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{i}}\int \text{d}^{3}v^{\prime }\,U_{ij}(\boldsymbol{u})\left[f(\boldsymbol{v}^{\prime },t)\frac{\unicode[STIX]{x2202}f(\boldsymbol{v},t)}{\unicode[STIX]{x2202}v_{j}}-f(\boldsymbol{v},t)\frac{\unicode[STIX]{x2202}f(\boldsymbol{v}^{\prime },t)}{\unicode[STIX]{x2202}v_{j}^{\prime }}\right], & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}f(\boldsymbol{v},t)}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x03C0}\left(\frac{3}{2}\right)^{3/2}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{i}}\int \text{d}^{3}v^{\prime }\,U_{ij}(\boldsymbol{u})\left[f_{0}(\boldsymbol{v}^{\prime })\frac{\unicode[STIX]{x2202}f(\boldsymbol{v},t)}{\unicode[STIX]{x2202}v_{j}}-f(\boldsymbol{v},t)\frac{\unicode[STIX]{x2202}f_{0}(\boldsymbol{v}^{\prime })}{\unicode[STIX]{x2202}v_{j}^{\prime }}\right]. & \displaystyle\end{eqnarray}$$

where $f$ is normalized such that $\int \text{d}^{3}v\,f(\boldsymbol{v})=n=1$ and $U_{ij}(\boldsymbol{u})$ is

(3.3) $$\begin{eqnarray}U_{ij}(\boldsymbol{u})=\frac{\unicode[STIX]{x1D6FF}_{ij}u^{2}-u_{i}u_{j}}{u^{3}},\end{eqnarray}$$

where $\boldsymbol{u}=\boldsymbol{v}-\boldsymbol{v}^{\prime }$ , $u=|\boldsymbol{u}|$ and the Einstein notation is introduced. In (3.1)–(3.2), and from now on, time is scaled to the inverse Spitzer–Harm frequency $\unicode[STIX]{x1D708}_{SH}^{-1}$ (Spitzer Reference Spitzer1956) and velocity to the particle thermal speed $v_{th}$ . Details about the numerical solution of (3.1)–(3.2) can be found in (Pezzi et al. Reference Pezzi, Valentini and Veltri2015a , Reference Pezzi, Valentini and Veltri2016a ). In (3.2), $f_{0}(\boldsymbol{v})$ is the three-dimensional Maxwellian distribution function associated with the initial condition of our simulations $f(\boldsymbol{v},t=0)$ and is built in such a way that density, bulk velocity and temperature of the two distributions $f(\boldsymbol{v},t=0)$ and $f_{0}(\boldsymbol{v})$ are equal. The two equations clearly differ because (3.2) is a linearized model of (3.1). The operator described in (3.2) has been in fact obtained by linearizing the coefficients of the Landau operator. Although this linear operator does not represent the exact linearization of the Landau operator, the procedure here utilized for linearizing the operator (i.e. simplifying only the Fokker–Planck coefficients) is commonly adopted. In the following, we will note that the simulations performed with the linearized operator thermalize to the same final VDF and produce also the same total entropy growth as the nonlinear simulations. This suggests that the term which is not included in the form collisional operator (whose Fokker–Planck coefficients depend on $(f-f_{0})(\boldsymbol{v}^{\prime })$ ) is not extremely relevant in the global thermalization of the system. This approximation corresponds to retain the gradients related to the out-of-equilibrium structures but to neglect their contribution to the integral in the $\boldsymbol{v}^{\prime }$ space. In other words, here we locally consider gradients but we neglect their contribution to the global Fokker–Planck coefficients.

When simulations are completed, we perform the following multi-exponential fit (Curtis, Berry & Bromander Reference Curtis, Berry and Bromander1970; Pezzi et al. Reference Pezzi, Valentini and Veltri2016a ) of the entropy growth $\unicode[STIX]{x0394}S$ to point out the presence of several characteristic times:

(3.4) $$\begin{eqnarray}\unicode[STIX]{x0394}S(t)=\mathop{\sum }_{i=1}^{K}\unicode[STIX]{x0394}S_{i}(1-\text{e}^{-t/\unicode[STIX]{x1D70F}_{i}}),\end{eqnarray}$$

$\unicode[STIX]{x1D70F}_{i}$ being the $i$ th characteristic time, $\unicode[STIX]{x0394}S_{i}$ the growth of entropy related to the characteristic time $\unicode[STIX]{x1D70F}_{i}$ and $K$ is evaluated through a recursive procedure. This procedure has been already adopted in Paper I to highlight the importance of fine velocity structures in the entropy growth. In the following subsections, we report and describe the results of the simulations performed with two different initial distribution functions, already adopted in Paper I. The first initial condition concerns the presence of non-Maxwellian signatures due to a strongly nonlinear wave – an electron acoustic wave (EAW) (Holloway & Dorning Reference Holloway and Dorning1991; Kabantsev, Valentini & Driscoll Reference Kabantsev, Valentini and Driscoll2006; Valentini, ONeil & Dubin Reference Valentini, ONeil and Dubin2006; Anderegg et al. Reference Anderegg, Driscoll, Dubin and ONeil2009a ,Reference Anderegg, Driscoll, Dubin, ONeil and Valentini b ; Johnston et al. Reference Johnston, Tyshetskiy, Ghizzo and Bertrand2009; Valentini et al. Reference Valentini, Perrone, Califano, Pegoraro, Veltri, Morrison and O’Neil2012) – in the core of the distribution function. The EAWs here excited are quite different from another type of electron acoustic fluctuations that occur in a plasma composed of two components at different temperatures (Watanabe & Taniuti Reference Watanabe and Taniuti1977) and can be also observed in the Earth’s magnetosphere (Tokar & Gary Reference Tokar and Gary1984; Lu, Wang & Wang Reference Lu, Wang and Wang2005). The EAWs here excited are undamped waves whose phase speed is close to the thermal speed. It has been shown that, in the usual theory of an equilibrium Maxwellian plasma, these waves are strongly damped; while they can survive if the distribution function is locally modified (and exhibits a flat region) around the wave phase velocity. To generate the nonlinearity in the distribution function and allow these waves to survive, external drivers are usually adopted to force the plasma. EAWs are also characterized by the presence of phase space Bernstein–Green–Kruskal (BGK) structures (Bernstein, Greene & Kruskal Reference Bernstein, Greene and Kruskal1957) in the core of the electron distribution function, associated with trapped particle populations. The second initial distribution is instead a typical VDF recovered in hybrid Vlasov–Maxwell simulations of solar wind decaying turbulence (Servidio et al. Reference Servidio, Valentini, Califano and Veltri2012; Valentini et al. Reference Valentini, Servidio, Perrone, Califano, Matthaeus and Veltri2014; Servidio et al. Reference Servidio, Valentini, Perrone, Greco, Califano, Matthaeus and Veltri2015). The two simulations from which we selected our initial VDF are quite different. Indeed, in the first case, the out-of-equilibrium structures present in the initial VDF are due to the wave–particle interaction with the EAW, which is an almost monochromatic (few excited wavenumbers) electrostatic wave. On the other hand, in the second case, the initial distribution function has been strongly distorted due to the presence of an electromagnetic, turbulent cascade.

3.1 First case study: wave–particle interactions and collisions

Figure 1. (a) Power spectral density of the electric energy $E_{E}(k_{z})$ as a function of the wavenumber $k_{z}$ . (b) Contour plot of $f_{e}(z,v_{z})$ at the time instant when the EAWs is fully developed. The red line represents the coordinate $z=z_{0}$ where the cut is performed. (c) Profile of $f_{e}(z_{0},v_{z})$ as a function of $v_{z}$ .

The first initial condition here adopted, which is a three-dimensional VDF that evolves according to (3.1)–(3.2) in the three-dimensional velocity space, has been designed as follows. We separately performed a $1D$ – $1V$ Vlasov-Poisson simulation of a electrostatic plasma composed of kinetic electrons and motionless protons whose resolution, in the $z-v_{z}$ phase space domain, is $N_{z}=256$ , $N_{v_{z}}=1601$ . In order to excite a large amplitude EAW, we forced the system with an external sinusoidal electric field, which has been adiabatically turned on and off to properly trigger the wave. Figure 1(a) reports the power spectral density of the electric energy $E_{E}(k_{z})$ as a function of the wavenumber $k_{z}$ , evaluated at the final time instant of the Vlasov–Poisson simulation (where the EAW is fully developed). Few wavenumbers are significantly excited and the EAW is almost monochromatic. The features of the electric fluctuations spectrum are reflected into the shape of the distribution function, which is locally distorted around the phase speed and presents a clear BGK hole, as reported in figure 1(b). Since the grid size in velocity space is quite small in the current simulation, relatively small velocity scales are dynamically generated during the simulation by wave–particle interaction.

Figure 2. Time history of $\unicode[STIX]{x0394}S$ in the case of the fully nonlinear Landau operator (black) and the linearized Landau operator (red). Blue diamonds indicate the time instants $t=T_{nl,1}=\unicode[STIX]{x1D70F}_{1}^{nl}$ , $t=T_{nl,2}=\unicode[STIX]{x1D70F}_{1}^{nl}+\unicode[STIX]{x1D70F}_{2}^{nl}$ and $t=T_{nl,3}=\unicode[STIX]{x1D70F}_{1}^{nl}+\unicode[STIX]{x1D70F}_{2}^{nl}+\unicode[STIX]{x1D70F}_{3}^{nl}$ ; the green triangles refer to $t=T_{lin,1}=\unicode[STIX]{x1D70F}_{1}^{lin}$ , $t=T_{lin,2}=\unicode[STIX]{x1D70F}_{1}^{lin}+\unicode[STIX]{x1D70F}_{2}^{lin}$ and $t=T_{lin,3}=\unicode[STIX]{x1D70F}_{1}^{lin}+\unicode[STIX]{x1D70F}_{2}^{lin}+\unicode[STIX]{x1D70F}_{3}^{lin}$ .

Then, we selected the spatial point $z_{0}$ in the numerical domain (red vertical line in figure 1 b), where this BGK-like phase space structure displays its maximum velocity width, and we considered the velocity profile $\hat{f}_{e}(v_{z})=f_{e}(z_{0},v_{z})$ , whose shape as a function of $v_{z}$ is reported in figure 1(c). $\hat{f}_{e}$ is highly distorted due to nonlinear wave–particle interactions and exhibits sharp velocity gradients (bumps, holes, spikes around the resonant speed). Finally, by evaluating the density $n_{e}$ , the bulk speed $V_{e}$ and the temperature $T_{e}$ of $\hat{f}_{e}$ , we built up the three-dimensional VDF $f(v_{x},v_{y},v_{z})=f_{M}(v_{x})f_{M}(v_{y})\hat{f}_{e}(v_{z})$ , which represents our initial condition, with $f_{M}$ the one-dimensional Maxwellian associated with $\hat{f}_{e}$ . We remark that this VDF does not exhibit any temperature anisotropy but it still exhibits strong non-Maxwellian deformations along $v_{z}$ , due to the presence of trapped particles, which make the system far from thermal equilibrium. The three-dimensional velocity domain is here discretized by $N_{v_{x}}=N_{v_{y}}=51$ and $N_{v_{z}}=1601$ grid points in the region $v_{i}=[-v_{max},v_{max}]$ , with $v_{max}=6v_{th}$ and $i=x,y,z$ , while boundary conditions assume that the distribution function is set equal to zero for $|v_{j}|>v_{max}$ .

Since no temperature anisotropies are present, the evolution of the total temperature and of the temperatures along each direction are trivial, since the total temperature is preserved. On the other hand, the evolution of the entropy variation $\unicode[STIX]{x0394}S=S(t)-S(0)$ ( $S=-\int f\ln f\,\text{d}^{3}v$ ) gives information about the approach towards equilibrium. The time history of $\unicode[STIX]{x0394}S$ obtained with the nonlinear Landau operator (black) and with the linearized Landau operator (red) are shown in figure 2. Since the initial condition and the equilibrium Maxwellian reached under the effect of collisions is the same for both operators, the total entropy growth $\unicode[STIX]{x0394}S$ is the same in the two cases. In other words, the free energy contained in the out-of-equilibrium structures in the initial VDF produces the same entropy growth in absolute terms but the growth occurs on different time scales in the two cases. Indeed, in the nonlinear operator case the entropy grows much faster ( $\simeq 1\unicode[STIX]{x1D708}_{SH}^{-1}$ ) compared to the linearized operator case ( $\simeq 4\unicode[STIX]{x1D708}_{SH}^{-1}$ ).

Figure 3. Distribution function $f(v_{x}=0,v_{y}=0,v_{z})$ as a function of $v_{z}$ , obtained in the case of the fully nonlinear Landau operator. Panels from (a) to (d) respectively display the time instants $t=T_{nl,1}=\unicode[STIX]{x1D70F}_{1}^{nl}$ (a), $t=T_{nl,2}=\unicode[STIX]{x1D70F}_{1}^{nl}+\unicode[STIX]{x1D70F}_{2}^{nl}$ (b), $t=T_{nl,3}=\unicode[STIX]{x1D70F}_{1}^{nl}+\unicode[STIX]{x1D70F}_{2}^{nl}+\unicode[STIX]{x1D70F}_{3}^{nl}$ (c) and $t=t_{fin}$ (d). Results for the linearized Landau operator are qualitatively the same with the primary difference being that $\unicode[STIX]{x1D70F}_{j}^{nl}$ is substantially smaller than $\unicode[STIX]{x1D70F}_{j}^{lin}$ for each value of $j$ .

To quantify the presence of several characteristic times, we perform the multi-exponential fit of (3.4) on the entropy growth curves reported in figure 2. The analysis of the growth recovered in the fully nonlinear Landau operator indicates that three different characteristic times are recovered in the entropy growth:

1. (i) $\unicode[STIX]{x1D70F}_{1}^{nl}=3.5\times 10^{-3}\,\unicode[STIX]{x1D708}_{SH}^{-1}\rightarrow \unicode[STIX]{x0394}S_{1}^{nl}/\unicode[STIX]{x0394}S_{tot}=13\,\%$ ,

2. (ii) $\unicode[STIX]{x1D70F}_{2}^{nl}=1.3\times 10^{-1}\,\unicode[STIX]{x1D708}_{SH}^{-1}\rightarrow \unicode[STIX]{x0394}S_{2}^{nl}/\unicode[STIX]{x0394}S_{tot}=42\,\%$ ,

3. (iii) $\unicode[STIX]{x1D70F}_{3}^{nl}=4.9\times 10^{-1}\,\unicode[STIX]{x1D708}_{SH}^{-1}\rightarrow \unicode[STIX]{x0394}S_{3}^{nl}/\unicode[STIX]{x0394}S_{tot}=40\,\%$ .

As discussed in Paper I, the presence of several characteristic times is associated with the dissipation of different velocity space structures. Figure 3 reports $f(v_{x}=v_{y}=0,v_{z})$ as a function of $v_{z}$ at the time instants $t=T_{nl,1}=\unicode[STIX]{x1D70F}_{1}^{nl}$ (a), $t=T_{nl,2}=\unicode[STIX]{x1D70F}_{1}^{nl}+\unicode[STIX]{x1D70F}_{2}^{nl}$ (b), $t=T_{nl,3}=\unicode[STIX]{x1D70F}_{1}^{nl}+\unicode[STIX]{x1D70F}_{2}^{nl}+\unicode[STIX]{x1D70F}_{3}^{nl}$ (c) and $t=t_{fin}$ (d), These time instants are displayed in figure 2 with blue diamonds. After the time $t=T_{nl,1}=\unicode[STIX]{x1D70F}_{1}^{nl}$ (a), steep spikes visible in figure 1(b) are almost completely smoothed out; then, at time $t=T_{nl,2}=\unicode[STIX]{x1D70F}_{1}^{nl}+\unicode[STIX]{x1D70F}_{2}^{nl}$ (b), the remaining plateau region is significantly rounded off, only a gentle shoulder being left; finally, after a time $t=T_{nl,3}=\unicode[STIX]{x1D70F}_{1}^{nl}+\unicode[STIX]{x1D70F}_{2}^{nl}+\unicode[STIX]{x1D70F}_{3}^{nl}$ (c), the collisional relaxation to equilibrium is completed for the most part. A small percentage $\simeq 5\,\%$ of the total entropy growth is finally recovered for larger times and corresponds to the final approach to the equilibrium Maxwellian (d).

By performing the same analysis for the linearized Landau operator case, three characteristic times are also recovered:

1. (i) $\unicode[STIX]{x1D70F}_{1}^{lin}=1.1\times 10^{-2}\,\unicode[STIX]{x1D708}_{SH}^{-1}\rightarrow \unicode[STIX]{x0394}S_{1}^{lin}/\unicode[STIX]{x0394}S_{tot}=11\,\%$ ,

2. (ii) $\unicode[STIX]{x1D70F}_{2}^{lin}=2.7\times 10^{-1}\,\unicode[STIX]{x1D708}_{SH}^{-1}\rightarrow \unicode[STIX]{x0394}S_{2}^{lin}/\unicode[STIX]{x0394}S_{tot}=23\,\%$ ,

3. (iii) $\unicode[STIX]{x1D70F}_{3}^{lin}=1.5\,\unicode[STIX]{x1D708}_{SH}^{-1}\rightarrow \unicode[STIX]{x0394}S_{3}^{lin}/\unicode[STIX]{x0394}S_{tot}=63\,\%$ .

These characteristic times are systematically larger than the times recovered in the nonlinear operator case. The shape of the distribution function after each characteristic time (not shown here) is quite similar to the shape recovered in the case of the fully nonlinear operator evolution. The process of dissipation of fine velocity structure is, hence, qualitatively similar if one adopts nonlinear or linearized operators. However, significant quantitative differences occur: similar profiles in velocity space are indeed reached at very different times, with the characteristic times recovered in the linearized case significantly larger (approximately 4 times) than the times recovered in the nonlinear operator case.

Therefore, from a qualitative point of view, both operators are able to recover the fact that fine velocity space structures are dissipated faster as their scale become finer (i.e. as the velocity space gradients become stronger). However, fine velocity structures are dissipated slower by linearizing the collisional operator. Moreover, it is also worth mentioning that the amount of entropy growth associated with each characteristic time slightly changes by ignoring nonlinearities. For example, in the case of the fully nonlinear Landau operator, approximately 55 % of the total entropy growth is produced when the initial spikes and the successive flat plateau are dissipated. On the other hand, in the linearized operator case, only approximately the 30 % of the total entropy growth is associated with these processes.

Figure 4. (a) Omni-directional power spectral densities (PSDs) of the magnetic energy $E_{B}(k_{\bot })$ (black line) and of the electric energy $E_{E}(k_{\bot })$ (red line) as a function of the perpendicular wavenumber $k_{\bot }$ . PSDs have been evaluated at the time instant where the turbulent activity is maximum. (b) Iso-surface of the initial distribution function. Red, green and blue axes refer to $v_{x}$ , $v_{y}$ and $v_{z}$ , respectively.

3.2 Second case study: kinetic turbulence and collisions

To support the scenario described in the previous section, here we focus on a second initial condition. This initial VDF has been selected from a $2D$ – $3V$ hybrid Vlasov-Maxwell numerical simulation of decaying turbulence in solar wind-like conditions (Valentini et al. Reference Valentini, Travnicek, Califano, Hellinger and Mangeney2007, Reference Valentini, Servidio, Perrone, Califano, Matthaeus and Veltri2014). The hybrid Vlasov–Maxwell simulation, whose resolution is $N_{x}=N_{y}=512$ and $N_{v_{x}}=N_{v_{y}}=N_{v_{z}}=51$ , is initialized with an out-of-the-plane background magnetic field. Then, magnetic and bulk speed perturbations at large, MHD scales are introduced. As a result of nonlinear couplings among the fluctuations, the energy cascades towards smaller kinetic scales. Hence, the particle VDF strongly departs from the thermal equilibrium due to the presence of kinetic turbulence and exhibits a potato-like shape similar to solar wind in situ observations (Marsch Reference Marsch2006). The omni-directional power spectral densities of the magnetic (black line) and electric energy (red line), evaluated at the time instant where the turbulent activity is maximum, are reported in figure 4(a): clearly a broadband spectrum is recovered. The iso-contour of the initial VDF, selected where non-Maxwellian effects are strongest (Servidio et al. Reference Servidio, Valentini, Perrone, Greco, Califano, Matthaeus and Veltri2015), is shown in figure 4(b). The VDF exhibits a hole-like structure in the upper part of the box and a thin ring-like structure in the bottom part of the box and is clearly elongated on the $v_{z}$ direction. Compared to the first case study, the current distribution function reflects the presence of a spectrum of excited wavenumbers and it contains several kinds of distortions, not only concentrated around the resonant speed as in the previous case.

In Paper I we showed that, once velocity space gradients are artificially smoothed out through a fitting procedure, the presence of several characteristic times associated with the dissipation of fine velocity structures is definitively lost. Here, we instead compare the evolution towards the equilibrium of this initial condition under the effect of the fully nonlinear Landau operator (3.1) and the linearized Landau operator (3.2). The velocity domain is here discretized with $N_{v_{x}}=N_{v_{y}}=N_{v_{z}}=51$ points. Note that, compared to the first case study, the resolution is here lower and it cannot be increased due to the computational cost of the hybrid Vlasov–Maxwell code. Therefore, the quite small velocity scales recovered in the first case study (spikes around the resonant speed etc.) are not present in this case.

Figure 5 reports the entropy growth obtained with the fully nonlinear Landau operator (black) and with its linearized version (red). The entropy growth is, also here, slower in the linearized operator case compared to the fully nonlinear operator case. To quantify the different evolution observed in figure 5, we perform the multi-exponential fit (Curtis et al. Reference Curtis, Berry and Bromander1970; Pezzi et al. Reference Pezzi, Valentini and Veltri2016a ) described in (3.4).

Figure 5. Time history of $\unicode[STIX]{x0394}S$ in the case of the fully nonlinear Landau operator (black) and the linearized Landau operator (red). Blue diamonds indicate the time instants $t=T_{nl,1}=\unicode[STIX]{x1D70F}_{1}^{nl}$ and $t=T_{nl,2}=\unicode[STIX]{x1D70F}_{1}^{nl}+\unicode[STIX]{x1D70F}_{2}^{nl}$ ; the green triangles refer to $t=T_{lin,1}=\unicode[STIX]{x1D70F}_{1}^{lin}$ and $t=T_{lin,2}=\unicode[STIX]{x1D70F}_{1}^{lin}+\unicode[STIX]{x1D70F}_{2}^{lin}$ .

The analysis performed in the case of the fully nonlinear operator indicates that the entropy grows with two characteristic times:

1. (i) $\unicode[STIX]{x1D70F}_{1}^{nl}=0.20\,\unicode[STIX]{x1D708}_{SH}^{-1}\rightarrow \unicode[STIX]{x0394}S_{1}^{nl}/\unicode[STIX]{x0394}S_{tot}=26\,\%$

2. (ii) $\unicode[STIX]{x1D70F}_{2}^{nl}=0.82\,\unicode[STIX]{x1D708}_{SH}^{-1}\rightarrow \unicode[STIX]{x0394}S_{2}^{nl}/\unicode[STIX]{x0394}S_{tot}=74\,\%$

As in the case described in the previous section, each characteristic time is associated with the dissipation of a different out-of-equilibrium features. Figure 6 reports the iso-surface of the particle VDF at the time $t=T_{nl,1}=\unicode[STIX]{x1D70F}_{1}^{nl}$ (a) and at the time $t=T_{nl,2}=\unicode[STIX]{x1D70F}_{1}^{nl}+\unicode[STIX]{x1D70F}_{2}^{nl}$ . At $t=T_{nl,1}$ , the initial hole-like structure and the slight ring-like signature have been significantly smoothed out. Then, at $t=T_{nl,2}$ , the approach towards the equilibrium is almost complete, the VDF shape becomes almost Maxwellian. Only a slight temperature anisotropy, which is finally thermalized in the late stage of the simulation, is still recovered. The approach towards the equilibrium confirms that small-scale gradients are dissipated quite fast, while the final approach towards the equilibrium – concerning also the thermalization of temperature anisotropy – occurs on larger characteristic times.

Figure 6. Iso-surface of the distribution function, obtained in the case of the fully nonlinear Landau operator. Panels (a) and (b) respectively display the time instants $t=T_{nl,1}=\unicode[STIX]{x1D70F}_{1}^{nl}$ and $t=T_{nl,2}=\unicode[STIX]{x1D70F}_{1}^{nl}+\unicode[STIX]{x1D70F}_{2}^{nl}$ . Red, green and blue axes refer to $v_{x}$ , $v_{y}$ and $v_{z}$ , respectively. The results for the linearized Landau operator are qualitatively the same with the primary difference being that $\unicode[STIX]{x1D70F}_{j}^{nl}$ is substantially smaller than $\unicode[STIX]{x1D70F}_{j}^{lin}$ for each value of $j$ .

In the linearized operator case, two characteristic times are also recovered:

1. (i) $\unicode[STIX]{x1D70F}_{1}^{lin}=0.54\,\unicode[STIX]{x1D708}_{SH}^{-1}\rightarrow \unicode[STIX]{x0394}S_{1}^{nl}/\unicode[STIX]{x0394}S_{tot}=16\,\%$

2. (ii) $\unicode[STIX]{x1D70F}_{2}^{lin}=2.60\,\unicode[STIX]{x1D708}_{SH}^{-1}\rightarrow \unicode[STIX]{x0394}S_{2}^{nl}/\unicode[STIX]{x0394}S_{tot}=84\,\%$

As described for the first case study, these recovered characteristic times are systematically larger (approximately three times) compared to the times recovered in the fully nonlinear operator case. The amount of entropy growth associated with each characteristic time is also different, a smaller amount of entropy growth is indeed associated with the fastest characteristic time when nonlinearities are neglected. The results here described confirm the insights described in the previous section. The evolution obtained with the two operators is qualitatively similar: in both cases, several characteristic times are recovered in the entropy growth and these characteristic times are associated with the dissipation of different velocity space structures. However, the observed evolutions are different from a quantitative point of view: the recovered characteristic times are very different in the two cases, and are significantly larger in the case of the linearized operator.

As described in Paper I, in the first case study, much smaller characteristic times are in general recovered compared to the second case study, probably since the numerical resolution in the second case study is approximately $30$ times lower compared to the first case study and the sharp velocity gradients present in the first case study (figure 1 c) are not accessible in the second case study (figure 4 b). The presence of finer velocity structures in the first case compared to the second case introduces smaller characteristic times.