Skew-Kappa distribution functions&whistler-heat-flux instability in the solar wind: the core-strahlo model

Electron velocity distributions in the solar wind are known to have field-aligned skewness, which has been characterized by the presence of secondary populations such as the halo and strahl. Skewness may provide energy for the excitation of electromagnetic instabilities, such as the whistler heat-flux instability (WHFI), that may play an important role in regulating the electron heat-flux in the solar wind. Here we use kinetic theory to analyze the stability of the WHFI in a solar-wind-like plasma where solar wind core, halo and strahl electrons are described as a superposition of two distributions: a Maxwellian core, and another population modeled by a Kappa distribution to which an asymmetry term has been added, representing the halo and also the strahl. Considering distributions with small skewness we solve the dispersion relation for the parallel propagating whistler-mode and study its linear stability for different plasma parameters. Our results show that the WHFI can develop in this system, and provide stability thresholds for this instability, as a function of the electron beta and the parallel electron heat-flux, to be compared with observational data. However, since different plasma states, with different stability level to the WHFI, can have the same moment heat-flux value, it is the skewness (i.e. the asymmetry of the distribution along the magnetic field), and not the heat-flux, the best indicator of instabilities. Thus, systems with high heat-flux can be stable enough to WHFI, so that it is not clear if the instability can effectively regulate the heat-flux values through wave-particle interactions.


INTRODUCTION
Space plasmas are magnetized systems that can be out of thermal equilibrium due to the low collision frequency between its constituent particles. Coulomb collisions are an efficient mechanism to relax particle populations to thermodynamic equilibrium, where the distribution functions reduce to Maxwellian profiles. Therefore, on collisionless systems the particles velocity distribution function can develop non-thermal features. These non-thermal characteristics represent free energy in the system, that can be emitted as electromagnetic radiation, such that the plasma relaxes to more stable states via non-collisional processes as wave particle interactions.
The solar wind shows several non-thermal features in the electron velocity distribution (eVDF).
Among them, the field-aligned skewness Scudder & Olbert 1979;Marsch et al. 1982;Pilipp et al. 1987;Salem et al. 2003;Marsch et al. 2004;Nieves-Chinchilla & Viñas 2008) is clearly observed. This asymmetry provides the energy to excite different modes, depending on the plasma parameters Shaaban et al. 2018a). One of these modes correspond to the so-called whistler heat-flux instability (WHFI), and the excitation of this branch is due to the free energy provided by the skewness or asymmetry of the eVDF. The WHFI has received attention throughout the years because the associated whistler waves are one of the main candidates for a non-collisional regulation of the electron heat-flux values in the solar wind (Abraham-Shrauner & Feldman 1977a;Gary et al. 1994;Gary & Li 2000;Scime et al. 1994Scime et al. , 2001Lacombe et al. 2014;Kuzichev et al. 2019;Shaaban et al. 2019b).
There is ample observational evidence that the electron heat-flux in the solar wind cannot be fully described by the collisional Spitzer-Harm theory (Spitzer Jr & Härm 1953). This model is able to adequately describe measurements of the heat-flux under slow solar wind conditions, but in many other cases the Spitzer-Harm law predicts higher values than those shown by in-situ measurements at 1 AU from the Sun (Bale et al. 2013). This phenomenon of depletion of the electron heat-flux below the values predicted by the collisional transport model has been studied for decades. For example, empirical calculations have been carried out in order to reproduce the measured heat-flux values through an ad-hoc reduction of the thermal conductivity (Cuperman et al. 1972). Theoretical models have been also proposed, considering different physical mechanisms that can potentially regulate the electron heat-flux through collisionless mechanisms (Forslund 1970;Hollweg & Jokipii 1972;Landi et al. 2012). The most accepted mechanism to explain this suppression corresponds to a non-collisional regulation due to kinetic process of wave-particle interaction Hollweg 1974;Perkins 1973). The main candidate for constraining the electron heat-flux values are the whistler waves excited by the heat-flux instability, however, the dominant wave mode involved in this non collisional regulation process is still under debate (Gary & Feldman 1977;Scime et al. 2001;Bale et al. 2013;Shaaban et al. 2018b;López et al. 2020a).
Most studies of the solar wind electron skewness assumed that the electron populations are composed of different sub-populations; each modelled by a Maxwellian, Bi-Maxwellian or Kappa distribution, which combined can form a skew non-thermal distribution. An example of this in the solar wind is the linear superposition of the drifting electron core, halo and strahl populations (Štverák et al. 2009;Saeed et al. 2016;Lazar et al. 2018;López et al. 2020a). Note that these typical distribution functions are symmetrical by themselves, and do not show any skewness. Under this context, in recent years studies have been developed where less common functions are used to model the electron sub-populations. The particularity of these "new" ad-hoc functions is that they are, in fact, asymmetrical. For example, in Horaites et al. (2018a) the authors analyze the kinetic stability of a plasma where the strahl population is described by an analytic function, which was derived from the collisional kinetic equation. Also, in Vasko et al. (2019) the authors modeled the strahl population by means of a bi-Maxwellian function to which extra parameters were added that allows the modification of its symmetry.
Along the same line, here we propose a new heuristic model for solar wind electrons, that can reproduce the behaviour of a core-halo-strahl representation but using only two sub-populations: a bi-Maxwellian core plus a modification to the Kappa distribution that introduces skewness, representing the halo and strahl electrons in a single skew distribution. This skew-Kappa distribution was first proposed by Beck (2000) in a study of fluid turbulence. In the original derivation, the author showed that the asymmetry of the VDF is related with the level of turbulence of the media, measured by the Reynolds number. The aim of this work is to study, using kinetic linear theory, the effect of nonthermal electrons described by a skew Kappa-like function over the excitation of parallel propagating whistler modes associated with the WHFI in a non-collisional, magnetized, solar wind like plasma.
We will show that the proposed eVDF reproduces the main field-aligned features of the eVDF as observed in the solar wind, potentially allowing simpler models of solar wind electrons modeled as a superposition of two sub-populations. In addition, considering a unified description of halo and strahl electrons may also be relevant for the understanding of the relevance of the electron nonthermal features for the dynamics of the heat transport by the solar wind (Bale et al. 2013), and also the kinetic physics governing the halo formation and its relation with the strahl (Vocks et al. 2005;Horaites et al. 2017Horaites et al. , 2018b.
The paper is organized as follows: in Section 2 we present our model, analyze the skew Kappa-like function and introduce it as a new distribution function for describing the halo and strahl electron populations. In section 3 we show the theoretical results of linear kinetic theory for the dispersion tensor of parallel propagating waves, in a plasma where the electron population is modeled by a core and skew kappa-like function. Then, in section 4 we systematize this analysis in order to obtain the the marginal stability thresholds for this distribution as a function of plasma beta and heat-flux, and present the best fit parameters for these contours. Finally, in Section 5 we summarize and discuss our results.

ELECTRON DISTRIBUTION: THE CORE-STRAHLO MODEL
We model the solar wind electrons distribution function f e as a superposition of two sub populations. (1) A bi-Maxwellian distribution (f c ) representing the core and a skew-Kappa function (f s ) to describe both, the halo and strahl electrons, that from now on we will name as the strahlo, and this representation of solar wind electrons as the core-strahlo model.
Under this model f s consists of a Kappa function to which an asymmetry term has been added. Namely, In both distribution functions, the subscripts and ⊥ are with respect to the background magnetic field, and n c , n s , denote the number density of the core and strahlo, respectively. In Eq.
(2) α ⊥ , α are the thermal speeds, and U c is the drift of the core. Also, in Eq. (3) A s is a normalization term such that f s dv = n s , θ and θ ⊥ are related to the thermal velocities as defined in Eqs. (A8) and (A6), respectively. Also, κ s is a measurement of the deviation of this function from a Maxwellian distribution, and δ s controls the field-aligned skewness. Note that when δ s = 0 we recover the well known Kappa distribution (Olbert 1968;Vasyliunas 1968;Scudder 1996;Maksimovic et al. 2005;Xiao et al. 2006;Lazar et al. 2016Lazar et al. , 2017Viñas et al. 2017).

Validity of the model
Depending on the value of the κ s and δ s parameters, Eq. (3) may become negative, complex or multi-valued, which imposes some caveats and limitations to the use of the skew-Kappa for the electron VDF. In particular, for an arbitrary value of δ s and κ s there is a particular value u = v /θ in which the skew-Kappa distribution diverges following a vertical asymptote. This value corresponds to the real solution of the following equation that always exists for real values of κ s and δ s . The dependency of u is strong and weak with respect to δ s and κ s , respectively. For example, for δ s = 0.1, the real solutions of Eq. (4) are u 30.1 for κ s = 3, and u 30.4 for κ s = 10, respectively. Furthermore, in the case of δ s = 0.2, the values are u 15.3 and u 15.7 (in units of the thermal speed of the strahlo) for κ s = 3 and κ s = 10, respectively. Moreover, due to the peak of the VDF at v 0 and the presence of the mentioned asymptote, the distribution always has a local minima at u min , with 0 < u min < u, given by the solution of the derivative of Eq. (4). Namely, a monotonically decreasing function of δ s , with u min 20.0 for δ s = 0.1, and u min 10.1 for δ s = 0.2.
Therefore, for an arbitrary value of δ s there is a speed regime in which the integrals necessary to build the moments of the VDF or the dispersion relation will present vertical asymptotes, branch cuts and poles so that the analytical continuation of the functions in the complex plane may become a quite complicated task. Even though we believe it may be possible to obtain a bounded reasonable solution, such calculation for any arbitrary parameters is beyond the scope of this article. Therefore, our skew-kappa model requires careful treatment when selecting the δ s values.
To avoid these issues, here we apply the heuristic core-strahlo model to situations in which the VDF has small skewness, and the asympote is far away from the main core in units of the thermal speed, so that we can expand all relevant integrals in a finite Taylor series around δ s = 0. It is important to mention that, since f s has a vertical asymptote at u, such Taylor series is not mathematically possible near v = uθ . For it to be allowable requires the first derivatives of f s with respect to v must exist, which it does not for those velocity values that are solutions of Eq. (4). However, this mathematical problem can be evaded if all relevant features of the distribution are contained at velocities within the |v |/θ < |u min | range, i.e., the asymptote of the VDF is far away from the main core. In this case, even though the Taylor series approximation will not be able to mathematically reproduce the exact VDF for all velocity values, calculations based on the approximated version of the VDF in the whole velocity domain will allow analytic calculations, keeping all of relevant physical properties of the skew-Kappa distribution, which subsequently will lead to a direct interpretation of the results and the relevance of each parameter. On the other hand, the general case with arbitrary skewness, when the asymptote may be closer to the main core of the VDF, remains to be solved. In such case the Taylor expansion approach may not be an adequate representation of the VDF near the singularity, and other functional expressions with more attractive properties in the complex plane could be a better option. Under this context, another way to approximate the initial distribution for arbitrary skewness may be the expansion of Eq. (3) on a different base. After a preliminary analysis it seems that the Padé approximant (Bender & Orszag 1999) may be a reasonable procedure for such endeavour, as this approximation does not present new singularities. We will leave this analysis to a future study. From now on we will consider small values of δ s , such that δ 3 s 1, and we will make use of a Taylor expansion of Eq.
(3) up to order δ 2 s (see Eq. (A1) in Appendix A). Figure 1 shows parallel cuts at v ⊥ = 0 of the electron eVDF considering isotropic sub-populations with n s /n e = 0.1, T s /T c = 7.0, and different choices of κ s and δ s . Top and bottom panels show the skew-Kappa strahlo (given by Eq. (3) and the total electron VDF (Eq. (1)), respectively, comparing the exact distribution with a Taylor expansion of up to order δ 2 s as shown in Eq. (A1) in Appendix A. In addition, vertical dotted lines indicate the value of v min given by Eq. 5. From the figure we can see that within the |v | < |v min | velocity range the exact and approximated versions of the VDF are mostly the same. In particular, all relevant features of the VDF, such as the skewness and suprathermal tails, can be clearly observed in both representations inside the |v | < |v min | velocity range (as we will see in Section 2.2). Thus, as all physical properties of the VDF are covered (shape, moments and dispersion properties), in the small skewness regime (δ 3 s 1), a second order approximation of electrons following a skew-Kappa distribution given by Eq. (3) can be reasonably represented by the Taylor expansion as shown in Eq. (A1). In this case all the dispersion functions are reduced to a superposition of standard integrals of the Kappa distribution in v similar to the Q integral given by Eq.(5) in Mace & Hellberg (1995) or Eq.(12) in Hellberg & Mace (2002), a regime already well investigated for integer (Summers & Thorne 1991) or arbitrary (Mace & Hellberg 1995;Hellberg & Mace 2002) values of the κ s parameter.
To further ascertain the validity of the expansion and the dispersion relation analysis results for the heat-flux instability, we have carried out a comparison of the dispersion properties between a corehalo model based upon drifting Maxwellian distributions and those of the skew-Kappa core-strahlo model, using the same parameters (see Figures 4b and 4c and their discussion in the next Section).
The dispersion results (shown on Figures 4b and 4c) demonstrate that both model essentially reproduce each other quite well. The real frequencies and growths profiles generated by both model are essentially the same.
Finally, it is worth mentioning that, since the Kappa functions behave as a power-law for large values of the velocity, depending on the value of the κ s parameter, the moments of the distribution may be divergent, which imposes restrictions for κ s . In the case of a standard Kappa VDF the pressure is well defined only for κ s > 3/2. In our case, to have real and finite values of the temperature and heat-flux moments, the values of kappa are restricted to κ s > 5/2 (see Appendix A for details). In summary, considering δ 3 s 1 as in the case of this study, up to second order in δ s the eVDF is real and positive for all real values of v , and the integrals in the velocity space share the same poles and branch cuts of Kappa distributions (Mace & Hellberg 1995;Hellberg & Mace 2002). Consequently, all the moments of the eVDF and also the dielectric tensor of the plasma are well defined within such caveats.

Properties of the core-strahlo model in the small skewness approximation
Even though the core-strahlo model has several free parameters, quasineutrality and zero-current conditions in the ions frame set some relationships between them. In particular, if the ions density is given by n p , to ensure quasineutrality we have n e = n c + n s = n p . In other words: Also, due to the particular shape of the skew-Kappa distribution, for δ s = 0 f s always has a fieldaligned drift U s = −δ s θ /4. Thus, the zero-current condition imposes the value of U c to satisfy Therefore, under this description, in the ions frame there will be a relative drift ∆U between core and strahlo populations given by ∆U = δ s θ (n e /n c )/4. Note that this relative drift is purely due to the skewness of the strahlo. When the electron distribution has no skewness (δ s = 0), then f e reduces to a symmetrical distribution with a quasi-thermal core and a non-thermal halo represented by a Kappa distribution (see e.g Pierrard et al. 2001;Nieves-Chinchilla & Viñas 2008). For more details, full expressions for the macroscopic parameters and the normalization constant of the strahlo distribution function can be found in Appendix A. It is important to recall that these neutrality and quasi-neutrality conditions are restricted to small skewness values. However, as we will show in Figs. 2 and 3, our approximation is able to describe thermal and non-thermal electrons in the solar wind, and the parallel cuts of the eVDF have remarkably similar shapes as previously reported using ISEE-1 (see e.g Fig. 1b in Scudder & Olbert 1979) or Wind (see e.g Fig. 6 in Nieves-Chinchilla & Viñas 2008) data.
To our knowledge, Beck (2000) was the first to propose this type of skew distributions in a study of fully developed hydrodynamic turbulent flows of skew flow velocity distributions via non-extensive statistical mechanics. Under this context the author showed that the asymmetry term, which we have denoted as δ s , is proportional to Re −1/2 , where Re is the Reynolds number of the media.
This model has been successfully applied to fit data from a turbulent jet experiment (Beck 2000) and environmental atmospheric turbulence (Rizzo & Rapisarda 2004). Here is noteworthy to mention that in both cases the adjusted velocity data lies between ± 10 thermal speeds as shown in Figure 2 of both studies, and that the obtained skewness parameters are small. In addition, the intrinsic mathematical issues of the skew-Kappa distribution discussed in Section 2.1 can be neglected. Therefore, for both cases the skew-Kappa model represent a useful tool to study the relevance and nature of skew velocity distributions in turbulent flows.
As turbulence is also present in plasma systems, this suggests that these distributions can be more than an ad-hoc function for the solar wind electrons. The δ s parameter can potentially be related to microscopic physical processes that allow the particles to exhibit a skew distributions. We strongly believe that this point of view should be further examined, and more rigorous theoretical works studying the underlying physics that allows particle distributions to present this non-thermal feature in plasma systems should be developed. However, such first principle description is beyond the objective of this paper. Here we focus on accepting and using heuristically this skew-Kappa distribution to describe the skewness and high energy tails of the eVDF within the aforementioned caveats. This choice allows us to model the eVDF with less free parameters, as an alternative to the usual "core-halo-strahl" models (Shaaban et al. 2019a;Štverák et al. 2009;Sarfraz et al. 2016;López et al. 2019), but at the same time allows us to mimic and study the effect of asymmetry and non-Maxwellian features of the electron population in a solar wind-like plasma, and study their effects on the whistler-heat-flux instability excitation. a feature that is more evident in the outermost contours. Also, we can see that for higher δ s values, the more skew is the distribution. In panels 3(c) and 3(d) we show how the distribution changes for three different values of κ s and fixed δ s = 0.15. In panel (c) we see that when we increase κ s , the high energy tails diminish. This feature is inherited from Kappa distributions, which are reduced to Maxwellian functions in the limit κ → ∞. However, unlike Kappas, skew-Kappa distributions never reduced to Maxwellian distributions because they maintain the skewness for all kappa values when δ s = 0. In panel 3(d) we see that the outer contour seem to shrink proportionally as κ s decreases Bottom panels consider fixed skewness (δ s = 0.15), and different kappa values: κ = 3 (blue), κ = 5 (green), and κ = 9 (red). All other parameters are the same as in Figure 2. In panels (b) and (d), from innermost to outermost, the levels plotted corresponds to f (v) = 5 × 10 −1 , 2 × 10 −3 , 3 × 10 −5 and 2 × 10 −6 , for all parameter combinations. In all panels parallel and perpendicular velocity components are expressed in unit of the thermal speed of the core (α ⊥ = α ).
while the core does not change, so that the overall contours' shape appears to remain the same.
Therefore, this feature suggests that κ s does not alter the distribution symmetry.
In summary, the combination of Maxwellian and skew-Kappa distributions can model three important non-thermal features observed in the eVDF in the solar wind, namely, quasi-thermal core, enhanced tails, and skewness. Therefore, its use may allow simpler solar wind models, where electrons are modeled as the superposition of core and strahlo, where the distribution skewness is controlled by only one parameter. This field-aligned skewness provides the energy for the excitation of the WHFI, on which we focus the analysis in this work.

LINEAR THEORY AND DISPERSION RELATION
We use linear kinetic theory to derive the dispersion relation of wave modes that can propagate in a magnetized, non-collisional, and initially uniform plasma. We perform this calculation in order to analyze the stability of the whistler mode associated with WHFI in a solar wind like plasma, where the core-strahlo model is used to describe the electron population. To obtain the dispersion relation we linearize the Vlasov-Maxwell system of equations. This is a well-known method (Stix 1962;Krall & Trivelpiece 1973) and assumes that the small amplitude perturbations of the relevant quantities are plane waves, allowing the Vlasov-Maxwell system to be rewritten in the form: where E k is the complex amplitude of the electric field perturbation and D(ω, k, f j ) is the dispersion tensor (which is associated with the dielectric tensor of the plasma). This tensor depends on the wave vector k, the complex wave frequency ω = ω r + iγ and the background distribution functions of the j species composing the plasma f j . The dispersion relation, ω = ω(k) is determined by the condition |D(ω, k, f j )| = 0, so that Eq. (8) has non-trivial solutions for E k .
As a first approximation to the problem, we focus our attention on wave modes that propagate parallel to the background magnetic field B 0 = B 0ẑ , so that k = kẑ. We make this restriction because the mathematical analysis is greatly simplified compared to the oblique case and because previous works have shown that the field-aligned WHFI has larger growth rates ) than the oblique case. We use the core-strahlo distribution given by Eq.
(1) as the background distribution f e for the electrons. As already mentioned, to perform the integrals involved, we assume that the electron skewness is small, i.e., δ 3 s 1. This approximation allows us to obtain expression for the dispersion tensor elements D i = D i (ω, k, pp) up to second order in δ s for the parallel propagating modes. These elements depend on the wavenumber k, the wave frequency ω and the macroscopic parameters of the initial distribution functions (number density, temperature, etc) here denoted as a whole by pp. Furthermore, as the distribution is a superposition of Maxwellian and skew-Kappa, the elements of D i depend on the Fried and Conte plasma dispersion function Z(ξ) (Fried & Conte 1961), and also on the modified dispersion function Z κ (ξ) (Hellberg & Mace 2002;Viñas et al. 2015;Moya et al. 2021). Full expressions of each element of the dispersion tensor can be found in Appendix B. Here we describe the results obtained in the analysis of the excitation of the parallel propagating whistler mode associated with the heat-flux instability.
To obtain the linear properties of the WHFI, we solve the complex dispersion relation using our own developed dispersion solver. We use the core-strahlo distribution (1) to model the eVDF, and the proton population is described by an isotropic Maxwellian such that quasi-neutrality and zero-current conditions are both fulfilled. Throughout this analysis we denote proton and electrons parameters with sub-indexes p and e respectively. We fix the proton distribution so that β p = 0.1, and T ⊥p /T p = 1.0, where β j = 8πn j k B T j /B 2 0 corresponds to the parallel plasma beta of the population j. Also, k B is the Boltzmann constant, and n j is the species number density. For the electrons, we fix the anisotropy of both components (core and strahlo) equal to one, i.e, T ⊥s /T s = T ⊥c /T c = 1.0 so that there is no free energy associated with the anisotropy of the eVDF. We also fix the density of the strahlo population to 10%. With this selection of parameters for protons and electrons, the only relevant non-thermal features of the electron distribution throughout this work, are enhanced tails represented by the κ s parameter, and the skewness represented by δ s . We perform the stability analysis of the WHFI for different values of δ s , κ s , β s and T s /T c , to study how the dispersion relation depends on these parameters. Further, as in the solar wind at 1 AU from the Sun, in all our calculations we have fixed the ratio between the electron plasma frequency (ω pe ) and gyrofrequency (Ω e ) to ω pe /|Ω e | = 200.  can see that the real frequency remains essentially the same when we modify κ s and the growth rates slightly decrease as κ s increases. The case κ s = 3.0 corresponds to the most unstable case. As shown in Figure 3(c) and 3(d), κ s does not control the symmetry of the distribution. Therefore, it is also reasonable the wave stability to slightly depend on κ s . Regarding the dependence on the skewness parameter, Figure 4(b) shows the dispersion relation of the whistler mode for different values of δ s , fixing κ s = 3.0. We can see that, in the wavenumber range shown, the real part of the frequency does not change considerably when we increase δ s . The imaginary part, however, depends more strongly on this parameter, and the wave becomes more unstable as δ s increases. Also, the wavenumber range in which the mode is unstable widens, and the k value corresponding to the maximum growth rate also increases with increasing δ s . Further, note that this relation is non-linear. For example, when δ s = 0.1 the maximum growth rate reaches a value γ max ∼ 5 × 10 −4 |Ω e |, whereas for δ s = 0.2 the maximum growth rate is γ max ∼ 2 × 10 −3 |Ω e |. This behavior is expected because δ s represents a measurement of the system's free energy associated with the distribution skewness. As δ s increases, the more skew the distribution becomes, as we saw in Figure 3(a) and 3(b). Therefore, it is also expected the relation between the maximum growth rate of the WHFI and the skewness parameter to be non-linear.  Figure 4(c) shows similar plots but considering two drifting isotropic Maxwellians given by Eq. (2), so that electrons follow a current-free core-halo model as in Gary et al. (1994).
Under this model, the normalized electron q e /q 0 heat-flux along the magnetic field is giving by q e /q 0 (5/3)(n c n h /n 2 e )(∆U ch /α c )(T h /T c − 1), where T h and T c are the temperature of the halo and core, respectively, and ∆U ch is the relative drift between core and halo (Gary et al. 1994).
As a comparison with the core-strahlo model, we consider core and halo with the same density, temperatures and plasma beta as core and strahlo. We also select ∆U ch such that for each δ s value shown in Figure 4(b) both models have the same heat-flux moment (we will present more details about the heat-flux of the core-strahlo model in Section 3.1). Comparing Figures 4(b) and 4(c), we can see that the dispersion relation obtained by using the skew-Kappa core-strahlo or Maxwellian core-halo models are qualitatively the same. Both models produce the same real part of the dispersion relation, and are similarly unstable to the WHFI. Nevertheless, differences do exist. In particular, for the same level of heat-flux moment the maximum growth of the Maxwellian core-halo model is shifted to larger wavenumber compared to the core-strahlo model. The slight differences between both results are due to the fact that, they are based on different mathematical functions, that have different shapes and velocity gradients (i.e. the dispersion relation depends on these gradients) in the valid domain, and perhaps also, due to the Taylor expansion to second order on the skewness parameter (not present in a Maxwellian description of the plasma). Furthermore, and to the best of our knowledge, these slight differences in the dispersion profiles between the two models have been postulated before e.g. by Abraham-Shrauner & Feldman (1977a,b) (and probably extend back to the early work of Bernstein modes; Bernstein 1958) in application to whistler and electromagnetic ion cyclotron waves in the solar wind, that indicated that wave dispersion characteristics are not only dependent on the physical moment parameters (e.g., density, temperature, drifts, heat-flux, etc.) but that they also depend on the shape of the distribution. In this case the differences lay on the lack of supra-thermal tails in the Maxwellian model or the fact that in the core-strahlo model the source of asymmetry is strongly dominated by the skewness parameter, not present in the Maxwellian core-halo approach. However, as shown by Figures 4(b) and 4(c) these are minor differences. Both models can adequately describe the WHFI in the small skewness regime. We can definitely conclude that the results of the core-strahlo model are a very good approximation to the heat-flux instability problem, since reproduces similar behavior as previously reported, reinforcing the validity of the model (see Section 2.1).
Finally, Figure 5 shows again the normalized real and imaginary frequencies (top and bottom respectively) of the whistler mode for different values of β s and T s /T c , for fixed δ s = 0.15 and κ s = 3.0. Figure 5(a) shows the dispersion relation of the whistler mode for different values of β s , and fixed T s /T c = 7.0. From the figure we can see that in this case the real part of the frequency slightly decreases as β e increases. On the other hand, as expected the imaginary frequency depends more strongly on the plasma beta. All cases shown in the plot have a range in which the growth rate is positive, so the plasma is unstable to the whistler mode under these conditions, and it is clear that for higher β s values, the maximum growth rate is also higher. Here, however, the wavenumber at which the growth rate crosses the axis from positive to negative values shifts to the left as β s increases; i.e., the range in which the growth rate becomes positive narrows. Therefore, as a general rule, we can say that for low values of k, the growth rates increase with β s and the amplitude of the waves grows faster. At the same time, as β s increases, the wave becomes stable for lower values of the wavenumber. This behavior is also expected since for higher β values the less magnetized the plasma is, meaning that the plasma is more susceptible to destabilize due to electromagnetic fluctuations (Viñas et al. 2015;Moya et al. 2021). Further, Figure 5 The growth rate also increases with the strahlo-to-core temperature ratio, a expected behaviour as a larger temperature of the only electron component providing the free energy for the instability should also result on a larger growth rate of the waves. However, the maximum value of the growth rate seems to saturate to γ max ∼ 10 −3 |Ω e | at kc ∼ 0.2 ω pe for T s /T c 5. This is an interesting result, however, as the main goal of our study is the analysis of the whistler heat-flux instability in the solar wind, from now on we will fix T s /T c = 7 and focus the analysis on the effect of the κ s and δ s parameters.

The effect of the heat-flux moment on the instability
We are particularly interested in understanding how the WHFI contributes to the electron heat-flux regulation through collisionless wave-particle interactions in the solar wind. Thus, we expand our stability analysis and study how the whistler wave changes as we modify the electron field-aligned heat-flux moment of the eVDF. The connection between the parallel electron heat-flux (q e ) and the parameters describing the electron distribution function (1) can be seen in Eq. (9). Namely, up to second order in δ s , q e is given by: where Ψ 6 and Ψ 7 are functions that depend only on κ s and µ j = T ⊥j /T j is the anisotropy of population j (see in Appendix A for details). Further, to express the heat-flux as a dimensionless quantity it is customary to normalize q e to the free-streaming or saturation heat-flux q 0 = (3/2) n e k B T c α (see e.g Gary et al. 1994) Taking this into consideration, we can write the normalized heat-flux as follows: In this expression it is clear that the normalized electron heat-flux increases linearly with δ s . In this case, (when all other parameters are fixed) as we increase δ s , q e /q 0 also increase and the plasma becomes more unstable to the whistler mode, just like we saw in figure 4(a). In other words, if the heat-flux increases linearly with δ s , a larger value of q e /q 0 corresponds to a more skew distribution, and therefore indicates a larger level of free energy to excite the WHFI. In contrast, functions Ψ 5 (κ s ) and Ψ 6 (κ s ) indicate that the heat-flux decreases as κ s increases. In this case, as we increase q e /q 0 , the stability of the plasma to the whistler mode remains essentially the same, just like figure 4b suggest, i.e. when the increase in heat-flux values is a consequence of changes in κ s , such increment is due to the modification of the high energy tails, which are enhanced when κ s decreases but do not changes the symmetry of the eVDF.
To further analyze the behavior of the heat-flux parameter, in Figure 6 we plot the normalized growth rates γ/|Ω e | of the whistler mode as a function of the normalized wave number for different values of the initial normalized electron heat-flux q e /q 0 . We calculate this parameter using different combinations of κ s and δ s , and fixing β s = 1.0 and T s /T c = 7. Again, we considered isotropic electron populations (µ c = µ s = 1.0) and a 10% density for the strahlo (n s /n e = 0.1). Blue lines correspond to q e /q 0 = 0.02, red lines to q e /q 0 = 0.03 and green lines to q e /q 0 = 0.04. The line styles differentiate the combinations of parameters used. We see that for the same value of q e /q 0 , the stability of whistler mode changes depending on the chosen parameters. In this plot, the combinations with higher δ s (dashed lines) are always more unstable to this mode: the growth rates are positive in a wider wavenumber range, and the maximum growth rate is higher as well. Another thing that should be noticed is that for higher κ s values higher δ s values are needed to achieve the same heat-flux value.
Accordingly, for a fixed value of κ s , the heat-flux parameter is a direct measurement of the distribution function's skewness and hence, of the plasma stability (see solid lines in Figure 6). In this case, we can safely say that the higher the initial heat-flux value, the more unstable the whistler mode will be. In contrast, for a fixed δ s , we cannot make the same straightforward association between the heat-flux and the plasma stability since changing the heat-flux value will not necessarily affect the stability of the whistler mode. Things get more interesting when we allow the variation of both parameters in calculating the initial heat-flux, because the same value of q e /q 0 can be achieved using different combinations of κ s and δ s . Since only the latter parameter significantly impacts the distribution's skewness, different combinations will have different stability for the whistler mode. In other words, systems whose distributions have different levels of asymmetry and, therefore, different stability to the whistler mode, can have the same heat-flux value. Hence, the heat-flux parameter can no longer be a direct measurement of this non-thermal feature (the asymmetry of the eVDF), which gives the plasma the free energy to radiate electromagnetic waves. In consequence, it is not definitive to assure that higher heat-flux values represent more unstable states.

WHFI INSTABILITY THRESHOLDS
As mentioned, one of the goals of this research is the understanding about under which plasma conditions the WHFI develops, and also to contrast the theoretical predictions with observational data. To do so, in this section, we systematize our analysis on the excitation of the WHFI, and calculate the instability thresholds for this wave mode in the q e /q 0 vs β s space. We calculate the normalized maximum growth rate of the WHFI, γ max /|Ω e |, as a function of the normalized electron heat-flux q e /q 0 and the electron beta parameter β s . Regarding the macroscopic plasma parameters, for protons we consider the same parameters as mentioned in previous section. For the electrons, as in the previous section we consider isotropic populations (T ⊥s /T s = T ⊥c /T c = 1.0) and a density of 10% for the strahlo population (n s /n e = 0.1). We also use κ s = 3.0 and T s /T c = 7.0. With this choice of parameters, q e /q 0 is a direct measure of the distribution skewness, as can be seen in Eq.
(10). Figure 7 shows a contour plot of the maximum growth rate of the WHFI for 0.1 ≤ β s ≤ 10, and 1.8 × 10 −3 ≤ q e /q 0 ≤ 0.45, which following Eq. (10) corresponds roughly to 0.001 ≤ δ e ≤ 0.25. From the figure we can see that the general behavior for γ max is to increase to the right and upwards in the plot. i.e., as expected, the waves become more unstable as β s and q e /q 0 increases. Therefore, we recover the behaviors seen in figures 4(a) and 5(a).  Table 1. Best fit parameters for the γ max /|Ω e | = 10 −3 threshold of the whistler heat-flux instability. The curve fitting for this thresholds was performed using the function shown in Eq. (11) for different κ s , and fixing µ c = µ s = 1.0 and T s /T c = 7, n s /n e = 0.1 we fit these stability thresholds using a generalized Lorentzian function given by: and we adjust Eq. (11) to the contour γ max /|Ω e | = 10 −3.0 . In table 1 we show the parameters A, B, β 0 and α of the best fit for every value of κ s . We show these parameters to allow easier comparison between these instability thresholds and solar wind observations. Figure 8 shows all of these fits and how the threshold change for different κ s values. Blue, red and green lines correspond to κ s = 3.0, κ s = 4.0, and κ s = 5.0, respectively. From the figure we can see that for a fixed q e the plasma, predominantly, becomes more unstable as κ e increases. This behavior is consistent with the results shown in Section 3.1, because as κ s increases, higher δ s values are needed in order to achieve the same heat-flux. In other words, as κ s increases, more skew distributions are needed to achieve a given q e /q 0 , so more free energy is available in the system to excite waves, which translate into higher growth rates, as we saw in figure 4(a). Therefore, all these results seem to strengthen our previous conclusion that it is not possible to make a direct relation between the heat-flux moment and the stability of the plasma to the WHFI.

SUMMARY AND CONCLUSIONS
Using linear kinetic theory, we have performed a stability analysis of the parallel propagating whistler mode associated with the heat-flux instability in a non-collisional, and magnetized plasma, considering parameters typical of the solar wind. Using the Vlasov-Maxwell system, we calculated the Notwithstanding the above, since the skew-Kappa distribution always has a singularity for any finite value of δ s , the use of the core-strahlo model is restricted to electron distributions exhibiting small skewness (i.e., δ 3 s 1), such that the singularity is far away from the main core of the distribution in units of the thermal speed. As shown in Section 2.1, in those cases all relevant features of the VDF can be represented through a Taylor series approximation around δ s = 0, and the singularity can be evaded. Using this model, in the case of small skewness we have studied the effects and sensitivity of different plasma parameters over the excitation of the whistler-heat-flux instability, such as κ s , β s , T s /T c and the skewness parameter δ s . We focused our attention on the effect that the macroscopic parameter q e has in the mode's stability since this is the parameter that has been customarily used as a measurement of the eVDF skewness, which is the non-thermal feature in the system that provides the energy for the excitation of the WHFI. We also characterized the β dependent marginal stability thresholds as a function of the parallel electron heat-flux parameter, and present threshold conditions for the instability that can be modeled to compare with observational data.
Our results showed that when δ s > 0 the plasma is unstable to the parallel propagating WHFI, and that the growth rates increase with increasing δ s , which is the parameter that controls the skewness of the eVDF. In addition, we also showed that κ s (the parameter controlling the extent of the highenergy power-law tails of the distribution) has a weak effect on the stability of this mode: as we increase its value, the mode becomes slightly more stable. Furthermore, we presented the analytical expression for the normalized electron heat-flux. For the isotropic case (µ e = 1) we showed that it is not possible to definitely predict how the growth rates will modify as we increase the electron heat-flux. This behavior results from the fact that a given q e /q 0 value can be achieved by multiple (κ s , δ s ) combinations. Therefore the stability of the whistler mode depends greatly on how q e /q 0 is calculated in terms of δ s and κ s . Finally, to allow comparison with observations, we presented the best-fit parameters, where a generalized Lorentzian has been used for the curve fitting of these stability thresholds. Considering that in this model only δ s controls the distribution skewness, and that high δ s values (rather than high q e /q 0 value) are consistent with more unstable states, our results suggest that studies regarding the excitation of the WHFI should be mostly focused on the distribution skewness (a purely kinetic property of the VDF) instead of the heat-flux moment, which is a fluid quantity of the plasma. Under this context, we expect these results to provide a new framework to study the role of the WHFI on the depletion of the field-aligned electron heat-flux below the values predicted by the collisional transport model and clearly observed by Bale et al. (2013). Moreover, we also expect our results to motivate the search for new methods to measure or estimate the skewness or asymmetry of the eVDF from observations, rather than the field-aligned heat-flux moment.
That being said, is worth mentioning that a skew distribution can also be unstable to other microinstabilities. In particular, as recently shown by López et al. (2020a,b), electron distribution functions composed by a core and a beam can be unstable to several instabilities. However, the electrostatic instability is the fastest growing mode only when the relative drift between core and beam is larger than the thermal speed. Furthermore, even when the electrostatic mode is faster than the electromagnetic instability, its saturation level is also faster but lower than the electromagnetic. Therefore the electromagnetic mode dominates in the non-linear regime. In the case of our study, by construction the relative drift between core and strahlo is smaller than the thermal speed. Thus, the electrostatic instability may be present, but the instability triggered by the skewness (the heat-flux instability) should dominate. Due to this reason, and also because of the particular interest of the community on the possible role of the WHFI on the regulation of solar wind heat-flux through wave-particle interactions, we have decided to focus only on the whistler heat-flux instability (WHFI). Nevertheless, we acknowledge the possible existence of more instabilities due to the different non-thermal properties of the electron distribution (temperature anisotropy, asymmetry, relative drifts, power-law tails, etc).
We believe it is important to address the coexistence and interplay of these instabilities, and we expect to perform such analysis in a subsequent study.
As previously mentioned, to our knowledge the skew-Kappa function has never been used in a space plasma context. In the original derivation of this type of distribution in the context of turbulent flows, Beck (2000) showed that the inverse of the asymmetry term 1/δ s , is proportional to the square root of the Reynolds number (Re 1/2 ). Also, a quick calculation shows that Re is inversely proportional to the Knudsen number Kn (Re ∝ Kn −1 ). As Kn is directly related with the heat-flux transport in a collisional plasma (Spitzer Jr & Härm 1953;Bale et al. 2013), this relation between δ s and Kn suggests that δ s can potentially be related to parameters relevant to the turbulence phenomenon and plasma collisionality, making a possible connection between the kinetic properties of the plasma (the skewness) and a fluid description of the media (the Knudsen number). Therefore, the core-strahlo model based on the skew-Kappa distribution can be more than an ad-hoc representation of solar wind electrons. Even though the detailed explanation of such relation is beyond the scope of these work, we strongly believe that theoretical studies that can link these parameters in plasma systems should be investigated further. For example, as mentioned in Bale et al. (2013), the transition between the collisional Spitzer-Härm heat-flux transport and the collisionless regimes occurs at Kn ∼ 0.3. Therefore, using Eq. 10, in terms of δ s our model predicts that such transition should occur at δ s ∼ 0.2 for κ s = 3, n s /n e = 0.1, T s /T c = 7.0, and µ c = µ s = 1.
In summary, through the introduction of the core-strahlo model we have showed that plasma states with the same initial heat-flux could have different stability to the WHFI. In other words, systems with high q e /q 0 values can be stable enough so that the WHFI cannot modify the electron heat-flux values effectively through wave-particle interactions. Thus, the heat-flux by itself does not seem to be the best indicator of the WHFI. The precise source of the a field-aligned heat-flux instability is the skewness of the distribution, and in the case of the skew-Kappa, such non-thermal feature is clearly represented by the skewness parameter δ s . The use of the skew-Kappa eVDF allows the distribution skewness to be controlled through just one parameter, and considerably reduces the space of free parameters to analyze. Furthermore, the core-strahlo model represents an alternative approach to the usual functions used to phenomenologically model the eVDF in terms of a superposition of three electron sub populations: core, halo and strahl. This new model to describe the eVDF greatly simplify the study of WHFI, and the role this instability plays in the electron thermal energy transport in the solar wind. Finally, the use of this distribution has a potential theoretical justification that should be explored as δ s could be related to the turbulence phenomenon in plasma systems, and also be relevant for the understanding of the relation between the strahl and the formation of the halo during the expansion of the solar wind from the outer solar corona to the heliosphere. Under this context, since the general case with arbitrary skewness remains to be solved, our results suggest that future research about the collisionless heat-flux transport in the solar wind should centered on the distribution skewness (encapsulated here by the δ s parameter), rather than just the heat-flux parameter. We expect this theoretical analysis inspired by observations to be relevant and provide valuable insights in the solar wind electron heat-flux regulation debate. Hopefully, the results shown here will be validated with experimental data in lights of the new Parker Solar Probe and Solar Orbiter missions.

ACKNOWLEDGMENTS
We thank the support of ANID, Chile through the Doctoral National Scholarship N (3). In order to perform the integrals in velocity space involved in these expressions, we assume a small skewness (i.e., δ 3 s 1) and compute these macroscopic parameters using a Taylor expansion of Eq (3) up to second order in δ s . As usual, in the following equations the subscripts and ⊥ are with respect to the background magnetic field.

B. DISPERSION TENSOR
In a non-colissional and uniform plasma, immersed in a background magnetic field B 0 = B 0ẑ , for parallel propagating waves k = kẑ, the dispersion tensor D can be written in the following form: The restriction |D(ω, k, f j )| determines the relation between the wave frequency ω = ω r + iγ and the wavenumber k for the parallel propagating modes. In Eq. (B15) the dispersion tensor elements are written in terms of the susceptibilities χ i (f j ), where the sums are carried on over all species j Λ 1 = 1 2 κ s + 1 κ s − 3