The Effects of Nonequilibrium Velocity Distributions on Alfvén Ion-cyclotron Waves in the Solar Wind

In this work, we investigate how the complex structure found in solar wind proton velocity distribution functions (VDFs), rather than the commonly assumed two-component bi-Maxwellian structure, affects the onset and evolution of parallel-propagating microinstabilities. We use the Arbitrary Linear Plasma Solver, a numerical dispersion solver, to find the real frequencies and growth/damping rates of the Alfvén modes calculated for proton VDFs extracted from Wind spacecraft observations of the solar wind. We compare this wave behavior to that obtained by applying the same procedure to core-and-beam bi-Maxwellian fits of the Wind proton VDFs. We find several significant differences in the plasma waves obtained for the extracted data and bi-Maxwellian fits, including a strong dependence of the growth/damping rate on the shape of the VDF. By applying the quasilinear diffusion operator to these VDFs, we pinpoint resonantly interacting regions in velocity space where differences in VDF structure significantly affect the wave growth and damping rates. This demonstration of the sensitive dependence of Alfvén mode behavior on VDF structure may explain why the Alfvén ion-cyclotron instability thresholds predicted by linear theory for bi-Maxwellian models of solar wind proton background VDFs do not entirely constrain spacecraft observations of solar wind proton VDFs, such as those made by the Wind spacecraft.


INTRODUCTION
Space and astrophysical plasmas are commonly weakly collisional.With infrequent collisions, the velocity distribution functions (VDFs) of these plasmas can significantly differ from the entropically-favored Maxwellian shape (Vasyliunas 1968;Gosling et al. 1981;Lui & Krimigis 1981;Marsch et al. 1982;Armstrong et al. 1983;Lui & Krimigis 1983;Williams et al. 1988).These departures from equilibrium can act as sources of free energy that drive instabilities.Instabilities play an important role in both the small-and large-scale plasma behavior, e.g. by pitch-angle scattering of particles acting as an effective viscosity (Kunz et al. 2011;Riquelme et al. 2018;Arzamasskiy et al. 2023) or by serving as channels to transport energy in space; see the review by Matteini et al. (2012).In the weakly collisional plasma that comprises the solar wind a variety of non-equilibrium features develop, including temperature anisotropy (Kasper et al. 2002), temperature disequi-librium between species (Neugebauer 1976), and beams (Alterman et al. 2018;Pilipp et al. 1987).
To better describe the non-equilibrium features in solar wind VDFs, it is common to fit observations to simplified model distributions such as bi-Maxwellian (Marsch 2006) or kappa (Summers et al. 1994) distributions.Bi-Maxwellian models have been extensively used to determine instability thresholds as a function of proton temperature anisotropy (T ⊥,p /T ∥,p where perpendicular/parallel are defined with respect to the mean magnetic field B) and the parallel proton plasma beta (β ∥,p = 8πn p kT ∥,p /B 2 where n j is the density of a plasma component j).Assuming a bi-Maxwellian shape of the background VDFs and a homogeneous background, different kinds of wave instabilities are commonly expected to be driven beyond these instability thresholds.The thresholds can account for either a single free energy source (e.g.Hellinger et al. 2006;Bale et al. 2009) or multiple sources (e.g.Matteini et al. 2013;Chen et al. 2016;Klein et al. 2018); reviews While the color scale is linear, gray contours show the shape of fp across 10 logarithmically-spaced levels between 0.01 and 0.9 of these thresholds are provided by Gary (1993); Verscharen et al. (2019).
A well-known open question in heliophysics, and the focus of this paper, concerns the Alfvén ion-cyclotron (AIC) instability, which occurs when the proton temperature is larger perpendicular to the mean magnetic field direction than parallel, T ⊥,p > T ∥,p (Kennel & Wong 1967;Davidson & Ogden 1975).The predictions of linear theory for the AIC instability threshold do not fully constrain the observations of solar wind proton VDFs by the Wind spacecraft when calculated for bi-Maxwellian proton background VDFs (Hellinger et al. 2006;Bale et al. 2009).Several explanations have been proposed for the failure of this approach to describe the observations, including inefficient energy extraction (Shoji et al. 2009) or the impact of other plasma populations such as α-particles (Maruca et al. 2012).Isenberg et al. (2013) suggest that the assumption of bi-Maxwellian proton distributions limits the applicability of the thresholds predicted by bi-Maxwellian linear theory.Our investigation focuses on understanding how the structure of realistic solar wind VDFs affects the onset and evolu-tion of the AIC instability in the solar wind compared to the bi-Maxwellian approximation.
The structure of a plasma VDF has physical consequences in terms of growth or damping of plasma waves, which are important carriers of energy in collisionless plasmas such as the solar wind.We are therefore interested in understanding how the properties of waves differ between those produced from core-and-beam bi-Maxwellian models of spacecraft observations of solar wind proton VDFs and the spacecraft data itself.
Linear plasma solvers (such as WHAMP, NHDS, or PLUME) designed to study the plasma waves associated with a given background VDF typically use a bi-Maxwellian model and solve the hot-plasma dispersion relation associated with this simplified description of the distribution (Roennmark 1982;Quataert 1998;Klein & Howes 2015;Verscharen & Chandran 2018).In our investigation of how deviations from an idealized bi-Maxwellian distribution affect the linear properties of the AIC instability, we use ALPS (the Arbitrary Linear Plasma Solver), a code that solves the hot-plasma dispersion relation by direct numerical integration of arbitrary background dis-tributions for the plasma; a detailed overview of ALPS is provided by Verscharen et al. (2018).
In section 2 we describe the VDFs used in this work and outline the numerical tools employed to calculate the linear plasma response.We present the plasma behavior of the Alfvén modes associated with these VDFs in section 3, along with a more detailed analysis of the connection between VDF structure and the stability of the Alfvén modes through calculation of the quasilinear diffusion operator.We summarize our results in section 4.

VDF Construction
The VDFs used in this paper are based on two representative distributions observed by the Solar Wind Experiment (SWE) aboard the Wind spacecraft (Ogilvie et al. 1995).The Wind spacecraft rotates with a period of ≈ 3 seconds and charge flux distributions, as a function of energy and incidence angle, are acquired over ≈ 93 seconds.This rotation allows the distribution function to be resolved over the full field of view of 4π steradians.The distributions used in this work are based on observations made over one acquisition period on 2009 December 19 from 02:49:41 UT and 2008 August 14 from 04:31:43.69UT for Interval 1 and Interval 2 respectively.The average magnetic field during each period is measured by the Magnetic Field Investigation (MFI) (Lepping et al. 1995) and the average magnitude of the magnetic field over this 93 second interval is 5.044 nT for Interval 1 and 6.065 nT for Interval 2. Based on a nonlinear bi-Maxwellian fit to the SWE proton measurements (Kasper et al. 2002), we find β ∥,p = 0.153, T ⊥,p /T ∥,p = 2.41 (Interval 1) and β ∥,p = 0.148, T ⊥,p /T ∥,p = 2.80 (Interval 2).Both intervals are therefore above the γ/Ω p = 10 −2 marginal stability threshold (where γ is the growth/damping rate, Ω p = q p B/m p is the proton cyclotron frequency, q p is the proton charge, and m p is the proton mass) for the AIC instability if the distribution function were modeled as a bi-Maxwellian, but below the same threshold for the mirror instability (Verscharen et al. 2016).We choose intervals in this region of β ∥ , T ⊥,p /T ∥,p space to allow us to characterize the role of non-equilibrium VDF structure on the stability of solar wind plasma that should be unstable to the AIC instability under the bi-Maxwellian assumption.These intervals exhibit a prominent proton beam population; such intervals have historically been treated with core-and-beam bi-Maxwellian fits.The α-particle density is below the detection threshold for the Wind non-linear analysis for Interval 1 and around n α /n p ∼ 0.007 for Interval 2. These low α-particle den-sities do not contribute significantly to the linear plasma response, and thus we will not consider the α-particles in this analysis.
The Wind SWE Faraday cup has a baseline resolution of 31 bins logarithmically distributed across voltages from 150 to 8000 V with energy resolution of ∆E/E ≈ 0.065 (Ogilvie et al. 1995) We further generalize the bi-Maxwellian VDF model typically extracted from these measurements, provided by Ogilvie et al. (2021), to better describe the Wind SWE Faraday Cup observations.The Faraday Cup response function is not uniquely or analytically invertible, so a numerical approach is required in order to recover details of the VDF that are not well-described by the bi-Maxwellian model.We choose to characterize the VDFs more completely by following a two-step approach: first, we capture the dominant features with a double bi-Maxwellian regression; second, we discretize that model and apply random gyrotropic perturbations to improve the match with the data.The resulting discrete gyrotropic models are suitable for analysis with ALPS.
We first choose a double bi-Maxwellian VDF model to capture both the core and secondary beam proton components that are evident in the data.The Faraday Cup data is natively expressed as ion charge fluxes measured as a function of voltage and spacecraft spin angle (Kasper et al. 2021).Following the protocol used in several previous works, (Alterman et al. 2018;Chen et al. 2016;Jiansen et al. 2018;Wicks et al. 2016;Gary et al. 2016), we fit the measured instrument response to the modeled instrument response for two bi-Maxwellian distributions.We superpose two instances of the Faraday Cup response to a single bi-Maxwellian, where the analytic form for the latter is described by Kasper et al. (2006) and the sensitive area as a function of angle is derived for each datum by interpolating on the effective area table provided in the metadata of Kasper et al. (2021).We then perform regressions to the data, with gyrotropic constraints applied to the bi-Maxwellian symmetry axes and relative drifts.The resultant bestfit parameters obtained using this method are hereafter referred to as "model 1," as they are labeled in Table 1a.
To explore non-Maxwellian variations on the fits, we first discretize each model 1 fit to a Cartesian (v ⊥ , v ∥ ) grid with 2 km/s x 2 km/s resolution and then recompute the Faraday Cup response numerically.Calculating the Faraday Cup response requires one to integrate the differential charge flux upon the sensor, multiplied by the effective sensitive area, over an appropriate domain.In this case, we take the differential charge flux associated with any point in phase space from the gridded model and retrieve the effective area from the table provided in Kasper et al. (2021).The integration domain for each datum is the conical frustum in phase space defined by the cup's full field of view and the bounds of the particular energy window.The integrations are performed using the IDL (Interactive Data Language) "int 3d" utility with twenty-point Gaussian quadrature.
For both spectra, we calculate the non-reduced chisquare statistic by comparing the numerical model to the data and then employ a Monte Carlo strategy.At each iteration, (1) a pseudorandom perturbation is superposed at each point in the discrete VDF model, (2) the Faraday Cup response to the perturbation is calculated, (3) the chi-square is recalculated for the perturbed model, and (4) the perturbation is accepted or rejected on the basis of whether the chi-square was reduced.Gaussian perturbations are generated in an ad hoc manner: the generator selects randomized magnitudes between -0.1 and 0.1 of the local VDF value and half-widths between 1 and 10 grid cells.To select phase space locations, the generator uses the VDF as a probability distribution function, with an ad hoc flattening and artificial enhancement out to six thermal widths.This is found to strongly favor low-energy solutions while still providing a clear and reliable signal at the tails.We carry out the Monte Carlo procedure until the chi-square statistic converges to 1 part in 10,000.The resulting modified spectrum models are hereafter referred to as the "inverted-interpolated data" or "I-I data" as they are labeled in Fig. 1.
We also use the Levenberg-Marquardt method to construct core-and-beam bi-Maxwellian fits to the I-I data.These best-fit models are hereafter referred to as "model 2," as they are labeled in Table 1a.As this work aims to understand the impacts of VDF structure on microinstabilities, we include both of these independentlyconstructed models to determine whether small differences between bi-Maxwellian models also yield different linear plasma response.The fit parameters for both model 1 and model 2 are shown in Table 1a.
Figure 1 shows the I-I data and two core-and-beam bi-Maxwellian models for both intervals, where fp is a re-scaling of f p to a maximum value of 1.The contours highlight how a two-component bi-Maxwellian fit is unable to capture the complex structure of the I-I data, at both small and large v ∥ .

Calculation of the Dispersion Relation
We employ ALPS to calculate the linear kinetic plasma behavior of the bi-Maxwellian model and I-I data VDFs.For full code details, see Verscharen et al. (2018).ALPS seeks non-trivial solutions to the wave equation by nu-merical calculation of the plasma susceptibility of each component species.Unlike dispersion solvers that simplify this process by modeling the distributions as bi-Maxwellians (e.g.PLUME, NHDS) or kappa distributions (e.g.DSHARK (Astfalk et al. 2015)), ALPS numerically integrates the susceptibility directly from a given gyrotropic background distribution f j (v ⊥ , v ∥ ).By numerically integrating the phase space density to solve the plasma wave equation rather than using the properties of special functions such as Maxwellians to find closedform solutions of the dispersion relation, ALPS allows us to analyze the I-I Wind data directly, without making the simplifying assumption that the VDF structure follows a closed analytical form.
To numerically solve the hot plasma dispersion relation, laid out by Stix (1992), ALPS uses the distributions f j for each species j defined on a discrete grid in perpendicular and parallel momentum space, with minimum and maximum values P min,j and P max,j and resolution defined by the number of steps N ⊥ or N ∥ between P min,j and P max,j in each direction.The calculation of the susceptibility χ j is given in equation 2.9 of Verscharen et al. (2018).A user defined parameter J max determines the order of resonances n max to include in the calculation of χ j such that the maximum amplitude of the Bessel function J nmax+1 < J max .Modes that resonantly interact with portions of the distribution outside the resolved numeric grid are unaffected by the particular shape of f j and have been shown to match fluid solutions.
For damped modes, the wave equation requires an analytic continuation, for which we apply the hybridanalytic continuation scheme described in section 3.2 of Verscharen et al. (2018).The Landau contour integration in this scheme is computed by decomposing the integral into a compound function, where numerical integration is used when possible and the residue component of the function is computed using an algebraic function fit to f j that can be selected by the user.Additional numeric methods are needed to integrate the poles; these are described in section 3.1 of Verscharen et al. (2018).For these integrations, t lim and M I are user defined parameters determining the regions over which these numerical procedures are implemented.
For our calculation of the dispersion relation using ALPS, we use the following parameters.For all cases, J max = 10 −50 , M I = 5, M p = 100, T lim = 0.01, and P min,⊥j = 0m p v A .For Interval 1, the momentum space resolution is N ⊥ = 189 and N ∥ = 206.For Interval 2, the momentum space resolution is N ⊥ = 178 and N ∥ = 203.The proton and electron momentum space ranges are presented in Table 1b.The electron VDFs used for the dispersion calculation in ALPS were con- The ω r solutions for the first (blue) and second (green) bi-Maxwellian models are qualitatively similar to the I-I data (red) solution, but the γ behavior differs significantly.For instance, the bi-Maxwellian models are unstable for the backward Alfvén solution for Interval 1 while the I-I data yields a stable solution.The forward Alfvén solution in Interval 2 is unstable in model 2 while model 1 and I-I data are stable.The ALPS scan is performed for k ∥ d p between 0.01 and 10, and we choose the limited range presented here to constrain the scan to the region where the resonant velocity is within the limits of the input distribution and γ/ω r < 1/e.structed using the observed electron temperature from Wind SWE (13 eV for Interval 1 and 12.5 eV for Interval 2), with relative drift and density chosen to preserve current-free and quasi-neutral conditions.

Dispersion Relation
Figure 2 presents the real and imaginary components of the linear mode frequencies, ω r and γ respectively, for the forward and backward Alfvén modes as a function of the normalized parallel wavenumber, k ∥ d p (where d p = v A /|Ω p | is the proton inertial length).We calculate the dispersion relation for the parallel-propagating Alfvén modes with normalized perpendicular wavenumber k ⊥ d p = 0.001 and k ∥ d p between 0.01 and 10.The k ∥ d p ranges shown in Fig. 2 start at k ∥ d p = 0.25 as this is the lower limit where the resonant velocity v res = (ω r − nΩ p ) /k ∥ (where n is an integer) exceeds the v ∥ bounds of the input distributions.For the upper limit of k ∥ , we require that γ/ω r < 1/e for consistency with our later calculations of the quasilinear diffusion operator.While the calculated real frequency ω r agrees for the data and bi-Maxwellian models for all modes (with slight differences emerging as the waves approach small scales k ∥ d p ≫ 1), there is significant disagreement in γ.The backward Alfvén mode for Interval 1 shows unstable behavior -a positive γ -for both bi-Maxwellian models, but the solution for the I-I data is stable for all k ∥ d p .The forward Alfvén mode for Interval 2 shows stable behavior in the data and model 1, but model 2 develops an instability.
Even in the cases for which all models predict an instability, the extent of the unstable wavenumber range in k ∥ d p differs.The forward Alfvén mode in Interval 1 ex- G is plotted against the right-hand y-axis and given by the cyan line with error bounds as described in the text, with the associated growth (solid green line) or damping (dashed) rates plotted against the left-hand y-axis.G is dimensionless but proportional to (n j /n 0 )/(v/v A ) 7 .Where G is positive, we find a corresponding instability (positive γ).For the forward Alfvén mode, gray contours are 8 logarithmically-spaced levels between |5e − 5| and |.01| for both positive (solid) and negative (dashed) values.For the backward Alfvén mode, contours are 8 logarithmically-spaced levels between -0.25 and -1e-2 and 8 logarithmically-spaced levels between 1e-5 and 0.05.hibits only a small region of instability from k ∥ d p = 0.33 to 0.41 for the I-I data, compared to both bi-Maxwellian model instabilities which range from k ∥ d p = 0.25 to 0.41 (model 1) and k ∥ d p = 0.25 to 0.57 (model 2).The backward Alfvén mode in Interval 2 has instabilities that span similar k ∥ d p ranges for the bi-Maxwellian models, with k ∥ d p = 0.25 to 1.04 (model 1) and k ∥ d p = 0.33 to 1.09 (model 2), but the instability in the I-I data spans only k ∥ d p = 0.25 to 0.79, terminating at a smaller k ∥ d p than either of the models.
To tie the differences in unstable behavior back to the structure of the input VDFs, we consider the resonant velocities associated with the wavenumber at which the instabilities terminate.These are indicated in Figure 1 by the vertical colored lines, with any v res < 0 corresponding to the forward Alfvén modes and v res > 0 corresponding to the backward Alfvén modes.In Interval 1, the forward Alfvén mode in the I-I data interacts with a region further from the core than either of the models, by which point the model VDF contours are quite smooth.The backward Alfvén modes for both models resonant with a region of smooth f p structure for v ∥ > 0 and large temperature anisotropy.The I-I data exhibit a more uneven structure and steeper drop-off of f p in this same region, and no instability in the backward Alfvén wave is found.For the Interval 2 forward Alfvén mode, we find an instability in model 2, but not model 1, with the width of the core population being the primary difference in VDF structure between the models.While the ALPS results for the I-I data exhibit an instability at -174 km/s, this is very close to the bounds of the input VDF and is likely an artifact of the ALPS scan moving into the resolved VDF range, so we consider this to be an overall stable mode.In the backward Alfvén mode, all three resonant velocities shown in Fig. 1 are close to the core population, with the I-I data instability terminating at a larger v res than either model, where the I-I data VDF has a much more varied structure than either model.

Diffusion Operator
To better understand what drives the differences in the linear plasma response between the I-I data and bi-Maxwellian models, we examine how the differences in the structure of the input distributions affect the sign and magnitude of γ.We use Eqn. 1, an analytical expression for the quasilinear growth rate γ j of a species j at a particular wavenumber.
where ω pj = 4πn j Z 2 q 2 j /m j is the plasma frequency (with Z the charge state) and the diffusion operator, Ĝ, is given by Ĝ A full derivation of Eqn. 1 from the Vlasov-Maxwell equations is provided by Kennel & Wong (1967).Eqn. 1 is valid when γ j ≪ ω r .The explicit dependence of γ j on the distribution function and the resonant velocity links VDF structure to the resultant stability or instability of a given plasma mode.The only component of Eqn. 1 that is not positive definite is Ĝf j , so this term determines whether a species contributes at the given wavenumber range through a positive or negative γ j , directly affecting the stability of a given mode based on the structure of f j (specifically on the v ⊥ and v ∥ gradients of f j seen in Eqns. 3 and 4).The electrons have a minor contribution and the intervals chosen have low α-particle density, so we consider only proton growth rate, γ p in our calculation.We neglect positive definite terms that do not impact the sign of the energy transfer, including the pre-factors before the sum over n and the wave energy density W . Evaluation of Eqn. 1 using the ALPS Alfvén solutions for the three VDFs allows us to better understand what aspects of the VDF structure lead to the strong variances in the growth rate.
The electric field interaction is captured in the ϕ n,k term: The argument of the Bessel function J is Assuming parallel propagating modes with k ⊥ = 0 limits our evaluation to n = 0, ±1 since J n̸ =0 (0) = 0.This approximation holds as long as k ⊥ ≪ k ∥ , as is the case for our calculations.We further focus on the Alfvénic solutions, which limits the calculation to n = ±1 for ω r /k ∥ ≶ 0 based on the polarization of the modes.We do not perform a full quasilinear calculation, but only use the quasilinear diffusion operator to determine the regions of phase space responsible for the growth and damping of waves.With these caveats, we evaluate Ĝf (k ∥ d p , v ⊥ /v A ) for Interval 2 using the real frequency calculated from the ALPS solutions for the two models and I-I data; these are plotted as color in Fig. 3.We define and plot this integrated value as a solid teal line in Fig. 3, with surrounding error range in cyan.Over the region in k ∥ d p where G is positive, an instability is present.The error range indicates the bounds for ±10% variation in the retrieved ALPS ω r for the Alfvén modes to account for any numerical discrepancies in the frequency resolution.The gray shaded regions are areas where γ p /ω r > 1/e.We find in both the bi-Maxwellian models shown in the top two rows of Fig. 3 that Ĝf has a smooth structure with defined lobes of positive and negative values.For the I-I data in the bottom row, patches of positive and negative regions intermingle and only upon calculation of G can the final stable or unstable behavior be discerned.The main difference between stable and unstable behavior of the I-I data for the Interval 2 VDF is the presence of a strong negative contribution at small v ⊥ (where v ⊥ /v A < 1) in the forward Alfvén mode.This originates from G ⊥ f , shown in the top-right panel of Fig. 4. Tracing this back to the distribution itself in Fig. 1, we see that over the region where Model 2 is unstable (at v ⊥ < -34 [km/s]), the I-I data VDF has more variation, particularly at low v ⊥ .This structure is not captured in either of the models, which show almost symmetric structures in Ĝ⊥ and Ĝ∥ in the first and second columns of Fig. 4. Similarly, in the backward Alfvén mode for Interval 2, the models display a smooth structure in both Fig. 3 and Fig. 4 whereas the I-I data structure is far more complex.Both the models and I-I data exhibit unstable behavior for this mode, indicating that the complex VDF structure and gradients does not necessarily suppress an AIC instability and that it may be more coincidence than accurate representation of the I-I data by the models that results in similar unstable behavior for all three distributions.
These features in the diffusion operator support the differences seen in the VDFs in Fig. 1, indicating that both the simple contour structure of f p and the underlying velocity gradients in the I-I data VDF are not being accurately captured by the two bi-Maxwellian models and that different regions and structures in v ⊥ and v ∥ space contribute to the resultant stable or unstable behavior of the plasma normal modes.Even in the case of the two bi-Maxwellian models, very slight differences in the density, temperature anisotropy, and drift speed lead to different behavior, as seen in Fig. 4. Slight differences in the relative values of the perpendicular and parallel contributions lead to an overall net negative γ p value, and thus a stable forward Alfvén mode for Model 1, but an unstable mode for Model 2.

DISCUSSION AND CONCLUSIONS
We utilize ALPS to solve the hot-plasma dispersion relation for inverted-interpolated solar wind proton VDFs based on two intervals observed by the Wind spacecraft and compared these to the dispersion relation for coreand-beam bi-Maxwellian models of these VDFs.These intervals are chosen with β ∥,p and T ⊥,p > T ∥,p such that linear theory predicts that the AIC instability will develop under a bi-Maxwellian assumption.We study each interval using three different VDFs.The first VDF, model 1, is the typical core-and-beam bi-Maxwellian fit to the measured instrument response.We next implement a novel method to extract a more realistic VDF from the Wind Faraday cup measurements; we use a Monte Carlo strategy to apply Gaussian perturbations to a discrete version of model 1.This procedure optimizes the chi-square statistic to produce a VDF that more closely matches the known Faraday cup response.This VDF, which we call the I-I data, better represents the underlying distribution.Our model 2 VDF is constructed from a core-and-beam bi-Maxwellian fit to the I-I data VDF rather than to the instrument response.
The more involved procedure used to construct the I-I data yields a distribution with many non-Maxwellian features that reproduces the measured instrument response more accurately than either of the model fits.
For the two selected intervals, a comparison between the model and I-I data VDFs shows that the bi-Maxwellian models result in fundamentally different plasma behaviors than those calculated for the I-I data they are based on, both by failing to find an instability when one is calculated from the I-I data or by finding significantly different k ∥ ranges over which the instability operates.We show that unstable behavior in Alfvén ion-cyclotron waves is sensitively dependent on small VDF structures, particularly in the region where v ⊥ /v A < 2. As bi-Maxwellian models are commonly used to study solar wind proton VDFs, these findings indicate that such simplifications may not accurately represent instabilities in the solar wind, and thus incompletely capture wave-particle energy transfer associated with these VDFs.
We use the quasilinear diffusion operator to demonstrate how the structural differences between the I-I data and model VDFs produce the dissimilar growth rates seen in the ALPS solutions.Even for modes where both the model and I-I data VDFs yield the same behavior in γ (such as the Interval 2 backward Alfvén mode where an instability is found in all three cases), calculation of the diffusion operator, Ĝf j , for the I-I data reveals a far more complex structure.These underlying structural differences between the I-I data and model VDFs are apparent even when considering the separate Ĝ⊥ and Ĝ∥ terms, revealing the failure of these models to truly capture the f p dependence on either v ⊥ or v ∥ that is found in the I-I data.
This work demonstrates the importance of including distribution structure when studying the development of microinstabilities.The AIC instability threshold, calculated from bi-Maxwellian models of solar wind protons, is frequently exceeded in solar wind observations.Such observations are unexpected as unstable AICs should reduce the temperature anisotropy towards equilibrium.When we consider a more complex VDF structure (the I-I data in this paper), entirely different stable or unstable AIC behavior is calculated compared to bi-Maxwellian models, potentially explaining this discrepancy.

Figure 1 .
Figure 1.Bi-Maxwellian core-and-beam models (left and center columns) and I-I data (right column) for Intervals 1 (top row) & 2 (bottom).Vertical lines correspond to the resonant velocity at which an instability terminates for the I-I data (red), first bi-Maxwellian model (blue), and second bi-Maxwellian model (green).Dashed and solid lines indicate resonant velocities associated with the forward and backward Alfvén modes respectively.The velocity range is trimmed in this figure to highlight the structure of the VDF, with the actual velocity ranges used in the ALPS calculation spanning v ∥ ∈ [−207, 205] km/s and and v ⊥ ∈ [0, 378] km/s for Interval 1 and v ∥ ∈ [−204, 202] km/s and v ⊥ ∈ [0, 356] km/s for Interval 2. While the color scale is linear, gray contours show the shape of fp across 10 logarithmically-spaced levels between 0.01 and 0.9

Figure 2 .
Figure2.ALPS dispersion solutions for real and imaginary frequencies ω r and γ.The forward (left column) and backward (right) Alfvén solutions are presented for Interval 1 (top panels) and 2 (bottom).The ω r solutions for the first (blue) and second (green) bi-Maxwellian models are qualitatively similar to the I-I data (red) solution, but the γ behavior differs significantly.For instance, the bi-Maxwellian models are unstable for the backward Alfvén solution for Interval 1 while the I-I data yields a stable solution.The forward Alfvén solution in Interval 2 is unstable in model 2 while model 1 and I-I data are stable.The ALPS scan is performed for k ∥ d p between 0.01 and 10, and we choose the limited range presented here to constrain the scan to the region where the resonant velocity is within the limits of the input distribution and γ/ω r < 1/e.

Figure 3 .
Figure 3. Diffusion operator Ĝf (Eqn.2) as a function of k ∥ and v ⊥ indicates regions where energy is transferred from the distribution to the fields (red) or vice versa (blue).Colorbars at bottom row apply to entire column.Ĝf is dimensionless but proportional to (n j /n 0 )/(v/v A ) 4 where n 0 ≡ d 3 f j .Ĝ calculations are shown for Interval 2 forward Alfvén mode (left) and backward Alfvén mode (right) for the first bi-Maxwellian model (top), second bi-Maxwellian model (middle), and I-I data (bottom).G is plotted against the right-hand y-axis and given by the cyan line with error bounds as described in the text, with the associated growth (solid green line) or damping (dashed) rates plotted against the left-hand y-axis.G is dimensionless but proportional to (n j /n 0 )/(v/v A ) 7 .Where G is positive, we find a corresponding instability (positive γ).For the forward Alfvén mode, gray contours are 8 logarithmically-spaced levels between |5e − 5| and |.01| for both positive (solid) and negative (dashed) values.For the backward Alfvén mode, contours are 8 logarithmically-spaced levels between -0.25 and -1e-2 and 8 logarithmically-spaced levels between 1e-5 and 0.05.

Figure 4 .
Figure 4. Perpendicular (top row) and parallel (right) diffusion operators, Eqn. 3 and Eqn. 4, for the bi-Maxwellian models 1 and 2 (left and middle columns) and I-I data (right) for the Interval 2 forward Alfvén solution.The color scheme follows that presented in Figure 3.The models have qualitatively similar structures, but exhibit different Ĝ⊥ f and Ĝ∥ f .This can be compared against the complexity seen in the I-I data, where the same broad structures are not observed.For model 1 and the I-I data, gray contours are 8 logarithmically-spaced levels between |5e − 4| and |0.03| for both positive (solid) and negative (dashed) values.For model 2, gray contours are 8 logarithmically-spaced levels between |5e − 3| and |0.3| for positive and negative values.

Table 1 .
Plasma and ALPS parameters used in this work.