Solar wind collisional heating

Abstract

To properly describe heating in weakly collisional turbulent plasmas such as the solar wind, interparticle collisions should be taken into account. Collisions can convert ordered energy into heat by means of irreversible relaxation towards the thermal equilibrium. Recently, Pezzi et al. (Phys. Rev. Lett., vol. 116, 2016a, 145001) showed that the plasma collisionality is enhanced by the presence of fine structures in velocity space. Here, the analysis is extended by directly comparing the effects of the fully nonlinear Landau operator and a linearized Landau operator. By focusing on the relaxation towards the equilibrium of an out of equilibrium distribution function in a homogeneous force-free plasma, here it is pointed out that it is significant to retain nonlinearities in the collisional operator to quantify the importance of collisional effects. Although the presence of several characteristic times associated with the dissipation of different phase space structures is recovered in both the cases of the nonlinear and the linearized operators, the influence of these times is different in the two cases. In the linearized operator case, the recovered characteristic times are systematically larger than in the fully nonlinear operator case, this suggesting that fine velocity structures are dissipated more slowly if nonlinearities are neglected in the collisional operator.

1 Introduction

Since the beginning of the last century, many theoretical efforts have been performed to model natural and laboratory plasmas. One of the first attempts to describe the interplanetary medium and its interaction with planetary magnetospheres was conducted by Chapman & Ferraro (1930, 1931), widely considered the fathers of the magnetohydrodynamics (MHD) theory. Their main intuition was to treat plasmas, approximated by neutral conducting fluids, as self-consistent media. One of the basic assumptions of this framework is that interparticle collisions are sufficiently strong to maintain a local thermodynamical equilibrium, e.g. the particle velocity distribution function (VDF) is close to the equilibrium Maxwellian shape. This approach is still widely adopted to analyse plasma dynamics at large scales and many models have been developed to study the features of MHD turbulence (Elsässer 1950; Chandrasekhar 1956; Iroshnikov 1964; Kraichnan 1965; Moffatt 1978; Parker 1979; Dobrowolny, Mangeney & Veltri 1980a , b ; Ng & Bhattacharjee 1996; Matthaeus et al. 1999; Verdini, Velli & Buchlin 2009; Bruno & Carbone 2013; Howes & Nielson 2013; Pezzi et al. 2017a , b ). One of the most studied natural plasmas is the solar wind, which is the high temperature, low density, supersonic flow emitted from the solar atmosphere. The solar wind is a strongly turbulent flow: the typical Reynolds number is approximately $Re\approx 10^{5}$ (Matthaeus et al. 2005); fluctuations are broadband and often exhibit a power-law spectra; several indicators of intermittency are routinely observed (Bruno & Carbone 2013; Matthaeus et al. 2015). Despite the fact that the solar wind is usually approached in terms of MHD turbulence, spacecraft in situ measurements reveal highly complex features, which go beyond the fluid MHD approach. Once the energy is transferred by turbulence towards smaller scales close to the ion inertial scales, kinetic physics signatures are often observed (Sahraoui, Galtier & Belmont 2007; Alexandrova et al. 2008; Gary, Saito & Narita 2010; Bruno & Carbone 2013). The particle VDF often displays a distorted out-of-equilibrium shape characterized by the presence of non-Maxwellian features such as temperature anisotropies, particle beams along the local magnetic field direction and ring-like structures (Marsch 2006; Kasper, Lazarus & Gary 2008; Maruca, Kasper & Bale 2011; Maruca et al. 2013; He et al. 2015). The principal models to take into account kinetic effects are based on the assumption that the plasma is collisionless, e.g. collisions are far too weak to produce any significant effect on the plasma dynamics (Valentini et al. 2005; Daughton et al. 2009; Parashar et al. 2009; Camporeale & Burgess 2011; Valentini et al. 2011a ; Valentini, Perrone & Veltri 2011b ; Greco et al. 2012; Perrone et al. 2012; Servidio et al. 2012; Valentini et al. 2014; Franci et al. 2015; Servidio et al. 2015; Valentini et al. 2016).

We would point out that, in order to comprehend the heating mechanisms of the solar wind, collisional effects should be considered. Indeed collisions are the unique mechanism able to produce irreversible heating from a thermodynamic point of view. Furthermore, to show that collisions can be neglected, the shape of the particle VDF is usually assumed to be close to the equilibrium Maxwellian (Spitzer 1956; Hernandez & Marsch 1985; Maruca et al. 2013). This approximation may be problematic for weakly collisional turbulent plasmas, where kinetic physics strongly distorts the particle VDFs and produces fine structure in velocity space. Collisional effects, which explicitly depend on gradients in velocity space, may be enhanced by the presence of these small-scale structures in velocity space (Pezzi, Valentini & Veltri 2016a ) (hereafter Paper I). Indeed, in Paper I we showed that the collisional thermalization of fine velocity structures occurs on much smaller times with respect to the usual Spitzer–Harm time (Spitzer 1956) $\unicode[STIX]{x1D708}_{SH}^{-1}$ (where $\unicode[STIX]{x1D708}_{SH}\simeq 8\times (0.714\unicode[STIX]{x03C0}ne^{4}\ln \unicode[STIX]{x1D6EC})/(m^{0.5}(3k_{B}T)^{3/2})$ , and $n$ , $e$ , $\ln \unicode[STIX]{x1D6EC}$ , $m$ , $k_{B}$ and $T$ are respectively the particle number density, the unit electric charge, the Coulomb logarithm, the Boltzmann constant and the plasma temperature). The smallest characteristic times may be comparable with the characteristic times of other physical processes. Therefore, collisions could play a significant role into the dissipation of strong gradients in the VDF, thus contributing to the plasma heating.

In this paper we focus on the importance of retaining nonlinearities in the collisional operators. In particular, by means of numerical simulations of a homogeneous force-free plasma, we describe the collisional relaxation towards the equilibrium of an initial VDF which exhibits strong non-Maxwellian signatures. Collisions among particles of the same species are here modelled through the fully nonlinear Landau operator and a linearized Landau operator. A detailed comparison concerning the effects of the two operators indicates that retaining nonlinearities in the collisional integral is crucial to assign proper importance to collisional effects. Indeed, both operators are able to highlight the presence of several characteristic times associated with the dissipation of fine velocity structures. However, the magnitude of these times is different if nonlinearities are neglected: in the linearized operator case, the characteristic times are systematically larger compared to the case of the fully nonlinear operator. This indicates that, when nonlinearities are not taken into account in the mathematical form of the collisional operator, fine velocity structures are dissipated much slower. Results herein described support the idea that to properly quantify the enhancement of collisional effects and, hence, to correctly compare collisional times with other dynamical times, it is important to adopt nonlinear collisional operators.

We would remark that, since the Landau operator is demanding from a computational perspective, self-consistent high-resolution simulations cannot currently be afforded and we are forced to restrict ourselves to the case of a force-free homogeneous plasma, where both force and advection terms have been neglected. This approximation represents a caveat of the work here presented and future studies will be devoted to the generalization of the results here shown to the self-consistent case.

The paper is organized as follows: in § 2 the solar wind heating problem is revisited in order to address and motivate our work. Then, in § 3 we give a brief description of the numerical codes and the adopted methods of analysis. Numerical results of our simulations are also reported and discussed in detail. Finally, in § 4 we conclude and summarize.

2 Solar wind heating: a huge problem

As introduced above, the solar wind is a weakly collisional, strongly turbulent medium (Bruno & Carbone 2013). 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. 1982; Goldstein et al. 1996; Marino et al. 2008; Cranmer et al. 2009). 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 (2013) 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. 2007; Marino et al. 2008; Sahraoui et al. 2009), 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 1995). 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 1956) 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 2014). 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 2006; Servidio et al. 2015), 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 1936), 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 1956), 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. 2016).

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. 2008; Chust et al. 2009; Servidio et al. 2012, 2015). 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 2002; Bobylev & Potapenko 2013; Pezzi et al. 2013, 2014; Escande, Elskens & Doveil 2015; Pezzi, Valentini & Veltri 2015b ; Pezzi et al. 2016a ; Banón Navarro et al. 2016; Hirvijoki et al. 2016; Tigik et al. 2016), 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 1960; Lenard 1960) or the Landau operator (Landau 1936; Akhiezer et al. 1986). 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 2002). 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 1936; Bobylev & Potapenko 2013).

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 1954; Livi & Marsch 1986) and the Dougherty operators (Dougherty 1964; Dougherty & Watson 1967; Pezzi, Valentini & Veltri 2015a ), 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 2007a , b ; Pezzi et al. 2013; Pezzi, Camporeale & Valentini 2016b ).

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 1956) 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. 2015a , 2016a ). 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 1970; Pezzi et al. 2016a ) 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 1991; Kabantsev, Valentini & Driscoll 2006; Valentini, ONeil & Dubin 2006; Anderegg et al. 2009a , b ; Johnston et al. 2009; Valentini et al. 2012) – 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 1977) and can be also observed in the Earth’s magnetosphere (Tokar & Gary 1984; Lu, Wang & Wang 2005). 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 1957) 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. 2012; Valentini et al. 2014; Servidio et al. 2015). 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. 2007, 2014). 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 2006). 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. 2015), 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. 1970; Pezzi et al. 2016a ) 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.

4 Conclusion

To summarize, here we discussed in detail the importance of considering collisions in the description of the weakly collisional plasmas. Collisions are enhanced by the presence of fine velocity space structures, such as the ones naturally generated by kinetic turbulence in the solar wind; therefore, they could play a role into the conversion of the VDF’s free energy into heat, by means of irreversible processes.

In particular, we focused on the importance of retaining nonlinearities in the collisional operator by performing a comparative analysis of the collisional relaxation of an out-of-equilibrium initial VDF. Collisions have been modelled by means of two collisional operators: the fully nonlinear Landau operator and a linearized Landau operator. Due to the demanding computational cost of the collisional integral, we restricted ourselves to the collisional relaxation in a force-free homogeneous plasma. Our results must be clearly extended to the more general, self-consistent case; however, performing a high-resolution collisional simulation cannot currently be afforded.

The cases of study here analysed indicate that both nonlinear and linearized collisional operators are able to detect the presence of several time scales associated with the collisional dissipation of small velocity scales in the particle VDF. A possible explanation of this behaviour is that the linearized operator also involves gradients in its structure while it does not describe the ‘second-order’ gradients related to the Fokker–Planck coefficients of the operator; therefore it is able to recover the presence of several characteristic times. The general message given in Paper I, namely that the presence of sharp velocity space gradients speeds up the entropy growth of the system, is confirmed also in the case of the linearized operator: indeed, the fastest recovered characteristic times are significantly smaller than the common Spitzer–Harm collisional time (Spitzer 1956).

However, we would point out that the importance of the fine velocity structures is weakened if nonlinearities are ignored in the collisional operator. In the case of a linearized collisional operator, slower characteristic times are systematically recovered with respect to the nonlinear operator case. This indicates that, when one neglects the nonlinearities of the collisional integral, fine velocity structures are dissipated slower. Therefore, to properly address the role of collisions and to attribute them the correct relevance with respect to other physical processes (Gary 1993; Matthaeus et al. 2014; Tigik et al. 2016), nonlinearities should be explicitly considered.