Differential kinetic dynamics and heating of ions in the turbulent solar wind

The solar wind plasma is a fully ionized and turbulent gas ejected by the outer layers of the solar corona at very high speed, mainly composed by protons and electrons, with a small percentage of helium nuclei and a significantly lower abundance of heavier ions. Since particle collisions are practically negligible, the solar wind is typically not in a state of thermodynamic equilibrium. Such a complex system must be described through self-consistent and fully nonlinear models, taking into account its multi-species composition and turbulence. We use a kinetic hybrid Vlasov-Maxwell numerical code to reproduce the turbulent energy cascade down to ion kinetic scales, in typical conditions of the uncontaminated solar wind plasma, with the aim of exploring the differential kinetic dynamics of the dominant ion species, namely protons and alpha particles. We show that the response of different species to the fluctuating electromagnetic fields is different. In particular, a significant differential heating of alphas with respect to protons is observed. Interestingly, the preferential heating process occurs in spatial regions nearby the peaks of ion vorticity and where strong deviations from thermodynamic equilibrium are recovered. Moreover, by feeding a simulator of a top-hat ion spectrometer with the output of the kinetic simulations, we show that measurements by such spectrometer planned on board the Turbulence Heating ObserveR (THOR mission), a candidate for the next M4 space mission of the European Space Agency, can provide detailed three-dimensional ion velocity distributions, highlighting important non-Maxwellian features. These results support the idea that future space missions will allow a deeper understanding of the physics of the interplanetary medium.


Introduction
The solar wind is a rarefied, highly variable flow of charged particles originating from the Sun. This medium is composed mainly by electrons and protons, with a small fraction of alpha particles (helium nuclei) and a lower percentage of heavier ions. Due to its low particle density, collisions are essentially negligible, and the solar wind plasma is in a state of non-thermodynamic equilibrium [1]. One of the most puzzling aspects of the solar wind dynamics consists in the empirical evidence [2] that it is hotter than expected from adiabatic expansion [3,2,4]. Understanding the mechanisms of energy dissipation and particle heating in such a collisionless system represents a real challenge for space plasma physics.
'In situ' spacecraft measurements reveal that the solar wind is in a state of fullydeveloped turbulence [5,6]. The energy injected by large scale solar dynamical features into the Heliosphere, in the form of long-wavelength fluctuations, cascades towards small scales via nonlinear interactions until it can be transferred to the plasma as heat. In the inertial range, the power spectrum of the solar wind fluctuations manifests a behavior reminiscent of the Kolmogorov phenomenology for fluid turbulence [7,8,9,10,11]. The turbulent cascade extends to smaller spatial scales, down to a range of wavelengths where kinetic effects start playing a non-negligible role. At the typical ion characteristic scales (Larmor radius and/or inertial length), different physical processes come into play, leading to changes in the spectral shape [12,13,14,15,16,17,18,19]. Here, the dynamics of the solar wind ions gets extremely complicated, being dominated by complex kinetic processes such as resonant wave-particle interactions, particle heating, generation of temperature anisotropy, production of beams of accelerated particles, etc. In this range of scales the particle velocity distribution is generally observed to deviate from the typical Maxwellian shape of thermodynamic equilibrium.
Several models have been developed to understand the kinetic scale phenomenology, focusing specifically on heating and acceleration processes. To explain the above phenomena, several physical mechanisms have been proposed. Among them, the ioncyclotron resonance can produce heating [20,21,22,23,24,25,26,27], where ions interact resonantly with ion-cyclotron waves generated along the cascade. An alternative model is the non-resonant stochastic heating [28,29], where the fluctuations at scales comparable with the ion gyro-scales produce significant distortions of the ion orbits that become stochastic in the plane perpendicular to the magnetic field, violating the magnetic moment conservation. Finally, it has been clearly observed that the interaction of particles with coherent structures can locally produce heating and non-Maxwellian features [30,31,32,33,34,35,36].
One of the most relevant aspects of the ion heating process in the solar wind is that the heavy ions are preferentially heated and accelerated with respect to protons. Measurements of solar wind hydrogen and helium temperatures [37] show compelling evidence of Alfvén-cyclotron dissipation mechanisms. A more recent observational study [38] shows that for the heavy ion component (A > 4 amu) the temperature displays a clear dependence on mass, probably reflecting the physical conditions in the solar corona. Again, the interpretation of these observations is mostly based on the physical processes of wave-particle interactions [39,40,41,42,43,44,45,46,47,48], stochastic ion heating [49], and interaction between particles and coherent structures [50,51]. Despite a significant theoretical effort, there is still no definitive solution to the problem of heavy ions differential heating. It is important to note that the physical processes of the differential heating of heavy ions are also important for laboratory plasmas research, aiming at heating a confined plasma to initiate fusion reactions.
In this paper, we study the kinetic dynamics of protons and alpha particles in typical conditions of the interplanetary medium, by employing multi-component Eulerian hybrid Vlasov-Maxwell (HVM) simulations [53,40]. The multicomponent version of the HVM code integrates the Vlasov equation numerically for the distribution function of both protons and alpha particles, while treating electrons as a massless, isothermal fluid. We present the numerical results of HVM simulations of decaying turbulence with guide field, in a 2D-3V phase space domain (two dimensions in physical space and three dimensions in velocity space). During the turbulent cascade, the ion species depart locally from the condition of thermodynamic equilibrium [33,50,54,55,56], exhibiting temperature anisotropy, differential kinetic behavior as well as preferential heavy ion heating. Interestingly, the preferential heating process occurs in thin filaments of the typical size of few inertial lengths, preferentially located in the proximity of vorticity structures. In the same regions, significant deviations from Maxwellian distributions are observed for ion species. Moreover, as we will discuss in detail below, we have analyzed the distribution of the numerical data in the fire-hose and mirror stability plane, as dependent on the multi-ion composition of the interplanetary medium, showing a good agreement between the numerical results and the observational evidences discussed in a recent paper by Chen et al. [57].
Finally, through a top-hat simulator [58] using the output of the HVM simulations, we show that the Cold Solar Wind (CSW) instrument [59] on board the Turbulence Heating ObserveR -THOR spacecraft [60], a candidate for the next M4 space mission of the European Space Agency, is able to provide adequate high resolution measurements of the ion velocity distributions for the study of turbulent ion heating. In this regard, the THOR mission will allow to measure the velocity distributions of the solar wind protons and alpha particles, with unprecedented phase space resolution.  Semi-logarithmic plot of the ratio between electric and magnetic energy. The vertical black dashed lines indicate the proton skin depth characteristic wavenumber. The Kolmogorov expectation, k −5/3 , is reported as a black dashed line, while the green diamond curve represents the contribution |ηJ k | 2 to the electric energy, due to the resistive term in the Ohm's law in Eq. (2).

Setup of multi-ion hybrid Vlasov-Maxwell simulations
To perform a numerical study of the microphysics of a multi-ion turbulent plasma, we employ an Eulerian hybrid Vlasov-Maxwell (HVM) numerical code [53,50] which solves the Vlasov equations for the distribution functions of protons and alpha particles. In this model electrons are treated as an isothermal massless fluid. Vlasov equations for ions are coupled to a generalized Ohm's law, in which both the Hall term and the electron pressure contribution are retained, and to the Ampére and Faraday equations. The displacement current is neglected and quasi-neutrality is assumed. The dimensionless HVM equations can be written as: In the above equations, masses and charges have been scaled by the proton mass m p and charge e respectively, time by the inverse proton-cyclotron frequency (Ω −1 cp = m p c/eB 0 , B 0 being the ambient magnetic field and c the speed of light), velocity by the Alfvén speed (V A = B 0 / √ 4πn 0,p m p , n 0,p being the equilibrium proton density), in which we only considered the contribution of the dominant proton species, and length by the proton skin depth (d p = V A /Ω cp ). In Eq.
(1), f s (r, v, t) is the ion distribution function (the subscript s = p, α refers to protons and alpha particles, respectively), E(r, t) and B(r, t) are the electric and magnetic fields, and ξ s is a constant that depends on the charge to mass ratio of each ion species (ξ p = 1 and ξ α = 1/2). In Eq. (2), the electron bulk velocity u e is defined as (Σ s Z s n s u s − j)/n e , where Z s is the ion charge number (Z p = 1 and Z α = 2), the density n s and the ion bulk velocity u s are zero-th and first order velocity moment of the ion distribution function, respectively. The electron density is derived from the quasi-neutrality equation n e = Σ s Z s n s . Furthermore, j = ∇ × B is the total current density and P e = n e T e is the electron pressure (T e is assumed to be constant in time and space). A small resistive term has been added as a standard numerical Laplacian dissipation (η = 2 × 10 −2 ), in order to remove any spurious numerical effects due to the generation of strong magnetic field gradients during the development of turbulence.
We solve the multi-ion hybrid Vlasov-Maxwell equations in a 2D-3V phase space domain. The system size in the spatial domain is L = 2π × 20d p in both x and y directions, where periodic boundary conditions have been implemented, while the limits of the velocity domain are fixed at v max,s = ±5v th,s in each velocity direction, v th,s = (k B T 0,s /m s ) 1/2 being the ion thermal speed at equilibrium and T 0,s the equilibrium temperature. In each direction in the velocity domain we set f s (|v| > v max,s ) = 0. The 5D numerical box has been discretized by 512 2 grid-points in the 2D spatial domain and 71 3 grid-points in both proton and alpha-particle 3D velocity domains. A higher velocity resolution has been adopted here with respect to previous works [50,51], in order to ensure a significantly improved conservation of the Vlasov invariants and to provide a better description of the details of the ion distribution function in velocity space. The time step, ∆t, has been chosen in such a way that the Courant-Friedrichs-Lewy condition for the numerical stability of the Vlasov algorithm is satisfied [61]. In order to control numerical accuracy, a set of conservation laws is monitored during the simulations: typical relative variations of the total energy is ≃ 0.6%, of the entropy is ≃ 0.4% and of the mass is ≃ 2.5 × 10 −3 % for protons and ≃ 0.35% for alpha particles.
At t = 0, both ion species have Maxwellian velocity distributions and homogeneous constant densities. To simulate physical conditions close to those of the pristine solar wind, at equilibrium we set the alpha particle to proton density ratio n 0,α /n 0,p = 5% and temperature of each species such that T 0,α = T 0,p = T e . The plasma is embedded in a uniform background out of plane magnetic field, B 0 (B 0 = B 0 e z ). We perturb the initial equilibrium with a 2D spectrum of Fourier modes, for both proton bulk velocity and magnetic field. The energy is injected with random phases and wave numbers in the range 0.1 < k < 0.3, where k = 2πm/L, with 2 ≤ m ≤ 6. Neither density disturbances nor parallel variances are imposed on the initial equilibrium, namely δn = δu z = δB z = 0. We performed three different simulations whose relevant parameters are summarized in Table 1. Below, we will discuss the results of RUN 1. RUN 2 and RUN 3 are qualitatively similar to RUN  significant statistics to perform a direct comparison with solar wind data.

Turbulence cascade and generation of temperature anisotropy
We numerically investigate the kinetic evolution of protons and alpha particles in a situation of decaying turbulence. The large scale fluctuations imposed on the initial equilibrium produce a turbulent cascade down to kinetic scales. In analogy with fluid models, by looking at the time evolution of the spatially averaged mean squared out of plane current density j 2 z , it is possible to identify an instant of time t * at which the turbulent activity reaches its maximum level [62]. The value of t * in units of Ω −1 cp is reported in the right column of Table 1, for each simulation. At this time t * we perform our analysis [33,50,54,56].
In panel (a) of Figure 1, we show the omni-directional magnetic |B k | 2 (blue solidsquare line) and electric |E k | 2 (red dashed-triangle line) spectra at t = t * for RUN 1. The Kolmogorov expectation, k −5/3 (black dashed line), is reported as a reference, while the vertical black dashed line indicates the proton skin depth characteristic wavenumber. In the range of small wavenumbers the magnetic activity is dominant, while the situation clearly changes for kd p > 2, where the electric energy becomes larger than the magnetic one. To highlight this point, in panel (b) of Figure 1 we report the electric to magnetic energy ratio as a function of the wavenumber. The break point between the two regimes is located around the proton skin depth, where dispersive Hall effects come into play [54]. These results present several similarities to spacecraft measurements in the solar wind [15,63,64].
As already noticed in previous works [33,50,54,56,51], ion temperature anisotropy is generated along the development of the turbulent cascade. We focus, at first, on the temperature anisotropy with respect to the direction of the local magnetic field T ⊥ /T for each ion species. The temperature of each ion species has been computed as the second-order velocity moment of the ion distribution function: In Figure 2, we show the 2D contours of the temperature anisotropy of protons (left) and alpha particles (right) at t = t * , for RUN 1 (the same behavior is observed for RUN 2 and RUN 3); as can be noticed from the comparison of both panels of this figure, the global bi-dimensional patterns of the temperature anisotropy display similar features, with larger values of anisotropy being concentrated in thin filaments (of the typical size of few ion skin depths) for both species. However, a significant differential behavior is observed, as the alpha particles appear consistently more anisotropic than the protons.
With the aim of showing that the numerical results of the ensemble of simulations in Table 1 reproduce closely the typical behavior observed in the solar wind multi-ion plasma, we consider here a recent paper by Chen et al. [57], in which solar wind ion temperature anisotropy from the Wind spacecraft at 1 AU has been analyzed. The authors proposed a detailed study of the plasma stability, as dependent on the simultaneous presence of protons, electrons and alpha particles, as main constituents of the solar wind plasma. Specifically, they made a comprehensive analysis based on three years of data and gave clear evidence that the fire-hose and mirror instability thresholds, calculated in Refs. [65,66] for a multi-species plasma system, well constrain the whole data distribution to the stable side and that the contours of the distribution follow the shape of the thresholds.
In order to compare with these observations, we have considered the conditions for the long-wavelength fire-hose and mirror instabilities in a multi-species plasma [65,66]: where s stands for protons, alphas and electrons respectively (s = p, α, e), β ⊥, = s β ⊥, s , ρ s is the mass density of the species s, ρ is the total mass density, v A is the local Alfvén speed, ∆v s is the difference between the bulk velocity of the species s and the center of mass velocity. In Figure 3 [panel (a)] we report a scatter plot of the first and second term of Λ f (logarithmic axes), for all the cases in Table 1.  Table 1. Solid (black) line represent the instability threshold. Panels (b) and (c): log 10 Λ f and log 10 Λ m , respectively, as a function of log 10 β.
Data are clearly constrained in the stable region (delimited by the solid line). These numerical results are in agreement with the solar wind analysis, as can be seen comparing our study with Figure 2 of Ref. [57]. This further confirms that the HVM model can describe the basic dynamics of the solar wind ion species. However, in this direct comparison one has to take into account the limitations of the approximation of isothermal (massless fluid) electrons, and limitations in the statistical convergence (data produced by three numerical runs cannot fully recover the statistics of three years of spacecraft measurements). In panels (b) and (c) of Figure 3, we show log 10 Λ f and log 10 Λ m respectively, as a function of log 10 β, β being the total plasma beta. Here, one can easily see that the instability thresholds (black-dashed lines) are reached for β ≥ 1,

Coherent structures, kinetic effects and differential heating
During the development of the turbulence cascade, temperature increase is observed for both ion species, alpha particles being preferentially heated with respect to protons. Moreover, the generation of ion temperature anisotropy discussed in the previous section is associated to the appearance of coherent structures, such as current sheets, filaments, strongly sheared flows, etc., observed in the 2D patterns in physical space. In particular, panel (a) of Figure 4 displays the shaded contours of the out of plane current density j z at t = t * for RUN 1, together with the contour lines of the magnetic potential A z of the in-plane magnetic field (B ⊥ = ∇A z × e z ). Magnetic flux tubes with different  Figure 5. Spatial profiles of the squared current density j 2 z (red dashed line, panels (a) and (b)) and of the squared vorticity ω 2 z for protons (green dashed line, panels (c) and (d)) and for alpha particles (green dot-dashed line, panels (c) and (d)). In all panels, the alpha particle to proton temperature ratio is also shown (black solid line). All profiles are evaluated along the y direction at two different fixed spatial points: x ∼ 110d p (panels (a) and (c)) and x ∼ 50d p (panels (b) and (d)).
polarizations (or magnetic islands in 2D) are identified by closed contours of A z , as it can be seen from panel (a) of this figure. As it is clear from this contour plot, current density becomes very intense in between adjacent magnetic islands. We point out that the mirror and fire hose unstable points in Figs. 3 (b)-(c) (obtained in RUN 3, with increased initial magnetic perturbations) are located in spatial positions very close to strong current sheets. In panel (b) of Figure 4 we show the contour map of the proton to alpha particle temperature ratio. Even though at t = 0 ions have been initialized with equal temperature (T 0,p = T 0,α ), during the evolution of turbulence, the alpha particles are heated preferentially with respect to protons. Moreover, the 2D map shows that the differential heating is not uniform, but is strongly inhomogeneous, displaying a pattern similar to the 2D map of j z and being significantly influenced by the topology of the magnetic field [51,33,54,56,35]. Moreover, the physical mechanism responsible for the ion heating process is evidently more efficient for the alpha particles than for the protons.
As the HVM code allows for a clean and almost noise-free description of the ion distribution function in phase space, we can analyze the kinetic dynamics of the ion species associated to the development of the turbulent cascade and to the differential ion heating. Several 'in situ' observations (see, for instance, Refs. [1,67]  that the ion distribution functions in the turbulent solar wind (especially in the less collisional fast wind) are typically far from thermodynamical equilibrium and that kinetic effects manifest through the appearance of field-aligned beams of accelerated particles, generation of ring-like modulations in the particle velocity distributions and, in general, through complicated non-Maxwellian deformations. In order to quantify the deviation of the ion velocity distributions from the Maxwellian shape in our simulations, we introduce the following non-Maxwellian measure [34], for each ion species: where g s is the associated equivalent Maxwellian distribution computed from the parameters of f s [34]. We point out that at t = 0, ǫ s is null, since both ion species have been initialized with a Maxwellian velocity distribution. We observe that, simultaneously with the evolution of turbulence, the ion distributions dramatically deviate from the Maxwellian shape. Panels (c) and (d) of Figure 4 display the contour maps of ǫ s for protons and alpha particles, respectively, at t = t * for RUN1. As it is clear from this figure, both ion species significantly deviate from thermodynamic equilibrium and the spatial distribution of ǫ s is not uniform, displaying features quite similar to the 2D contours of the out of plane current density [panel (a)] and of the ion temperature ratio [panel (b)]. In particular, kinetic effects seem to be localized near the peaks of the current density in the shape of thin filaments. It is worth noting that deviations from Maxwellian are significantly stronger for the alpha particles than for the protons; that is, larger values of ǫ α are recovered with respect to ǫ p and are achieved in spatial positions close to peaks of the current density, where differential heating is also observed. Ion differential heating usually occurs close to thin current sheets. Panels (a) and (b) of Figure 5 report the y profile of j 2 z (red-dashed curves), in two different ranges y ∈ [3, 7]d p and y ∈ [116, 122]d p , computed at x ∼ 50d p and x ∼ 110d p , respectively. Two current sheets of typical width of a few proton skin depths are visible here. The black-solid curves in the same panels represent the corresponding y profile of T α /T p . As can be seen in the two panels, differential heating is observed close to each strong current sheet. However, the behavior of the alpha to proton temperature ratio looks different in the two cases: in panel (a), alpha particles are preferentially heated in the center of the current sheet, while in panel (b) the differential heating occurs on both edges of the current sheet. This is a quite general behavior, in the sense that, close to each coherent current structure identified in the spatial domain, the alpha to proton temperature ratio displays a peak in the middle of the current sheet (as in panel (a)) or two peaks at its boundaries (as in panel (b)). Panels (c) and (d) of the same figure show similar plots for the squared vorticity component z associated to the proton velocity, ω 2 z,p , and to the alpha particle velocity, ω 2 z,α where ω = ∇ × v, at the same two positions as for the current j 2 z . The correlation observed for the vertical current is also clearly present for the vorticity.
The robustness of the correlation observed between the quantities shown in Figure 5 can be visualized by means of the joint probability distributions P (T α /T p , j z ) and similarly for the vorticity P (T α /T p , ω z ), as shown in Figure 6. In the left panel, it is evident that there is actually no statistical correlation between the current and the differential heating, because of the spatial shift discussed above. On the contrary, an higher correlation is visible between the proton vorticity and the differential heating (central panel of Figure 6). Even a stronger correlation is observed with the alphaparticle vorticity (right panel). A more quantitative assessment can be obtained by estimating the Spearman correlation coefficients between the alpha-proton temperature ratio and the current, C S (T α /T p , j z ), and similarly for the two vorticity fields, as reported in each panel of Figure 6. This observation suggests that the vortical motion of the alpha particles may play a significant role in, or may be significantly affected by, the differential heating process [70]. This correlation can be understood also in terms of magnetic reconnection. Topologically, magnetic reconnection might form quadrupolar vortical structures located near the current peak, as clearly shown in [71]. In these regions, both magnetic and velocity shears can trigger non-Maxwellian and thermal processes. A deeper analysis of the relationship of current and vortex structures with collisionless plasma heating can be found in Ref. [72].
At this point, it is important to better quantify the connection between the differential particle heating and the level of kinetic activity. In Figure 7 we show the joint distributions of ǫ α and ǫ p (left panel) and of ǫ α and T α /T p (right panel) for RUN 1 at t = t * . Both panels show that (i) alpha and proton velocity distributions tend to deviate from Maxwellian at similar spatial locations, and (ii) such a deviation is correlated with the differential heating. In Figure 8, the joint distribution of ǫ p and ǫ α has been conditioned to specific ranges of values of the temperature ratio T α /T p , divided into five bins as indicated in each panel. Each distribution is normalized to its peak value. As the temperature ratio increases (from top to bottom in the figure), an evident shift of the distribution towards larger values of both ǫ p and ǫ α is observed; this shift can be clearly appreciated by looking at the location of the average value of T α /T p , indicated by a red dot in each panel of this figure. These results indicate that a larger deviation from a Maxwellian distribution is observed for larger temperature ratios T α /T p . On the other hand, the conditioned joint distributions are relatively symmetric, so that there is no evidence of different increase in ǫ p or ǫ α when the alpha particles are preferentially heated.
Temperature variations of the ions seem to be related to their kinetic dynamics and, therefore, to the complex deformations of the particle velocity distributions. Future space missions, such as THOR [60], designed to provide measurements of the threedimensional ion distribution functions with very high phase space resolution, will be able to capture the non-Maxwellian deviations of the velocity distributions, so as to provide relevant clues on the long-standing problem of particle heating.
As discussed previously, close to the sites of enhanced magnetic activity, the non-Maxwellian measure ǫ s achieves large values for both ion species. In order to show what the three-dimensional ion velocity distributions look like in turbulence, in Figure 9 we report the iso-surfaces of the proton [panel (a)] and alpha particle [panel (b)] velocity distributions, computed in the spatial point where ǫ s is maximum. The red tube in both panels indicates the direction of the local magnetic field. The proton velocity distribution displays a field-aligned beam and the generation of ring-like modulations in planes perpendicular to the direction of the local field, presumably signatures of cyclotron resonance. On the other hand, for the alpha particles we recovered more evident deformations, with the appearance of multiple blobs and elongations; the alpha particle velocity distribution seems to lose any property of symmetry with respect to the direction of the local magnetic field. It is clear from these plots that kinetic effects, working along the turbulent cascade at typical ion scales, can efficiently drive the plasma away from thermodynamic equilibrium. Moreover, when the particle velocity distributions are not constrained by limiting assumptions on their shape in velocity space, they can be distorted in a very complicated way and lose their symmetry. Figure 9. Iso-surfaces of the proton (panel (a)) and alpha particle (panel (b)) velocity distribution in a spatial point where ǫ s is maximum. The direction of the local magnetic field vector is reported as a red tube.

Deviations from Maxwellian distributions
An alternative way to describe the deviations from thermodynamic equilibrium consists in defining the so-called minimum variance frame of the three-dimensional velocity distribution. We compute the preferred directions of the velocity distribution of each ion species in velocity space [33], in each spatial position, from the stress tensor: This tensor can be diagonalized by computing its eigenvalues {λ 1 , λ 2 , λ 3 } (ordered in such way that λ 1 > λ 2 > λ 3 ) and the corresponding normalized eigenvectors {ê 1 ,ê 2 ,ê 3 }. We point out that λ i are proportional to temperatures andê i are the anisotropy directions of the velocity distribution.
From this latter figure it can be noticed that, in the range of small ǫ s , the PDFs (for both protons and alphas) have a peak close to unity (they are not exactly centered around 1, since the minimum value of ǫ s is not zero at t = t * ). As ǫ s increases, high tails appear in the PDFs suggesting that, in the case of significant deviations from Maxwellian, the ion velocity distribution loses its properties of isotropy and gyrotropy. This observation suggests that the use of reduced models, based on restrictive approximations on the shape of the ion velocity distribution, might not be appropriate, and more complete models, describing the evolution of the velocity distributions in a full 3D velocity space, should be employed.
With the aim of characterizing the nature of the deformation of the particle velocity distributions and to identify the spatial regions which are the sites of kinetic activity, we computed for each ion species two indices of departure from Maxwellian, i. e. the temperature anisotropy index and the gyrotropy index, in two different reference frames, namely the MVF and the local magnetic field frame (LMF). We define the anisotropy measures ζ = |1 − λ 1 /λ 3 | (MVF) and ζ * = |1 − T ⊥ /T | (LMF), where T ⊥ and T are the temperatures with respect to the local magnetic field, and the gyrotropy indicator in the MVF η = |1 − λ 2 /λ 3 |. The gyrotropy index in the LMF η * can be computed by using the normalized Frobenius norm of the nongyrotropic part N of the full pressure tensor Π [73]: where N ij are the components of the tensor N, and T r(N) = 0. It is worth to point out that all indices defined above are identically zero if the particle velocity distribution is Maxwellian. At the maximum of the turbulent activity, these indices significantly differ from zero for both ion species, meaning that the ion velocity distributions are anisotropic and non-gyrotropic in both MVF and LMF.
In Figure 11, we show the spatial profiles of all these indices inside the region where the current sheet shown in Figures 5 (b) is located (indicated by vertical black-dashed lines). The indices for protons are reported in the left panels, while for alpha particles in the right panels. As it is evident from this figure, the non-Maxwellian measures of both protons and alpha particles are not null in the regions close to the peak of the current density, in both reference frames; specifically, the ion velocity distributions become anisotropic and non-gyrotropic, as can be seen from both the MVF and LMF indices. It is worth noting that the indices for the alpha particles are significantly larger than for the protons, additional evidence that alpha particles are more efficiently driven out of equilibrium with respect to protons.

Optimized measurements for the study of protons and alpha particles kinetic dynamics
In the previous sections we have shown that kinetic effects produce highly non Maxwellian distributions, particle beams and ring-like modulations. An instrument designed to study these effects should have proper phase space (and temporal) resolution, since a poor resolution could impede the observations of such fine features. In near-Earth space it is possible to perform high resolution in situ plasma observations. Measurements of particle velocity distribution functions are usually performed by top-hat electrostatic analyzers [74]. Such instruments are composed of two concentric hemispheres set at different voltage, with an aperture on the outer hemisphere. The electric field between the plates allows for particles within a specific energy-per-charge (E/q) ratio to pass through the gap and reach the detector. Varying the voltage permits to sort particles according to their E/q. Parallel incident particles will be focused on a specific sector of the detector at the analyzer exit, each sector identifying the particle velocity azimuthal direction. Analyzers of this type measure incoming particles with a 360 • disk shaped field of view. In order to sample the entire 4π solid angle, either the spacecraft rotation is used or electrostatic polar deflectors are employed. Usually analyzers for solar wind measurements are characterized by a restricted field of view, since the solar wind at 1 AU is a cold beam in velocity space. A schematic illustration of such a top-hat analyzer is given in figure 12.
As already said, the energy-angular resolution of the top-hat analyzer is extremely important when aiming to study the signatures of kinetic-scale processes, since low phase space resolution would smooth out the fine structures of the particles distribution functions. The Cold Solar Wind (CSW) [59] on board the THOR mission is a tophat analyzer conceived to measure solar wind protons and alpha particles distribution functions with unprecedented temporal and phase space resolutions [58,59]. In particular, it has an angular and energy resolution that will permit to observe the fine structure characterizing the particle distribution functions related to the effects of turbulence. In this respect, we use a top-hat simulator to find out how the proton and alpha velocity distributions characterized by large ǫ s , and shown in Figure 9, would be detected by CSW. Such distribution functions are meaningful examples since they are characterized by non Maxwellian shapes and a high degree of complexity as has been shown in the previous sections. Moreover, we compare the moments and the ǫ s of the HVM and CSW distributions in proximity of large ǫ s sites, showing that CSW measurements capture the complexity present in the HVM distributions.
CSW will detect particles in 96 energy-per-charge intervals with ∆E/E = 7%. It will have a restricted field of view of 48 • in both azimuth and elevation, focused on the cold solar wind population, and will have an unprecedentedly high angular resolutions of 1.5 • , in both azimuth and elevation. The CSW response is further determined by its geometric factor G = 2.2 × 10 −5 cm 2 sr and the time needed for the accumulation of particle counts in one energy channel and elevation angle T acc ∼ 0.25ms.  Figure 14. HVM (red lines) and CSW (black lines) density, velocity and temperature, computed for protons (left) and alpha particles (right) distributions along the y and x direction at x 0 ≃ 60d p and y 0 ≃ 115d p , respectively.
As a first step, the HVM distributions are expressed in physical units by rescaling the normalized variables of the simulation to physical quantities. In order to adapt to a real solar wind, a proton density of 8.2 cm −3 , a proton thermal velocity of 35 km/s, which corresponds to a thermal energy of 12.8 eV, and a T α /T p = 4 have been used. Moreover, a bulk velocity of 400 km/s has been taken into account. Then the distributions functions are reassigned to the field of view of CSW by means of an averaging process. At this point, the counts that the sensors would detect in each energy-angular interval for each species are computed [75]: where i, j, k refer to the energy, elevation and azimuth intervals, v i is the velocity in the interval, f the distribution function, T acc and G the accumulation time and the geometric factor, respectively. f is normalized in such a way that its integral in velocity space gives the particle number density. Afterwards, we determine the corresponding moments.
In Figure 13 proton (top panels) and alpha (bottom panels) velocity distributions are shown at the spatial point where ǫ is maximum, as recovered from the HVM simulations (left panels) and as sampled by CSW (right panels). In particular, slices of the distribution functions in the V x V z plane are presented.
Looking at the protons, it can be seen that all the features that characterize the HVM distribution are revealed by CSW. In fact, the ring-like modulations, that can be identified in the HVM distribution along the V x direction, are present also in the CSW distribution, although they are less readily identifiable. Moreover, this distribution function spans only over one order of magnitude and will be easily observed by CSW for regular or dense solar wind.
Looking at the alpha particles, it has to be noted that the alphas population appears in the CSW distribution at √ 2 times the bulk velocity recovered from the simulation, since electrostatic analyzers like CSW detect particles with respect to their energy-percharge. Regarding the features characterizing the HVM and CSW alpha distributions, it can be seen that the shape of the main alpha population, at about (V x , V z ) = (−340, 30), is correctly reproduced. In fact, the elongation visible in the HVM is clearly unveiled also in the CSW distribution, and also discernible is the faint widening on the left part of the core population. Moreover, the faint blob at negative V z in the HVM distribution, is recognizable in the CSW distribution as well, and its shape is well captured.
In Figure 14 we show the HVM (red lines) and CSW (black lines) density, velocity and temperature, for protons (left panel) and alpha particles (right panel), computed for distribution functions along the y and x direction at x 0 ≃ 60d p and y 0 ≃ 115d p , respectively. Along these directions, the maximum ǫ p and ǫ α are encountered. Results show that the velocity, and to a lesser extent the density, for HVM and CSW are extremely well in agreement. Small differences can be noted between HVM and CSW temperature. The different behaviors of the temperature and lower order moments was already observed in [58]. In general, larger errors are found between the HVM and CSW moments for alpha particles than for protons. This can be explained in terms of the alpha particles lower counts statistics with respect to the protons.
Finally, the non Maxwellian measure ǫ s , computed for the HVM (red line line) and CSW (black line) distributions, is presented in Figure 15. The ǫ s evolution is the same for the HVM and CSW distributions. In particular, the two HVM ǫ s peaks, of the protons and alphas cases, are well pronounced in the CSW ǫ s evolution as well. The small difference between the peaks of the HVM and CSW ǫ α can probably be explained in terms of the alphas lower counts statistics, like in the case of moments.
It must be noted here that, starting from the HVM simulation, the moments for the alpha particles can be readily computed. In the real case, although the alpha population can be easily identified in the CSW distribution functions, part of the protons and alpha particles populations can overlap, especially when the plasma temperature is high, and the moments computation can be hard (we recall that CSW selects particles according to their energy-per-charge). In such cases it is extremely useful to have a particle instrument that is able to identify the different ion species. Such an instrument will be on board THOR: it is the Ion Mass Spectrometer (IMS), a top-hat analyzer complemented with a time of flight section.

Summary and conclusions
We have employed a multi-component Eulerian hybrid Vlasov-Maxwell code to analyze the differential kinetic dynamics of the dominant ions in the turbulent solar wind. We have performed simulations of decaying turbulence in a five-dimensional phase space, varying relevant parameters, such as the level of turbulence and the plasma beta at equilibrium, in order to explore a large portion of the parameter space and to compare our numerical results with direct solar wind measurements. As the numerical results do not change qualitatively for the runs summarized in Table 1, after showing that our simulations reproduce with a good degree of realism the phenomenology recovered in solar wind observations, we restricted our analysis to a single simulation.
We focused in particular on the departure of the plasma from thermodynamic equilibrium along the cascade and on the role of kinetic effects coming into play when the cascading energy reaches the typical ion scales. As turbulence develops, we observed the generation of ion temperature anisotropy with respect to the direction of the local magnetic field as well as a significant heating for both ion species. Also, the non-Maxwellian index ǫ s , which is a measure of the deviations from Maxwellian of each ion species, gets different from zero; at the same time, the generation of coherent structures (current sheets, vortices, etc.) of the typical width of few ion skin depths is observed in the contour map of the out of plane current density and of the ion vorticity. Generation of temperature anisotropy, particle heating and deformation of the ion velocity distributions clearly occur in the form of thin filaments at scales of the order of the ion kinetic scales.
A detailed statistical analysis of the simulation results provides evidence that the distribution of the numerical data is constrained to the stable region with respect to the fire-hose and mirror stability thresholds, evaluated by taking into account the multi-ion composition of the solar wind plasma. These results appear in good agreement with the observational evidences, recently discussed in Ref. [57].
The generation of temperature anisotropy and the process of particle heating are significantly more efficient for the alpha particles than for the protons and also the non-Maxwellian index gets larger values for the heavy ions. This means that not only a differential heating occurs, but also a differential kinetic behavior is observed for the dominant ion species in the solar wind: in fact, as shown in Figure 9, the velocity distribution of the alpha particles appears much more distorted than that of the protons. Interestingly, the process of differential ion heating takes place close to the peaks of the ion vorticity and in spatial regions where the ion velocity distributions appear significantly distorted and far from thermodynamic equilibrium.
In order to characterize the deviations from Maxwellian, we analyzed the ion velocity distribution both in its minimum variance frame and in the local magnetic field frame, by defining the anisotropy and the gyrotropy indicators in both frames. The study of these indicators reveals that, when kinetic effects are at work and ǫ s gets large, the shape of the velocity distributions of each ion species loses its symmetry, developing anisotropic and non-gyrotropic features. Moreover, close to the peaks of ǫ s , sharp velocity gradients and complicated distortions are recovered.
If, as suggested by our numerical results, temperature variations of the ion species are related to the kinetic dynamics at short spatial scales, then the details of the velocity distribution, which keep the memory of the particle kinetic dynamics, are clearly crucial ingredients for the understanding of the complex process of particle heating. In fact, as we discussed above, different ion species, even though experiencing the same turbulent fields, exhibit a deeply different response, being heated and shaped differentially. By feeding a virtual top-hat analyzer with the ion velocity distributions obtained from the kinetic simulations discussed here, we have demonstrated that the energy and angular resolutions of the Cold Solar Wind instrument on board the THOR mission should allow to obtain unprecedentedly high resolution measurements of both proton and alpha velocity distributions and will therefore provide important information for unveiling the puzzling aspects of solar wind heating. The results discussed in this paper support the idea that future space missions, designed to provide measurements of particle velocity distributions and fields with very high resolution, will allow to gain relevant insights into the dynamics of the solar wind plasma at kinetic scales.