<figure><img src="image_1.png" /></figure>

# Boundary of the Distribution of Solar Wind Proton Beta versus Temperature Anisotropy

P. H. Yoon¹, M. Lazar²,³, C. Salem⁴, J. Seough⁵, M. M. Martinović⁶,⁷, K. G. Klein⁶, and R. A. López⁸

¹Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742-2431, USA; yoonp@umd.edu

²Centre for Mathematical Plasma Astrophysics, Celestijnenlaan 200B, B-3001 Leuven, Belgium

³Institut für Theoretische Physik, Lehrstuhl IV: Weltraum- und Astrophysik, Ruhr-Universität Bochum, D-44780 Bochum, Germany

⁴Space Sciences Laboratory, University of California, Berkeley, CA 94720, USA

⁵Korea Astronomy and Space Science Institute, Daejeon 34055, Republic of Korea

⁶Lunar & Planetary Laboratory, University of Arizona, Tucson, AZ 85721-0092, USA

⁷LESIA, Observatoire de Paris, Meudon, France

⁸Research Center in the intersection of Plasma Physics, Matter, and Complexity (P²inc), Comisión Chilena de Energía Nuclear, Casilla 188-D, Santiago, Chile

Received 2024 January 9; revised 2024 May 2; accepted 2024 May 2; published 2024 July 3

## Abstract

The frequency distribution of solar wind protons, measured in the vicinity of Earth's orbit, is customarily plotted in $(\beta_\|, T_\perp/T_\|)$ phase space. Here, $T_\perp/T_\|$ is the ratio of perpendicular and parallel temperatures, and $\beta_\| = 8\pi n T_\|/B^2$ is the ratio of parallel thermal energy to background magnetic field energy, the so-called "parallel beta," with $\perp$ and $\|\|$ denoting directions with respect to the ambient magnetic field. Such a frequency distribution, plotted as a two-dimensional histogram, forms a peculiar rhombic shape defined with an outer boundary in the said phase space. Past studies reveal that the threshold conditions for temperature anisotropy-driven plasma instability partially account for the boundary on the high-$\beta_\|$ side. The low-$\beta_\|$ side remains largely unexplained despite some efforts. Work by Vafin et al. recently showed that certain contours of collisional relaxation frequency, $\nu_{pp}$, when parameterized by $T_\perp/T_\|$ and $\beta_\|$, could match the overall shape of the left-hand boundary, thus suggesting that the collisional relaxation process might be closely related to the formation of the left-hand boundary. The present paper extends the analysis by Vafin et al. and carries out the dynamical computation of the collisional relaxation process for an ensemble of initial proton states with varying degrees of anisotropic temperatures. The final states of the relaxed protons are shown to closely match the observed boundary to the left of the $(\beta_\|, T_\perp/T_\|)$ phase space. When coupled with a similar set of calculations for the ensemble in the collective instability regime, it is found that the combined collisional/collective effects provide the baseline explanation for the observation.

*Unified Astronomy Thesaurus concepts: Solar wind (1534); Interplanetary physics (827); Space plasmas (1544); Heliosphere (711); Plasma physics (2089)*

*Supporting material: animations*

## 1. Introduction

The origin and dynamics of solar wind is a subject of ongoing research, especially in view of contemporary inner heliospheric missions, specifically Parker Solar Probe (Fox et al. 2016; Bale et al. 2023; Raouafi et al. 2023) and Solar Orbiter (Müller et al. 2020; Opie et al. 2022). The charged particles in the solar wind (protons, electrons, or alpha particles) are customarily organized as frequency distributions, or a two-dimensional (2D) histogram of physical states, in the phase space defined by parallel plasma beta, $\beta_{\|,s}$, and the ratio of perpendicular and parallel temperatures (defined with respect to the ambient magnetic field **B**), $T_{\perp,s}/T_{\|,s}$, where $T_{\perp,s}$ and $T_{\|,s}$ are perpendicular and parallel temperatures for each species (labeled *s*). The parallel plasma beta for species *s* is defined by $\beta_{\|,s} = 8\pi n_s T_{\|,s}/B^2$, where $n_s$ and *B* stand for the ambient density for species *s* and magnetic field intensity, respectively. A substantial body of literature exists that discusses the properties of charged particles in the solar wind as they relate to their distributions in $(\beta_{\|,s}, T_{\perp,s}/T_{\|,s})$ phase space (Kasper et al. 2002, 2003; Marsch et al. 2006; Štverák et al. 2008; Bale et al. 2009; Maruca et al. 2011, 2012; He et al. 2013; Matteini et al.

2013; Hellinger & Trávníček 2014; Chen et al. 2016; Klein et al. 2018, 2019, 2021; Perrone et al. 2019; Huang et al. 2020). For the rest of the present paper, we will only concern ourselves with the proton temperatures, densities, and betas. Henceforth, we omit the subscript *s*.

An interesting feature associated with the data distributions in $(\beta_\|, T_\perp/T_\|)$ phase space is that the outer boundaries in the high-beta regime, that is, the boundaries to the right, are partially defined by various threshold conditions for temperature anisotropy instability (Gary et al. 1997, 2001; Marsch et al. 2004; Hellinger et al. 2006; Matteini et al. 2007; Yoon 2017; Zhao et al. 2019; Sun et al. 2019, 2020; Verscharen et al. 2019; Martinović 2021; Xiang et al. 2021) with the peak broadly located in quasi-isotropic conditions, $T_\perp \sim T_\|$. The outer boundaries in the low-beta regime, that is, the boundaries to the left, remain largely unexplained, although several suggestions have been put forward. These include that based on the propagation effect along the spiraling magnetic field of the solar wind (Matteini et al. 2012), Alfvén-cyclotron resonant heating (Chandran et al. 2011; Cranmer 2014; Ozak et al. 2015), the effects of a relative drift in the proton–proton system (Ibscher & Schlickeiser 2013, 2014; Ibscher et al. 2013; Chen et al. 2016), and the result of driven turbulence (Servidio et al. 2014). Note that the nature of the secondary drifting ion populations (including both protons and alphas) has recently been overviewed (Alterman et al. 2018).

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
It is interesting to note that the threshold conditions for the temperature (or more precisely, the pressure) anisotropy instability of the type discussed in the context of the solar wind are also pertinent to astrophysical plasma problems, such as in the intracluster medium dynamics, astrophysical dynamo and turbulence generation processes, and for astrophysical charged-particle acceleration processes (Santos-Lima et al. 2016; St-Onge & Kunz 2018; Bott et al. 2021; Kim et al. 2021; Arzamasskiy et al. 2023). The notion of various plasma instabilities regulating the temperature anisotropy actually began with observations made in the Earth's magnetosheath in the first place (Anderson et al. 1991; Gary et al. 1993a, 1993b; Phan et al. 1994; Tan et al. 1998). Later, this notion became quite popular in the context of solar wind research, as already overviewed.

While the solar wind near Earth's orbit is largely collision-free, the influence of binary collisions can nevertheless be nonnegligible, particularly on the relaxation of temperature anisotropies (Hellinger & Trávníček 2014; Kasper et al. 2017). Hellinger & Trávníček (2009, 2010) computed the temperature relaxation transport coefficient based upon the Landau collision integral, assuming that the underlying particle velocity distribution function is given by a bi-Maxwellian model. Vafin et al. (2019) made use of the proton collisional relaxation frequency, $\nu_{pp}$, derived in the above-referenced papers and demonstrated that contours of $\nu_{pp}$, when parameterized in $T_{\perp}/T_{\parallel}$ and $\beta_{\parallel}$, could match the overall shape of the left-hand boundary, thus suggesting that the collisional relaxation process might be closely related to the formation of the left-hand boundary. However, they did not actually compute the dynamical collisional relaxation process itself.

The purpose of the present paper is to build upon the work by Vafin et al. (2019) and to actually compute the dynamical progression of how the collisional relaxation of bi-Maxwellian protons proceeds. To complement the analysis, we also carry out the collisionless relaxation of bi-Maxwellian protons in the high-beta regime by excitation of instabilities and quasi-linear relaxation. The organization of the present paper is as follows. In Section 2, we overview the theory of Vafin et al. (2019), where we propose a simpler and approximate description of the same problem. Section 3 discusses the numerical example of temperature relaxation processes by collisions. Section 4 expands upon the temperature relaxation problem by including the collisionless processes, that is, quasi-linear wave excitation and saturation processes. Section 5 puts both the collisional and collective relaxation calculations together and plots the combined results against the background of the 2D data distribution, thus highlighting the essential findings of the present study.

## 2. Proton-Proton Collisional Relaxation Frequency in ($\beta_{\parallel}$, $T_{\perp}/T_{\parallel}$) Phase Space

Taking the collisional frequency for proton-proton interactions as derived by Hellinger & Trávníček (2009), which is based upon the Landau collision integral under the assumption of a drifting bi-Maxwellian plasma distribution function, Vafin et al. (2019) proceeded to write down the equation for proton temperature relaxation. Note that Jubeh & Barghouthi (2017) have also carried out a similar calculation to what Hellinger & Trávníček (2009) did. In the absence of relative drift in the bi-Maxwellian plasma, the collision frequency derived by Hellinger & Trávníček (2009) and Jubeh & Barghouthi (2017)

is actually identical to the well-known formula available in the Naval Research Laboratory (NRL) Plasma Formulary (Richardson 2019), which in turn is taken from Kogan (1961). Note that ignoring the secondary protons and relative drift is an idealization. For actual solar wind, the secondary protons and alphas that are often seen to be drifting with respect to the background protons are quite common, as recently overviewed by Alterman et al. (2018). Herewith, we restrict ourselves to the relatively simple case of a single proton species characterized by the temperature anisotropy only. We therefore assume such a simplifying configuration and choose to express the collisional relaxation equation following the convention in the NRL Plasma Formulary,

$$ 
\begin{aligned}
\frac{dT_{\perp}}{dt} &= \nu_{pp}(T_{\parallel} - T_{\perp}), \\
\frac{dT_{\parallel}}{dt} &= 2\nu_{pp}(T_{\perp} - T_{\parallel}), \\
\nu_{pp} &= \frac{2\pi^{1/2}ne^4 \ln \Lambda}{m_p^{1/2} T_{\perp}^{3/2}} F,
\end{aligned}
\quad (1) $$

where $\nu_{pp}$ is the collisional frequency, $m_p$ is the proton mass, $e$ is the unit electric charge, $\ln\Lambda = \ln(4\pi n\lambda_D^3)$ is the Coulomb logarithm, $\lambda_D = T_e^{1/2}/(4\pi n e^2)^{1/2}$ is the Debye length, and the quantity $F$ is given by

$$ F = \frac{(A + 1)^{3/2}}{A^2} \left( (A + 3) \frac{\arctan \sqrt{A}}{\sqrt{A}} - 3 \right), \quad (2) $$

if $A > 0$,

$$ F = \frac{(A + 1)^{3/2}}{A^2} \left( (A + 3) \frac{\arctanh \sqrt{-A}}{\sqrt{-A}} - 3 \right), \quad (3) $$

if $A < 0$, and

$$ F = \frac{4}{15} \quad (4) $$

if $A = 0$, where

$$ A = \frac{T_{\perp}}{T_{\parallel}} - 1. \quad (5) $$

While this formula is exact and well known, the actual derivation is not easily found in the literature. For the sake of completeness, we present a brief derivation in Appendix A, the purpose of which is not only to recount the actual derivation, but also to point out that Equations (1)–(5) can be obtained from a more general Balescu–Lenard collision operator (Balescu 1960; Lenard 1960), as compared to the reduced Landau collision integral (Landau 1936), which Hellinger & Trávníček (2009) and Jubeh & Barghouthi (2017) adopted as the starting point of their analyses.

We find that the exact definition of $\nu_{pp}$ is not necessary, but for the present purpose, it is sufficient to employ an approximate form,

$$ \nu_{pp} = \frac{8\pi^{1/2}ne^4 \ln \Lambda}{15m_p^{1/2} T_{\perp}^{3/2}}, \quad (6) $$

which follows from the approximation $F \approx 4/15$, valid when $|A| = |T_{\perp}/T_{\parallel} - 1|$ is sufficiently small—Equation (4). We have
arrived at this conclusion by comparing the numerical solutions based upon the full expression for F against the above approximation, and having found that the solutions are virtually identical for a large range of input conditions. We will return to this issue briefly when we discuss the dynamical solutions later. For the moment, let us rewrite Equation (1) in terms of dimensionless quantities:

$$ 
\begin{aligned}
\beta_{\perp, \parallel} &= \frac{8\pi n T_{\perp, \parallel}}{B^2}, \\
\tau &= \Omega_p t, \\
g &= \frac{\sqrt{2}}{\pi^{3/2}} \frac{1}{n(c/\omega_{pi})^3} \frac{c^4}{v_A^4},
\end{aligned} 
\quad (7) $$

and by making use of the simplified collision frequency (6). The resulting normalized temperature (or beta) relaxation equation is

$$ 
\begin{aligned}
\frac{d\beta_{\perp}}{d\tau} &= \mu(\beta_{\parallel} - \beta_{\perp}), \\
\frac{d\beta_{\parallel}}{d\tau} &= 2\mu(\beta_{\perp} - \beta_{\parallel}), \\
\mu &= \frac{\nu_{pp}}{\Omega_p} = \frac{g \ln \Lambda}{15 \beta_{\perp}^{3/2}},
\end{aligned} 
\quad (8) $$

where

$$ 
\begin{aligned}
\Omega_p &= \frac{eB}{m_p c}, \\
\omega_{pp}^2 &= \frac{4\pi n e^2}{m_p}, \\
v_A^2 &= \frac{B^2}{4\pi n m_p}, \\
\Lambda &= 4\pi n \lambda_D^3 = \frac{T_e^3}{(4\pi n)^{1/2} e^3},
\end{aligned} 
\quad (9) $$

denote the proton cyclotron frequency, the square of the proton plasma oscillation frequency, the square of the Alfvén speed, and the argument of the Coulomb logarithm, in that order. In Equation (7), $\beta_{\perp, \parallel}$ is the perpendicular/parallel plasma beta, $\tau = \Omega_p t$ is the dimensionless time variable, and $g$ is a quantity that measures the collisionality of the plasma ($g=0$ being the collisionless limit). We should note that $g$ defined in Equation (7) is somewhat different from a similarly defined quantity in the case of collisional processes in magnetized plasmas. For example, the same quantity in Yoon et al. (2019), where the magnetized plasma formalism is adopted, is defined by $g = (n\rho_p^3)^{-1}(c/v_A)^4$, where $\rho_p = [2T_{\parallel}/(m\Omega_p^2)]^{1/2}$ denotes the proton thermal gyroradius. By contrast, the quantity $g$ in Equation (7), which is derived from the unmagnetized collision operator—see Appendix A, is defined in terms of the proton inertial length, $d_p = c/\omega_{pp}$.

Following Vafin et al. (2019), we consider a quantity of interest,

$$ R(\beta_{\parallel}, A) = \frac{1}{\beta_{\perp}^{3/2}} \left| \frac{\beta_{\perp}}{\beta_{\parallel}} - 1 \right| = \frac{1}{\beta_{\parallel}^{3/2} (A + 1)^{3/2}} \quad (10) $$

<figure><img src="image_3.png" /><figcaption>Figure 1. The contours of the function $R(\beta_{\parallel}, A) = \beta_{\parallel}^{-3/2}|A|(A+1)^{-3/2}$, which is closely related to the proton-proton collision frequency $\nu_{pp}$, plotted against the distribution of the background solar wind data in $(\beta_{\parallel}, T_{\perp}/T_{\parallel})$ phase space. Numerical values indicated by yellow boxes denote the contour levels for $R$.</figcaption></figure>

Clearly, $R(\beta_{\parallel}, A)$ is proportional to $\nu_{pp}$ or its normalized version $\mu$. Vafin et al. (2019) plotted a quantity related to $R(\beta_{\parallel}, A)$ and noticed that certain contours coincide with the distribution of solar wind proton data in 2D phase space $(\beta_{\parallel}, T_{\perp}/T_{\parallel}) = (\beta_{\parallel}, A+1)$. To illustrate this point, we show in Figure 1 the contour plot of the function $R$ against the backdrop of the distribution of 2D solar wind proton data, where contour levels are indicated by numerical values enclosed in yellow boxes. The background proton data distribution is derived in Faraday cup observations in the Solar Wind Experiment on the Wind spacecraft (Ogilvie et al. 1995), to which a nonlinear least-squares fitting algorithm is applied to determine moments of the solar wind proton velocity, namely density, bulk speed, temperature, and temperature anisotropy (Maruca & Kasper 2013). As Figure 1 demonstrates, certain contours of the function $R$ closely match the left-hand outer boundaries of the data distribution, especially for the $T_{\perp} < T_{\parallel}$ region (the lower left boundary). Vafin et al. (2019) state that the “Coulomb collisions are an interesting candidate for explaining” the left-hand boundaries, and also note that other processes could contribute as well. They conclude that their “results may provide the lower limits” or the baseline explanation for the observation. While they did not actually carry out the dynamical process of collisional relaxation, they did provide the justification for a particular value of the contour, namely, that which corresponds to $\nu_{pp} = 2.8 \times 10^{-6} \text{ s}^{-1}$, which they chose to fit the boundary. Their value, of course, makes sense from a back-of-the-envelope calculation: taking 1 au as the travel distance, and 400 km s$^{-1}$ as the solar wind speed gives a transit time of $3.73 \times 10^5 \text{ s}$, whose reciprocal is $\sim 2.6 \times 10^{-6} \text{ s}^{-1}$.

### 3. Proton Temperature Relaxation by Collisions: Dynamical Process

To provide a more solid basis for determining the importance of the collisional relaxation process as it relates to the formation of the low-$\beta_{\parallel}$ boundaries, we further analyze the problem by
actually solving for the dynamical equations over a hypothetical ensemble of “initial” solar wind proton states. Our implicit assumption is that this initial ensemble of protons is related to its source close to the corona and propagates to the near-Earth vicinity, but it should be emphasized that we are not actually attempting to model the realistic solar wind conditions close to the solar source (Huang et al. 2020; Mostafavi et al. 2024). Rather, our purpose is just to cover the entire phase space uniformly by choosing a hypothetical ensemble of protons for the initial conditions. As we normalize the physical quantities according to the scheme outlined in Equation (8), it is necessary to estimate the strength of the magnetic field, the density, and the electron temperature. Upon making use of the well-known NRL Plasma Formulary (Richardson 2019),

$$ 
\begin{aligned}
\Omega_p &= 9.58 \times 10^3 B_{[\text{Gauss}]} \text{ rad s}^{-1}, \\
\frac{c}{\omega_{\text{pp}}} &= \frac{2.2727 \times 10^7}{n_{[\text{pcc}]}^{1/2}} \text{ cm}, \\
\frac{c}{\nu_A} &= \frac{0.1376 n_{[\text{pcc}]}^{1/2}}{B_{[\text{Gauss}]}}, \\
\Lambda &= \frac{5.16 \times 10^9 T_{[\text{eV}]}^{3/2}}{n_{[\text{pcc}]}^{1/2}},
\end{aligned} 
\quad (11) $$

we find that the crucial quantities g and lnΛ necessary for the numerical computation of ensemble dynamics are given by

$$ 
\begin{aligned}
g &= 7.7565 \times 10^{-27} \frac{n_{[\text{pcc}]}^{5/2}}{B_{[\text{Gauss}]}^4}, \\
\ln\Lambda &= \ln \frac{5.16 \times 10^9 T_{[\text{eV}]}^{3/2}}{n_{[\text{pcc}]}^{1/2}}.
\end{aligned} 
\quad (12) $$

If we use the typical 1 au value for the plasma density, $n \sim 8 \text{ cm}^{-3}$, the magnetic field intensity, $B \sim 5 \text{ nT} = 5 \times 10^{-5} \text{ G}$, and the electron temperature, $T_e \sim 140,000 \text{ K} \sim 12 \text{ eV}$ (Wilson et al. 2018; Klein & Veck 2019), then we have $\Omega_p \sim 0.4790 \text{ rad s}^{-1}$ and $g \ln\Lambda/15 \sim 3.7520 \times 10^{-7}$. If we consider that 1 au equals $1.496 \times 10^8 \text{ km}$ and also that the typical solar wind speed is $\sim 300-600 \text{ km s}^{-1}$, then the time it takes for a typical solar parcel to travel from its solar source to 1 au is $\sim 2.5 \times 10^5-5 \times 10^5 \text{ s}$. This implies that the maximum value for dimensionless time is $\Omega_{pt} \sim 10^4-10^5$. In this first-order approximation, we work with a mean value of the magnetic field at 1 au, and consider only collisions between particles of the same species, in particular, proton-proton collisions, as in Vafin et al. (2019).

On the basis of such estimations, we now carry out the numerical evolution of an ensemble of initial data points presumed to originate from a solar source to 1 au. We thus solve the dynamical Equation (8) with input parameters

$$ g' = \frac{g \ln \Lambda}{15} = 10^{-7}, \quad \tau_{\text{max}} = 5 \times 10^4. \quad (13) $$

Note that these estimates are approximately the same as Vafin et al.’s value. In Figure 2, we show the initial ensemble (the hypothetical solar wind conditions near the solar source) by showing each point as an open circle. The ensemble points, which uniformly cover the low-beta region of the phase space, is plotted in the background of the quantity R, which is related

<figure><img src="image_4.png" /><figcaption>Figure 2. The initial ensemble points representing a hypothetical set of solar wind proton conditions near the solar source, indicated by open circles, are plotted against the backdrop of a function defined by $R = g'\beta_{\parallel}^{-3/2}|A|(A+1)^{-3/2}$, which is closely related to the proton-proton collision frequency $\nu_{pp}$.</figcaption></figure>

to that used by Vafin et al. (2019), that is, Equation (10), except that we have multiplied R by $g' = g \ln\Lambda/(15\beta_{\parallel}^{3/2})$. Thus, the gray scale indicates the magnitude of the quantity (here, we remind readers that we work with the normalized collisional frequency $\mu = \nu_{pp}/\Omega_p$)

$$ \frac{g \ln \Lambda}{15 \beta_{\parallel}^{3/2}} \frac{|A|}{(A + 1)^{3/2}} \quad (14) $$

in logarithmic value. This is an appropriate place to return to the issue of the approximation $F \sim 4/15$. As readers may appreciate, some ensemble points are characterized by relatively large $T_{\perp}/T_{\parallel}$ values (e.g., 0.1 or 10), for which one might justifiably question the validity of the approximation. However, we find that the net result, which we will discuss subsequently, is indistinguishable from the approximated solution. The reason is as follows: while for such extreme initial ensembles, the relatively early dynamical processes indeed proceed at different rates, over a longer time, these ensemble points all eventually settle down close to marginally collisional states (where the collision frequency $\mu$ is extremely low). To put it another way, regardless of the initial states, the final states of the ensemble points are indistinguishable, whether they are computed from exact F or its approximated form. From this, we have arrived at the conclusion that the approximation $F \approx 4/15$ is valid for a large range of initial ensemble states, but it is particularly applicable for marginally collisional states, $\mu \ll 1$.

Figure 3 shows an intermediate result when the ensemble calculations have proceeded to normalized time corresponding to $\tau = \Omega_{pt} = 10^3$. The dynamical paths of initial ensemble points in $(\beta_{\parallel}, T_{\perp}/T_{\parallel}) = (\beta_{\parallel}, A+1)$ phase space are shown with magenta-colored lines, and the state of ensemble points at $\Omega_{pt} = 10^3$ is indicated by black dots. As one can see, certain ensemble points evolve rather swiftly and settle down to the final states as indicated by black dots, while some other ensemble points move rather sluggishly,
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The dynamic paths of initial ensemble points in $(\beta_\|, T_\perp/T_\|) = (\beta_\|, A)$ phase space are shown with magenta-colored lines, and the state of ensemble points at $\Omega_p t = 10^3$ is indicated by black dots.</figcaption></figure>

or almost not at all, such that for such states, the “final” ensemble points, indicated by black dots, almost coincide with the initial ensemble points shown with open circles. As such, the open circles are hidden behind the black dots. It turns out that the ensemble points that do not evolve happen to correspond to the state of solar wind protons that are well within the rhombic 2D data distribution, while those points that move rather promptly happen to be those that are well outside the same distribution. Thus, when we solve the dynamic Equation (8) up to the maximum normalized time period, namely $\tau_{max} = \Omega_p t_{max} = 5 \times 10^4$, i.e., the characteristic transit time for the protons to move from the solar source to 1 au, then the result is that the initial ensemble points would have all moved close to the left-hand (or low-$\beta_\|$) boundary of the 2D data distribution, while those points that are already within the data distribution hardly move.

To illustrate this finding, we show in Figure 4 the initial ensemble points as open dots, dynamic paths in $(\beta_\|, T_\perp/T_\|) = (\beta_\|, A+1)$ phase space as blue curves, and the final states (at the end of computation corresponding to $\Omega_p t = 5 \times 10^4$) as black dots, against the backdrop of the 2D histogram. (Note that we have also computed the same ensemble equations with exact $F$ and found a virtually identical result, which is shown in Appendix B for the sake of completeness and to prove our earlier point that the present approximate approach is sufficiently accurate.) It is quite obvious that the “initial” ensemble points that are sufficiently far away from the rhombic data distribution all move quickly toward and settle around the outer perimeter to the left within the time frame, $\Omega_p t = 5 \times 10^4$, which is estimated to be comparable to the transit time between the solar source and 1 au. By contrast, those ensemble points that are already within the data distribution barely move. In short, Figure 4 provides an additional, and more solid, justification for the baseline explanation, which Vafin et al. (2019) put forward, for the left-hand boundary of the proton data in terms of the collisional relaxation. Recall that Vafin et al. (2019) chose one particular contour for the collision frequency, $\nu_{pp} = 2.8 \times 10^{-6} \text{ s}^{-1}$, which is derived based upon physical arguments, to fit the

<figure><img src="image_6.png" /><figcaption>**Figure 4.** The initial ensemble points, shown with open dots, dynamic paths in $(\beta_\|, T_\perp/T_\|)$ phase space, and the final states at $\Omega_p t = 5 \times 10^4$, shown with black dots, are plotted against the background of the 2D histogram. In the accompanying animation, the entire dynamic paths for the ensemble points, subject to collisional processes, can be seen.</figcaption></figure>

(An animation of this figure is available.)

outer boundary to the left, but nevertheless they did not carry out a calculation of dynamical relaxation to prove that the collisional relaxation mechanism indeed leads to an almost exact overlap between the final states of the ensemble calculations and the actual boundary of the distribution of the solar wind proton data. It is interesting to note that Hellinger & Trávníček (2014) have also investigated the collisional relaxation problem by reanalyzing the Wind spacecraft data. In Figure 2 of their paper one finds that the authors have organized the data as a function of collisionality. Their results are broadly consistent with the model in the present paper. In the video file accompanying this paper, we also present the entire course of ensemble evolution.

4. Evolution of Proton Temperature Anisotropy by
Collective Relaxation

The physics of the right-hand boundary associated with the 2D data distribution is well known—such boundaries are the result of various temperature anisotropy instabilities and saturation. In spite of this well-known fact, it is nevertheless interesting to repeat the demonstration of an instability-driven relaxation process. Such a dynamical computation not only serves the purpose of complementing the analysis solely based on the collisional relaxation—that of the previous section—but also it is informative in the sense that the numerical demonstration reveals how fast the collective relaxation is when compared to the collisional processes. As is well known, a number of temperature anisotropy instabilities are considered in the literature, but for the sake of simplicity, we confine ourselves to the excitation of collective modes propagating in the direction parallel (and antiparallel) to the ambient magnetic field. The evolution equations for temperatures subject to parallel proton-cyclotron (C) and parallel firehose (F) modes
are given by (Seough & Yoon 2012)

$$ \frac{dT_{\perp}}{dt} = -\frac{B^2}{2\pi n} \int dk \left[ \gamma_k^C W_k^C \left( 1 + \frac{\omega_k^C \omega_{pp}^2}{\Omega_p c^2 k^2} \right) \right. \\ \left. \qquad + \gamma_k^F W_k^F \left( 1 - \frac{\omega_k^F \omega_{pp}^2}{\Omega_p c^2 k^2} \right) \right], \\ \\ \frac{dT_{\parallel}}{dt} = \frac{B^2}{2\pi n} \int dk \left[ \gamma_k^C W_k^C \left( 1 + \frac{2\omega_k^C \omega_{pp}^2}{\Omega_p c^2 k^2} \right) \right. \\ \left. \qquad + \gamma_k^F W_k^F \left( 1 - \frac{2\omega_k^F \omega_{pp}^2}{\Omega_p c^2 k^2} \right) \right], \quad (15) $$

where $\omega_k^{C,F}$, $\gamma_k^{C,F}$, and $W_k^{C,F}$ denote the wave frequency, growth rate, and the wave intensities for the cyclotron/firehose mode. The wave intensity $W_k^{C,F}$ is determined at each time according to the wave kinetic equation, $\partial W_k^{C,F}/\partial t = 2\gamma_k^{C,F} W_k^{C,F}$, and the wave frequency $\omega_k^{C,F}$ and the growth rate $\gamma_k^{C,F}$ are computed at each time step by solving the instantaneous dispersion relation,

$$ \begin{aligned} 0 &= \frac{c^2 k^2}{\omega_{pp}^2} \pm \frac{\omega}{\Omega_p} - A - \eta Z(\zeta), \\ \eta &= \frac{(A + 1)\omega \mp A\Omega_p}{k\alpha_{\parallel}}, \quad \zeta = \frac{\omega \mp \Omega_p}{k\alpha_{\parallel}}, \end{aligned} \quad (16) $$

where $\omega = \omega_k^{C,F} + i\gamma_k^{C,F}$ is the complex frequency, $A = T_{\perp}/T_{\parallel} - 1$ is the temperature anisotropy factor (which we have already encountered), $\alpha_{\parallel} = \beta_{\parallel}^{1/2} v_A$ corresponds to the parallel proton thermal speed, and the upper/lower sign applies for C/F modes. This set of equations can be derived from the self-consistent quasi-linear kinetic equation with the assumption of a bi-Maxwellian proton distribution function (Seough & Yoon 2012; Yoon 2017). Such a scheme, which can be termed the velocity moment quasi-linear theory, has been used for a number of applications (Seough & Yoon 2012; Yoon & Seough 2012; Seough et al. 2013, 2023; Yoon & Sarfraz 2017; Yoon et al. 2019, 2023), and has compared quite well with the particle-in-cell or hybrid simulations (Seough et al. 2014, 2015; Sarfraz et al. 2022; López et al. 2023).

Let us rewrite Equation (15) by making use of the same dimensionless quantities that we have already introduced, namely, $\tau = \Omega_p t$, $\beta_{\perp,\parallel}$, $A = T_{\perp}/T_{\parallel} - 1 = \beta_{\perp}/\beta_{\parallel} - 1$, as well as in terms of additional normalized variables,

$$ z_{C,F} = \frac{\omega_k^{C,F}}{\Omega_p}, \quad \kappa = \frac{ck}{\omega_{pp}} = \frac{kv_A}{\Omega_p}. \quad (17) $$

The result is

$$ \begin{aligned} \frac{d\beta_{\perp}}{d\tau} &= -4 \sum_{+, -} \int d\kappa z_i^{\mp} W_{\kappa}^{\mp} \left( 1 \pm \frac{z_r^{\mp}}{\kappa^2} \right), \\ \frac{d\beta_{\parallel}}{d\tau} &= 4 \sum_{+, -} \int d\kappa z_i^{\mp} W_{\kappa}^{\mp} \left( 1 \pm \frac{2z_r^{\mp}}{\kappa^2} \right), \\ \frac{\partial W_{\kappa}^{\mp}}{\partial \tau} &= 2z_i^{\mp} W_{\kappa}^{\mp}, \\ 0 &= \kappa^2 + z_{\mp} \mp A - \eta_{\mp} Z(\zeta_{\mp}), \\ \eta_{\mp} &= \frac{(A+1)z_{\mp} \mp A}{\kappa \sqrt{\beta_{\parallel}}}, \quad \zeta_{\mp} = \frac{z_{\mp} \mp 1}{\kappa \sqrt{\beta_{\parallel}}}, \end{aligned} \quad (18) $$

<figure><img src="image_7.png" /><figcaption>Figure 5. Open circles denote the initial ensemble points representing unstable conditions subject to the anisotropy-driven plasma instabilities. These dots are plotted against the background, showing the contours of the maximum growth rate for the cyclotron instability and the same for the parallel firehose instability.</figcaption></figure>

where upper/lower (or $\mp$) signs correspond to C and F modes, respectively.

Figure 5 displays the initial ensemble points subject to the anisotropy-driven plasma instabilities. These points, which uniformly cover the high-beta region of the phase space, are plotted against the background, in which the contours of maximum growth rate for the cyclotron instability, $\gamma_{\max}^C/\Omega_p$, and the same for the parallel firehose instability, $\gamma_{\max}^F/\Omega_p$, are indicated.

We have solved Equation (16) for these ensemble points up to $\Omega_p t = 20$, which is much shorter than the maximum normalized time for the collisional relaxation case. These calculations do not include the collisional relaxation, but only the quasi-linear relaxation by instability saturation. As Figure 6 exhibits, the relaxation of temperature anisotropy by instability excitation and saturation has already progressed substantially toward the marginally stable state. This shows how much more efficient the collective (that is, instability-induced) relaxation process can be when compared with the collisional relaxation process. It is also interesting to note that those ensemble points that are already located within the stable region do not move, which is as expected. In Figure 6, the dynamical paths of the initial ensemble points subject to the quasi-linear relaxation are shown with dotted magenta lines. The final states are depicted with black dots. Because of the fact that the ensemble points that are stationed within the stability boundary move very slowly or not at all, the final states depicted with black dots cover the initial open circles, but readers should be aware that there are open circles underneath these black dots. In the video accompanying this article, we also show the entire dynamics of ensemble evolution over the time period corresponding to $0 < \Omega_p t < 70$.

## 5. Summary: Combined Collisional and Collective Relaxations in ($\beta_{\parallel}$, $T_{\perp}/T_{\parallel}$) Phase Space

Let us combine both the collisional relaxation and instability relaxation that we have considered thus far, and
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The dynamical paths of the initial ensemble points subject to the quasi-linear relaxation are shown with dotted magenta lines. The final states are depicted with black dots. In the accompanying animation, the entire dynamic paths for the ensemble points that are subject to instability excitation and quasi-linear relaxation, can be seen.</figcaption></figure>

(An animation of this figure is available.)

present the findings in a single, final plot. The result is displayed in Figure 7. In the case of collisional relaxation, the computation is carried out up to $\Omega_{pt} = 5 \times 10^4$, which is identical to Figure 4. For the collective (instability-induced quasi-linear) relaxation, the quasi-linear Equation (18) for the velocity moment is solved up to $\Omega_{pt} = 70$. That is, we have combined two separate sets of ensemble calculations, one in which only collisional processes are taken into account and the other in which only collective processes are considered. One could, in principle, combine both, but for the present study, we have chosen to separate the two processes. For both cases, the dynamical paths of the initial ensemble points are shown with dotted magenta lines, and the final states are depicted with blue dots. These results are plotted in the background of the proton distribution in 2D ($\beta_\|$, $T_\perp/T_\|$) phase space.

Figure 7 thus constitutes the summary of the present discussion, and the interpretation thereof is as follows: If we compute the collisional relaxation process over a typical transit time of the solar wind parcel from its solar source to 1 au, or equivalently up to $\Omega_{pt} = 5 \times 10^4$ or so, then we find that the solar wind protons all accumulate close to the low-$\beta_\|$ boundary. Likewise, if, for some reason, the solar wind protons in the vicinity of the Earth's orbit happen to be in a state subject to instability excitation, then such a state will rather quickly relax back to the marginally stable state. Those solar wind protons that are already located within the interior of the 2D data distribution do not evolve any further and remain stationary, unless other effects are taken into account such as the inhomogeneity-driven adiabatic expansion, the effects of radial expansion, Alfvén-cyclotron heating, driven turbulence, the influence of secondary drifting species, etc.

The significance of our finding is that we have clarified the issue of what defines the outer boundary of the solar wind

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Combined results of the dynamical ensemble calculation. For the collisional relaxation, the computation is carried out up to $\Omega_{pt} = 5 \times 10^4$, which is identical to Figure 4. For the collective (instability) relaxation, the quasi-linear Equation (18) for the velocity moment is solved up to $\Omega_{pt} = 70$. For both cases, the dynamical paths of the initial ensemble points are shown with dotted magenta lines, and the final states are depicted with blue dots. These results are plotted in the background of the proton distribution in 2D ($\beta_\|$, $T_\perp/T_\|$) phase space.</figcaption></figure>

proton distribution in ($\beta_\|$, $T_\perp/T_\|$) phase space. That the high-$\beta_\|$ boundary to the right is defined by instabilities had been well known, but the low-$\beta_\|$ boundary to the left had not been clearly defined despite some earlier works cited previously. However, the work by Vafin et al. (2019), which motivated our research, has shown for the first time that the collisional relaxation might be the primary mechanism that accounts for the presence of the low-$\beta_\|$ boundary. In the present work, we have provided an additional justification for Vafin et al.'s idea in that we have computed the dynamical progression of hypothetical proton ensemble points toward the low-$\beta_\|$ boundary. We believe that the combined collisions and instabilities provide the best and most reasonable baseline explanation. Finally, we remind readers that our baseline model computation does not take into account other effects such as the radial expansion or heating by pre-existing turbulence. Nor does it contain the dynamical influence of multiple charged-particle species, including the relative drifts among them. Investigating the influence of these and other factors is a task for future research.

## Acknowledgments

This research was supported by NASA Award 80NSSC19K0827 (grant NNH18ZDA001N-HSR). This material is also based upon work funded by the Department of Energy (DE-SC0022963) through the NSF/DOE Partnership in Basic Plasma Science and Engineering. This research was also partially supported by NSF grant 2203321 to the University of Maryland. These results were also obtained (M. L.) in the framework of the projects WEAVE project—G.0025.23N/FI 706/31-1 (FWO-Vlaanderen/DFG-Germany) and SIDC Data Exploitation (ESA Prodex-12). J.S. was supported by basic research funding from the Korea
Astronomy and Space Science Institute (KASI2024185002). The authors have no conflicts to disclose.

## Data Availability

This paper deals with largely theoretical aspects except for the solar wind proton data that were used as the backdrop. The data used in this project were accessed through the SPFD CDA Web service at https://cdaweb.gsfc.nasa.gov. No new data are generated herewith.

# Appendix A
## Derivation of Proton-Proton Collision Frequency

In Hellinger & Trávníček (2009) and also in a similar analysis by Jubeh & Barghouthi (2017), the foundational equation is that of the Landau collision integral (Landau 1936), which suffers from divergences in the limits of both infinitely short and long wavelengths. The earlier work by Kogan (1961) is also based upon the same Landau collision integral. However, a more general collisional operator is that of Balescu–Lenard type (Balescu 1960; Lenard 1960), which is well behaved, at least in the infinitely long wavelength regime. Thus, the general collision integral, which we adopt, is that of the Balescu–Lenard kinetic equation,

$$ 
\begin{aligned}
\frac{\partial f_a(\mathbf{v})}{\partial t} = & -\frac{e_a^2}{m_a^2} \int d\mathbf{k} \int d\omega \left(\frac{\mathbf{k}}{k} \cdot \frac{\partial}{\partial \mathbf{v}}\right) \mathrm{Im} \frac{1}{\omega - \mathbf{k} \cdot \mathbf{v} + i0} \\
& \times \left[ \frac{m_a}{2\pi^3 k |\epsilon(\mathbf{k}, \omega)|^2} f_a(\mathbf{v}) + \langle \delta E^2 \rangle_{k, \omega} \left( \frac{\mathbf{k}}{k} \cdot \frac{\partial f_a(\mathbf{v})}{\partial \mathbf{v}} \right) \right], \\
\langle \delta E^2 \rangle_{k, \omega} = & -\frac{2}{\pi^2 k^2 |\epsilon(\mathbf{k}, \omega)|^2} \sum_a e_a^2 \int d\mathbf{v} \mathrm{Im} \frac{1}{\omega - \mathbf{k} \cdot \mathbf{v} + i0} f_a(\mathbf{v}), \\
\epsilon(\mathbf{k}, \omega) = & 1 + \sum_a \frac{4\pi e_a^2}{m_a k^2} \\
& \times \int d\mathbf{v} \frac{1}{\omega - \mathbf{k} \cdot \mathbf{v} + i0} \left( \mathbf{k} \cdot \frac{\partial}{\partial \mathbf{v}} \right) f_a(\mathbf{v}).
\end{aligned}
\quad (A1) $$

Without loss of generality, we may take the **k** vector to lie in the xz plane, **k** = k⊥**x** + k‖**z**. Let us take the velocity moments of Equation (A1) by assuming the bi-Maxwellian distribution, f<sub>a</sub>(**v**) = n<sub>a</sub>(π<sup>3/2</sup>α<sub>⊥a</sub><sup>2</sup>α<sub>‖a</sub>)<sup>-1</sup> exp(-v<sub>⊥</sub><sup>2</sup>/α<sub>⊥a</sub><sup>2</sup> - v<sub>‖</sub><sup>2</sup>/α<sub>‖a</sub><sup>2</sup>), where T<sub>⊥a</sub> = m<sub>a</sub>α<sub>⊥a</sub><sup>2</sup>/2 and T<sub>‖a</sub> = m<sub>a</sub>α<sub>‖a</sub><sup>2</sup>/2 are perpendicular and parallel temperatures. Here, the notion of “perpendicular” and “parallel” refers to directions with respect to the ambient magnetic field, which is assumed to lie along the z-axis. Note that the Balescu–Lenard Equation (A1), in and of itself, has no explicit dependence on the ambient magnetic field, but the extension of such an equation to include an explicit magnetic field dependence is available in the literature (Dong et al. 2016, 2023; Yoon 2016). For the sake of simplicity, however, we choose to work with the classical version of the Balescu–Lenard Equation (A1). In any event, taking the velocity

moments,

$$ 
\begin{aligned}
\frac{dn_a T_{\perp a}}{dt} &= \frac{m_a}{2} \int d\mathbf{v} v_{\perp}^2 \frac{\partial f_a(\mathbf{v})}{\partial t}, \\
\frac{dn_a T_{\parallel a}}{dt} &= m_a \int d\mathbf{v} v_{\parallel}^2 \frac{\partial f_a(\mathbf{v})}{\partial t},
\end{aligned}
\quad (A2) $$

and inserting the right-hand side of Equation (A1) into the right-hand side of Equation (A2), with the assumption of a bi-Maxwellian distribution, we arrive at the following result:

$$ 
\begin{aligned}
\frac{dn_a T_{\perp a}}{dt} &= \frac{n_a e_a^2}{\pi^{3/2} m_a \alpha_{\perp a}^2 \alpha_{\parallel a}} \int \frac{d\mathbf{k}}{k^2} \int d\omega \\
&\quad \times \mathrm{Im} \left[ \frac{m_a}{2\pi^3} \frac{\mathrm{Im} \, \epsilon(\mathbf{k}, \omega)}{|\epsilon(\mathbf{k}, \omega)|^2} I_1 - 2 \langle \delta E^2 \rangle_{k, \omega} \left( \frac{I_2}{\alpha_{\perp a}^2} + \frac{I_3}{\alpha_{\parallel a}^2} \right) \right], \\
\frac{dn_a T_{\parallel a}}{dt} &= \frac{2 n_a e_a^2}{\pi^{3/2} m_a \alpha_{\perp a}^2 \alpha_{\parallel a}} \int \frac{d\mathbf{k}}{k^2} \int d\omega \\
&\quad \times \mathrm{Im} \left[ \frac{m_a}{2\pi^3} \frac{\mathrm{Im} \, \epsilon(\mathbf{k}, \omega)}{|\epsilon(\mathbf{k}, \omega)|^2} I_4 - 2 \langle \delta E^2 \rangle_{k, \omega} \left( \frac{I_3}{\alpha_{\perp a}^2} + \frac{I_5}{\alpha_{\parallel a}^2} \right) \right],
\end{aligned}
\quad (A3) $$

where

$$ 
\begin{pmatrix} I_1 \\ I_2 \\ I_3 \\ I_4 \\ I_5 \end{pmatrix} = \int d\mathbf{v} \begin{pmatrix} k_{\perp} v_{\perp} \\ k_{\perp}^2 v_{\perp}^2 \\ k_{\parallel} k_{\parallel} v_{\parallel} \\ k_{\parallel} v_{\parallel} \\ k_{\parallel}^2 v_{\parallel}^2 \end{pmatrix} \frac{\exp(-v_{\perp}^2 / \alpha_{\perp a}^2 - v_{\parallel}^2 / \alpha_{\parallel a}^2)}{\omega - k_{\perp} v_{\perp} - k_{\parallel} v_{\parallel} + i0}. \quad (A4) $$

The velocity integrals can be carried out, and the result is expressed in terms of the plasma dispersion function,

$$ 
\begin{aligned}
I_1 &= \frac{\pi^{3/2} \alpha_{\perp a}^2 \alpha_{\parallel a}}{2} \frac{k_{\perp}^2 \alpha_{\perp a}^2}{\omega^2} \zeta_a^2 Z'(\zeta_a), \\
I_2 &= -\frac{\pi^{3/2} \alpha_{\perp a}^2 \alpha_{\parallel a}}{2} \frac{k_{\perp}^2 \alpha_{\perp a}^2}{\omega} \zeta_a \left( Z(\zeta_a) + \frac{k_{\perp}^2 \alpha_{\perp a}^2}{2\omega^2} \zeta_a^2 Z''(\zeta_a) \right), \\
I_3 &= -\frac{\pi^{3/2} \alpha_{\perp a}^2 \alpha_{\parallel a}}{2} \frac{k_{\perp}^2 k_{\parallel}^2 \alpha_{\perp a}^2 \alpha_{\parallel a}^2}{2\omega^3} \zeta_a^3 Z''(\zeta_a), \\
I_4 &= \frac{\pi^{3/2} \alpha_{\perp a}^2 \alpha_{\parallel a}}{2} \frac{k_{\parallel}^2 \alpha_{\parallel a}^2}{\omega^2} \zeta_a^2 Z'(\zeta_a), \\
I_5 &= -\frac{\pi^{3/2} \alpha_{\perp a}^2 \alpha_{\parallel a}}{2} \frac{k_{\parallel}^2 \alpha_{\parallel a}^2}{\omega} \zeta_a \left( Z(\zeta_a) + \frac{k_{\parallel}^2 \alpha_{\parallel a}^2}{2\omega^2} \zeta_a^2 Z''(\zeta_a) \right),
\end{aligned}
\quad (A5) $$

where $\zeta_a = \omega/(k_{\perp}^2 \alpha_{\perp a}^2 + k_{\parallel}^2 \alpha_{\parallel a}^2)^{1/2}$ and $Z(\zeta) = \pi^{-1/2} \int_{-\infty}^{\infty} dx e^{-x^2}(x-\zeta)^{-1}$ is the plasma dispersion (or Fried–Conte) function, the prime denoting the derivative with respect to the argument.
Substituting this result into the temperature Equations (A3), we have

$$
\begin{align}
\frac{dT_{\perp a}}{dt} = \frac{\pi^{1/2} e_a^2}{m_a} &\int dk \frac{k_{\perp}^2}{k^2} \int \frac{d\omega}{\omega} \left[ -\frac{\operatorname{Im} \epsilon(\mathbf{k}, \omega)}{\pi^3 |\epsilon(\mathbf{k}, \omega)|^2} \frac{\zeta_a^2}{\omega} T_{\perp a} \right. \nonumber \\
&\left. + \langle \delta E^2 \rangle_{\mathbf{k}, \omega} \left( 1 - \frac{k^2 \alpha_{\perp a}^2}{\omega^2} \zeta_a^2 (1 - 2\zeta_a^2) \right) \right] \zeta_a e^{-\zeta_a^2}, \nonumber \\
\frac{dT_{\parallel a}}{dt} = \frac{2\pi^{1/2} e_a^2}{m_a} &\int dk \frac{k_{\parallel}^2}{k^2} \int \frac{d\omega}{\omega} \left[ -\frac{\operatorname{Im} \epsilon(k, \omega)}{\pi^3 |\epsilon(k, \omega)|^2} \frac{\zeta_a^2}{\omega} T_{\parallel a} \right. \nonumber \\
&\left. + \langle \delta E^2 \rangle_{\mathbf{k}, \omega} \left( 1 - \frac{k^2 \alpha_{\parallel a}^2}{\omega^2} \zeta_a^2 (1 - 2\zeta_a^2) \right) \right] \zeta_a e^{-\zeta_a^2}, \tag{A6}
\end{align}
$$

where we have made use of $Z'(\zeta) = -2 - 2\zeta Z(\zeta)$, $Z''(\zeta) = 4\zeta - 2Z(\zeta) + 4\zeta^2 Z(\zeta)$, and $\operatorname{Im} Z(\zeta_a) = \pi^{1/2} e^{-\zeta_a^2}$.

The quantities of relevance, namely, $\langle \delta E^2 \rangle_{\mathbf{k},\omega}$ and $\epsilon(\mathbf{k}, \omega)$, can easily be computed on the basis of a bi-Maxwellian distribution. The intermediate steps are straightforward and thus omitted. The final result is

$$
\begin{align}
\langle \delta E^2 \rangle_{\mathbf{k}, \omega} &= \frac{1}{k^2 \omega |\epsilon(\mathbf{k}, \omega)|^2} \sum_b \frac{2n_b e_b^2}{\pi^{3/2}} \zeta_b e^{-\zeta_b^2}, \nonumber \\
\epsilon(\mathbf{k}, \omega) &= 1 - \sum_b \frac{4\pi n_b e_b^2}{m_b \omega^2} \zeta_b^2 Z'(\zeta_b), \tag{A7}
\end{align}
$$

where $\zeta_b$ is the same as $\zeta_a$, except that the charged-particle species label $a$ is to be replaced by $b$. By combining various terms that we have derived, we may explicitly write the temperature equations as follows:

$$
\begin{align}
\frac{dT_{\perp a}}{dt} &= \sum_b \frac{2n_b e_a^2 e_b^2}{\pi m_a} \int dk \frac{k_{\perp}^2}{k^2} \int \frac{d\omega}{\omega^2} \frac{1}{|\epsilon(\mathbf{k}, \omega)|^2} \nonumber \\
&\quad \times \left( \frac{1}{k^2} - \frac{2T_{\perp a}}{m_a} \frac{\zeta_a^2}{\omega^2} + \frac{4T_{\perp a}}{m_a} \frac{\zeta_a^4}{\omega^2} - \frac{4T_{\perp a}}{m_b} \frac{\zeta_a^2 \zeta_b^2}{\omega^2} \right) \zeta_a \zeta_b e^{-\zeta_a^2 - \zeta_b^2}, \nonumber \\
\frac{dT_{\parallel a}}{dt} &= \sum_b \frac{4n_b e_a^2 e_b^2}{\pi m_a} \int dk \frac{k_{\parallel}^2}{k^2} \int \frac{d\omega}{\omega^2} \frac{1}{|\epsilon(\mathbf{k}, \omega)|^2} \nonumber \\
&\quad \times \left( \frac{1}{k^2} - \frac{2T_{\parallel a}}{m_a} \frac{\zeta_a^2}{\omega^2} + \frac{4T_{\parallel a}}{m_a} \frac{\zeta_a^4}{\omega^2} - \frac{4T_{\parallel a}}{m_b} \frac{\zeta_a^2 \zeta_b^2}{\omega^2} \right) \zeta_a \zeta_b e^{-\zeta_a^2 - \zeta_b^2}. \tag{A8}
\end{align}
$$

At this point, we note that if we replace $\epsilon(\mathbf{k}, \omega)$ by unity—which amounts to the approximation that leads to the Landau collisional integral—then the $\omega$ integration in Equation (A8) can be performed exactly. In fact, Hellinger & Trávníček (2009) and Jubeh & Barghouthi (2017) have done just that, and obtained the result in terms of generalized hypergeometric functions. However, in this case, the $k$ integral diverges at both ends, $k=0$ and $k=\infty$. Hence, the upper and lower limits must be replaced by the cutoffs prescribed by Landau. This leads to the Coulomb logarithm approximation. For a general situation, the $\omega$ integral cannot be carried out exactly. If we suppose, however, that the Z function that appears within the definition for $\epsilon(\mathbf{k}, \omega)$ can be approximated by a small argument expansion, then further progress can be made. This is a

reasonable assumption since the $\omega$ integral has the largest contribution from the $\omega \sim 0$ region as it is weighted by the exponential factor $e^{-(\zeta_a^2+\zeta_b^2)}$, which peaks at $\omega=0$. If we make such an assumption then we may approximate $Z(\zeta_a) \sim 0$, which leads to

$$
\epsilon(\mathbf{k}, \omega) \approx 1 + \frac{1}{k^2 \lambda_D^2}, \qquad \lambda_D^2 = \frac{T_e}{4\pi n e^2}. \tag{A9}
$$

Then we have

$$
\begin{align}
\frac{dT_{\perp a}}{dt} &= \sum_b \frac{4n_b e_a^2 e_b^2 s_{ab}^{1/2}}{\pi^{1/2} m_a (\alpha_{\parallel a}^2 + \alpha_{\parallel b}^2)^{3/2}} \nonumber \\
&\quad \times \int dk \frac{\lambda_D^4}{(1 + k^2 \lambda_D^2)^2} \frac{k_{\perp}^2}{(k_{\perp}^2 + s_{ab} k_{\parallel}^2)^{1/2}} \nonumber \\
&\quad \times \left( \frac{T_{\parallel a}}{m_a} - \frac{T_{\perp a}}{m_a} \frac{s_{ab} k^2}{k_{\perp}^2 + s_{ab} k_{\parallel}^2} + \frac{T_{\parallel b}}{m_b} - \frac{T_{\perp b}}{m_b} \frac{s_{ab} k^2}{k_{\perp}^2 + s_{ab} k_{\parallel}^2} \right), \nonumber \\
\frac{dT_{\parallel a}}{dt} &= \sum_b \frac{8n_b e_a^2 e_b^2 s_{ab}^{1/2}}{\pi^{1/2} m_a (\alpha_{\parallel a}^2 + \alpha_{\parallel b}^2)^{3/2}} \nonumber \\
&\quad \times \int dk \frac{\lambda_D^4}{(1 + k^2 \lambda_D^2)^2} \frac{k_{\parallel}^2}{(k_{\parallel}^2 + s_{ab} k_{\perp}^2)^{1/2}} \nonumber \\
&\quad \times \left( \frac{T_{\parallel a}}{m_a} - \frac{T_{\perp a}}{m_a} \frac{s_{ab} k^2}{k_{\perp}^2 + s_{ab} k_{\parallel}^2} + \frac{T_{\parallel b}}{m_b} - \frac{T_{\perp b}}{m_b} \frac{s_{ab} k^2}{k_{\perp}^2 + s_{ab} k_{\parallel}^2} \right), \nonumber \\
&\qquad s_{ab} = \frac{\alpha_{\parallel a}^2 + \alpha_{\parallel b}^2}{\alpha_{\perp a}^2 + \alpha_{\perp b}^2}. \tag{A10}
\end{align}
$$

We write the $k$ integral in the spherical coordinate system, $\int dk = 2\pi \int_0^\infty dk \, k^2 \int_{-1}^1 dx$, where $x = \cos\theta$ is the cosine of the polar angle. We also express $k_\perp = k(1-x^2)^{1/2}$ and $k_\parallel = kx$. We then carry out the $k$ integral by replacing the upper limit by $k_{\max} = T_e/e^2 = 4\pi n \lambda_D^2$, that is, $\int_0^\infty dk \to \int_0^{4\pi n \lambda_D^2} dk$. Then we obtain

$$
\begin{align}
\int_0^\infty dk \frac{k^3 \lambda_D^4}{(1 + k^2 \lambda_D^2)^2} &\to \int_0^{4\pi n \lambda_D^2} dk \frac{k^3 \lambda_D^4}{(1 + k^2 \lambda_D^2)^2} \nonumber \\
&= \frac{1}{2} \left( \ln(1 + \Lambda^2) - \frac{\Lambda^2}{1 + \Lambda^2} \right) \approx \ln\Lambda, \tag{A11}
\end{align}
$$

where $\Lambda = 4\pi n \lambda_D^3$. Note that, unlike the Landau collision integral, one does not have to replace the lower limit with a finite cutoff value.

In what follows, let us focus on the proton–proton collisions. Using $T_\perp$ and $T_\parallel$ to represent the proton perpendicular and parallel temperatures (that is, we drop the species label $a, b = p$ for the sake of simplicity), we arrive at

$$
\begin{align}
\frac{dT_{\perp}}{dt} &= \frac{2\pi^{1/2}ne^4 \ln\Lambda}{m^{1/2}T_{\perp}^{3/2}} F(T_{\parallel} - T_{\perp}), \nonumber \\
\frac{dT_{\parallel}}{dt} &= \frac{4\pi^{1/2}ne^4 \ln\Lambda}{m^{1/2}T_{\perp}^{3/2}} F(T_{\perp} - T_{\parallel}), \nonumber \\
F &= \int_{-1}^{1} dx \frac{(1-x^2)x^2}{[1+(T_{\parallel}/T_{\perp}-1)x^2]^{3/2}}. \tag{A12}
\end{align}
$$
The integral *F*, which is the same quantity as defined in Equation (2), can be obtained in closed form, and it is as given specifically in Equations (2)–(4); this was first obtained by Kogan (1961) and also found in the NRL Plasma Formulary. This concludes the derivation of the proton–proton collision frequency.

## Appendix B

### Exact Ensemble Calculation of Collisional Relaxation

The main body of the present paper adopted the approximate collision frequency (6), which is based upon the approximation of function *F* by 4/15. Figure 4 is the result of ensemble calculation based on Equation (6), and we have noted therewith that the result is not too different from that based on an exact definition for *F*, as specified in Equations (2)–(4). To concretely demonstrate this, we display in Figure 8 the actual ensemble calculation based on exact *F*. As Figure 8 shows, the result is quite similar to that of Figure 4. One might wonder why we adopted the approximate collision frequency in the first place since the exact solution is available. Aside from the approximate form of collision frequency being simple, thereby facilitating the numerical analysis, the more important reason is for future reference. As part of the ongoing project, the present analysis will be extended to charged-particle species of multiple components in the solar wind. For such a problem, the exact collision frequency becomes too cumbersome. The present proton–proton collision problem, on the other hand, provides the justification for employing the approximated collisional transport and relaxation theory for such a more complex problem in the future.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** The exact ensemble calculation, which can be contrasted to Figure 4, which is based on the approximate collision frequency. As one may appreciate, the result is quite similar to that of Figure 4.</figcaption></figure>

## ORCID iDs

P. H. Yoon https://orcid.org/0000-0001-8134-3790

M. Lazar https://orcid.org/0000-0002-8508-5466

C. Salem https://orcid.org/0000-0002-6536-1531

J. Seough https://orcid.org/0000-0002-1723-5944

M. M. Martinović https://orcid.org/0000-0002-7365-0472

K. G. Klein https://orcid.org/0000-0001-6038-1923

R. A. López https://orcid.org/0000-0003-3223-1498

## References

Alterman, B. L., Kasper, J. C., Stevens, M. L., & Koval, A. 2018, ApJ, 864, 112

Anderson, B. J., Fuselier, S. A., & Murr, D. 1991, GeoRL, 18, 1955

Arzamasskiy, L., Kunz, M., Squire, J., Quataert, E., & Schekochihin, A. A. 2023, PhRvX, 13, 021014

Bale, S. D., Drake, J. F., McManus, M. D., et al. 2023, Natur, 618, 252

Bale, S. D., Kasper, J. C., Howes, G. G., et al. 2009, PhRvL, 103, 211101

Balescu, R. 1960, PhFl, 3, 52

Bott, A. F. A., Arzamasskiy, L., Kunz, M. W., Quataert, E., & Squire, J. 2021, ApJL, 922, L35

Chandran, B. D. G., Dennis, T. J., Quataert, E., & Bale, S. D. 2011, ApJ, 743, 197

Chen, C. H. K., Matteini, L., Schekochihin, A. A., Stevens, M. L., et al. 2016, ApJL, 825, L26

Cranmer, S. R. 2014, ApJS, 213, 16

Dong, C., Zhang, W., Cao, J., & Li, D. 2023, RvMPP, 7, 19

Dong, C., Zhang, W., & Li, D. 2016, PhPl, 23, 082105

Fox, C., Velli, M. C., & Bale, S. D. 2016, SSRv, 204, 7

Gary, S. P., Fuselier, S. A., & Anderson, B. J. 1993a, JGR, 98, 1481

Gary, S. P., McKean, M. E., & Winske, D. 1993b, JGR, 98, 3963

Gary, S. P., Skoug, R. M., Steinberg, J. T., & Smith, C. W. 2001, GeoRL, 28, 2759

Gary, S. P., Wang, J., Winske, D., & Fuselier, S. A. 1997, JGR, 102, 27159

He, J., Tu, C., March, E., Bourouaine, S., & Pei, Z. 2013, ApJ, 773, 72

Hellinger, P., & Trávniček, P. M. 2009, PhPl, 16, 054501

Hellinger, P., & Trávniček, P. M. 2010, JCoPh, 229, 5432

Hellinger, P., & Trávniček, P. M. 2014, ApJL, 784, L15

Hellinger, P., Trávniček, P., Kasper, J. C., & Lazarus, A. J. 2006, GeoRL, 33, L09101

Huang, J., Kasper, J. C., Vech, D., et al. 2020, ApJS, 246, 70

Ibscher, D., Lazar, M., Michno, M. J., & Schlickeiser, R. 2013, PhPl, 20, 012103

Ibscher, R., & Schlickeiser, R. 2013, PhPl, 20, 042121

Ibscher, R., & Schlickeiser, R. 2014, PhPl, 21, 022110

Jubeh, W., & Barghouthi, I. A. 2017, PhPl, 24, 122104

Kasper, J. C., Klein, K. G., Weber, T., et al. 2017, ApJ, 849, 126

Kasper, J. C., Lazarus, A. J., Gary, S. P., & Szabo, A. 2003, in AIP Conf. Proc. 679, Solar Wind Ten, ed. M. Velli et al. (Melville, NY: AIP), 538

Kasper, J. C., Lazarus, A. J., & Gary, S. P. 2002, GeoRL, 29, 1839

Kim, S., Ha, J.-H., Ryu, D., & Kang, H. 2021, ApJ, 913, 35

Klein, K. G., Alterman, B. L., Stevens, M. L., Vech, D., & Kasper, J. C. 2018, PhRvL, 120, 205102

Klein, K. G., Martinović, M., Standsby, D., & Horbury, T. S. 2019, ApJ, 887, 234

Klein, K. G., & Veck, D. 2019, RNAAS, 3, 107

Klein, K. G., Verniero, J. L., Alterman, B., et al. 2021, ApJ, 909, 7

Kogan, V. I. 1961, in Plasma Physics and the Problem of Controlled Thermonuclear Reactions, ed. M. A. Leontovich (New York: Pergamon), 153

Landau, L. D. 1936, Phys. Z. Sowjetunion, 10, 154

Lenard, A. 1960, AnPhy, 10, 390

López, R. A., Yoon, P. H., Viñas, A. F., & Lazar, M. 2023, ApJ, 954, 191

Marsch, E., Ao, X. Z., & Tu, C. Y. 2004, JGRA, 109, A04102

Marsch, E., Zhao, L., & Tu, C. Y. 2006, AnGeo, 24, 2057

Martinović, M. M., Klein, K. G., Durovcová, T., & Alterman, B. L. 2021, ApJ, 923, 116
Maruca, B. A., & Kasper, J. C. 2013, AdSpR, 52, 723

Maruca, B. A., Kasper, J. C., & Bale, S. D. 2011, PhRvL, 107, 201101

Maruca, B. A., Kasper, J. C., & Gary, S. P. 2012, ApJ, 748, 137

Matteini, L., Hellinger, P., Goldstein, B. E., et al. 2013, JGRA, 118, 2771

Matteini, L., Hellinger, P., Landi, S., Trávníček, P. M., & Velli, M. 2012, SSRv, 172, 373

Matteini, L., Landi, S., Hellinger, P., et al. 2007, GeoRL, 34, L20105

Mostafavi, P., Allen, R. C., Jagarlamudi, V. K., et al. 2024, A&A, 682, A152

Müller, D., St. Cyr, O. C., Zouganelis, I., et al. 2020, A&A, 642, A1

Ogilvie, K. W., Chornay, D. J., Fritzenreiter, R. J., et al. 1995, SSRv, 71, 55

Opie, S., Verscharen, D., Chen, C. H. K., Owen, C. J., & Isenberg, P. A. 2022, ApJ, 941, 176

Ozak, N., Ofman, L., & Viñas, A. F. 2015, ApJ, 799, 77

Perrone, D., Stansby, D., Horbury, T. S., & Matteini, L. 2019, MNRAS, 483, 3730

Phan, T.-D., Paschmann, G., Baumjohann, W., Sckopke, N., & Lühr, H. 1994, JGR, 99, 121

Raouafi, N. E., Matteini, L., Squire, J., et al. 2023, SSRv, 219, 8

Richardson, A. S. 2019, NRL Plasma Formulary (Washington, DC: Naval Research Laboratory), 33

Santos-Lima, R., Yan, H., de Gouveia Dal Pino, E. M., & Lazarian, A. 2016, MNRAS, 460, 2492

Sarfraz, M., López, R. A., Ahmed, S., & Yoon, P. H. 2022, MNRAS, 509, 3764

Seough, J., & Yoon, P. H. 2012, JGR, 117, A08101

Seough, J., Yoon, P. H., & Hwang, J. 2014, PhPl, 21, 062118

Seough, J., Yoon, P. H., & Hwang, J. 2015, PhPl, 22, 012303

Seough, J., Yoon, P. H., Kim, K.-H., & Lee, D.-H. 2013, PhRvL, 110, 071103

Seough, J., Yoon, P. H., Nariyuki, Y., & Salem, C. 2023, ApJ, 953, 8

Servidio, S., Osman, K. T., Valentini, F., et al. 2014, ApJL, 781, L27

St-Onge, D. A., & Kunz, M. W. 2018, ApJL, 863, L25

Štverák, Š., Trávníček, P., Maksimovic, M., et al. 2008, JGRA, 113, A03103

Sun, H., Zhao, J., Liu, W., Xie, H., & Wu, D. 2020, ApJ, 902, 59

Sun, H., Zhao, J., Xie, H., & Wu, D. 2019, ApJ, 884, 44

Tan, L. C., Fung, S. F., Kessel, R. L., et al. 1998, GeoRL, 25, 587

Vafin, S., Riazantseva, M., & Pohl, M. 2019, ApJL, 871, L11

Verscharen, D., Klein, K. G., & Maruca, B. A. 2019, LRSP, 16, 5

Wilson, L. B., III, Stevens, M. L., Kasper, J. C., et al. 2018, ApJS, 236, 41

Xiang, L., Lee, K. H., Wu, D. J., Yu, H. W., & Lee, L. C. 2021, ApJ, 916, 30

Yoon, P. H. 2016, JGRA, 121, 665

Yoon, P. H. 2017, RvMPP, 1, 4

Yoon, P. H., López, R. A., & Zaheer, S. 2023, ApJ, 950, 131

Yoon, P. H., & Sarfraz, M. 2017, ApJ, 835, 246

Yoon, P. H., & Seough, J. 2012, JGRA, 117, A08102

Yoon, P. H., Seough, J., Salem, C. S., & Klein, K. G. 2019, PhRvL, 123, 145101

Zhao, G. Q., Li, H., Feng, H. Q., et al. 2019, ApJ, 884, 60