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

# A Parametric Study of Solar Wind Properties and Composition Using Fluid and Kinetic Solar Wind Models

Paul Lomazzi<fnref n="1" />, Alexis P. Rouillard<fnref n="1" />, Michael A. Lavarra<fnref n="1" />, Nicolas Poirier<fnref n="2" />, Pierre-Louis Blelly<fnref n="1" />, Jean-Baptiste Dakeyo<fnref n="3" />, Viviane Pierrard<fnref n="4" />,<fnref n="5" />, Victor Réville<fnref n="1" />, Christian Vocks<fnref n="6" />, and Simon Thomas<fnref n="1" />

## Abstract

The physical processes in the solar corona that shape the solar wind remain an active research topic. Modeling efforts have shown that energy and plasma exchanges near the transition region play a crucial role in modulating solar wind properties. Although these regions cannot be measured in situ, plasma parameters can be inferred from coronal spectroscopy and ionization states of heavy ions, which remain unchanged as they escape the corona. We introduce a new solar wind model extending from the chromosphere to the inner heliosphere, capturing thermodynamic coupling across atmospheric layers. By including neutral and charged particle interactions, we model the transport and ionization processes of the gas through the transition region and corona and into the solar wind. Instead of explicitly modeling coronal heating, we link its spatial distribution to large-scale magnetic field properties. Our results confirm that energy deposition strongly affects wind properties through key mechanisms involving chromospheric evaporation, thermal expansion, and magnetic flux expansion. For sources near active regions, the model predicts significant solar wind acceleration, with plasma outflows comparable to those inferred from coronal spectroscopy. For winds from large coronal holes, the model reproduces the observed anticorrelation between charge state and wind speed. However, the predicted charge state ratios are overall lower than observed. Inclusion of a population of energetic electrons enhances both heavy ion charge states and solar wind acceleration, improving agreement with observations.

*Unified Astronomy Thesaurus concepts: Solar wind (1534); Solar corona (1483); Solar chromosphere (1479); Ionization (2068); Astronomical simulations (1857)*

## 1. Introduction

The mechanisms in the corona that control the bulk properties of the solar wind are still under debate. Previous theoretical studies have provided evidence that energy and plasma flows through the transition region play a central role in the mass flux escaping the solar corona (e.g., V. H. Hansteen & E. Leer 1995). These energy exchanges in turn depend on the heating processes at work in the lower solar atmosphere, which are poorly understood.

The source region of the solar wind plasma is still collisional in the corona, transitioning to a collisionless state at ~1 $R_⊙$ above the photosphere. In the collisional region, heavy ions undergo collisions with hot electrons and become increasingly ionized with height, and this continues until electron densities are so low and their temperatures so high that collisional processes vanish and the charge state of an ion can no longer change and becomes “frozen in” (J. Geiss et al. 1995).

Ion charge states measured remotely or in situ therefore provide information on solar wind source conditions and, in particular, electron source temperatures. Previous observational studies using data from Ulysses and the NASA Advanced Composition Explorer (ACE; E. C. Stone et al. 1998) have shown that the charge states in fast and slow solar wind streams

differ greatly (J. Geiss et al. 1995; L. Zhao et al. 2009). In particular, they showed that the charge state ratios (CSRs) $O_7^+/O_6^+$ and $C_6^+/C_5^+$ tend to be higher in slow winds than in fast winds, suggesting a hotter coronal origin of the slow wind.

While coronal holes are widely accepted as the source region of the Alfvénic fast wind, the source location of slow wind is still under debate and is thought to have multiple source components including dynamic solar streamers (e.g., A. P. Rouillard et al. 2009; A. Rouillard et al. 2010), boundaries of coronal holes (Y.-M. Wang et al. 2009), and low-latitude (LL) coronal holes (R. D'Amicis & R. Bruno 2015; Y.-M. Wang & Y.-K. Ko 2019). The temperatures of the dynamic sources of the slow solar wind originating in streamers are complex to grasp because they involve periodic plasma exchanges between loops and the open field (e.g., V. Réville et al. 2022).

Assuming that the fast and slow wind components originate in the more quiescent coronal holes, previous theoretical studies based on fluid wind models have tried to reproduce the anticorrelation between wind speed and charge state (S. R. Cranmer et al. 2007; S. R. Cranmer 2014; R. Oran et al. 2015). Although they have been successful in obtaining an anticorrelation, the resulting charge states are orders of magnitude lower than those measured in situ.

The sources of these discrepancies are still unknown. Y.-K. Ko et al. (1997) and R. Esser & R. J. Edgar (2001) proposed that differential flows between minor and major species could lead to increases in the charge state. A extra population of ionizing suprathermal electrons was also invoked

<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](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Université de Toulouse, UPS-OMP; IRAP; Toulouse, France CNRS IRAP; 9 Av. colonel Roche, BP 44346, France; arouillard@irap.omp.eu</footnote>

<footnote n="2">² Rosseland Centre for Solar Physics, University of Oslo, P.O. Box 1029, Blindern, NO-0315 Oslo, Norway</footnote>

<footnote n="3">³ Space Sciences Laboratory, University of California, Berkeley, CA, USA</footnote>

<footnote n="4">⁴ Royal Belgian Institute for Space Aeronomy, Space Physics and STCE, 3 av. circulaire, B-1180 Brussels, Belgium</footnote>

<footnote n="5">⁵ Université Catholique de Louvain, ELI-C, 3 place Louis Pasteur, B-1348 Louvain-la-Neuve, Belgium</footnote>

<footnote n="6">⁶ Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany</footnote>

<footnote>Received 2025 February 14; revised 2025 March 17; accepted 2025 March 17; published 2025 May 9</footnote>
to explain the ionization rates of minor ions (Y.-K. Ko et al. 1996; R. Esser & R. J. Edgar 2000; S. R. Cranmer 2014; R. Oran et al. 2015).

An accurate model of the ion charge states in the solar wind requires a detailed transport model not only of thermal particles accounting for the thermodynamic coupling of the different layers of the low solar atmosphere but also of the transport of nonthermal particles and their interactions with the thermal population. This is a challenging task that ideally should account for the detailed evolution of the particle distribution function of all species.

In the present paper, we use a more tractable approach by using the IRAP Solar Atmospheric Model (ISAM) developed by M. Lavarra (2022) in combination with an exospheric kinetic solar wind model (V. Pierrard et al. 2023) to study the different processes that could influence the ionization states of heavy ions in the different types of solar winds. This approach allows us to model the ionizing effects of energy deposition, nonthermal particles, and differential flows for different conditions at the sources of the fast and slow solar winds. The model is constrained by several observations, including SOlar and Heliospheric Observatory (SOHO; V. Domingo et al. 1995) spectroscopy and in situ data from the Solar Orbiter (D. Müller et al. 2020) and Parker Solar Probe (PSP; J. Kinnison et al. 2020) missions.

This paper is organized as follows. After a presentation of the models used in Section 2, we begin by solving the transport of neutral and charged particles from the upper chromosphere to the solar wind in Sections 4 and 6. For this, we use the results of the recent solar magnetic connectivity study of J.-B. Dakeyo et al. (2024) to specify the input parameters necessary to model the fast and slow solar winds in Section 3. This provides first estimates of the charge state of ions for different corona-wind conditions. In Section 7, the thermal model is used in conjunction with an exospheric solar wind model that provides the radial evolution of nonthermal electrons. This allows us to model the additional ionizing effects on ions of these nonthermal electrons.

## 2. Description of the Numerical Model

### 2.1. The Conservation Equations

To model the thermal component of the coupled chromosphere-corona-wind system, we use ISAM, which dynamically solves the transport and interactions of neutrals and charged particles. The model uses a 16-moment set of transport equations described in detail in P. L. Blelly & W. Schunk (1993). A first version of this model, the IRAP Plasmasphere-Ionosphere Model, was developed to simulate plasma transport in Earth's upper atmosphere (A. Marchaudon & P.-L. Blelly 2015) and has been thoroughly checked against a broad range of terrestrial observations (A. Marchaudon & P. L. Blelly 2020) and adapted to other planetary atmospheres (P. L. Blelly et al. 2019). It was later modified to model the solar atmosphere by solving neutrals together with charged particles to simulate the variable low solar atmosphere for open field configurations (M. Lavarra 2022) and magnetic loops (C. Poirier 2022).

Other implementations of the 16-moment approach to model the solar atmosphere are detailed in X. Li (1999) and later ø. Lie-Svendsen et al. (2001), Å. M. Janse et al. (2007), and H. S. Byhring et al. (2011). A comparison between ISAM and

these earlier applications of the 16-moment models is discussed in M. Lavarra (2022).

The density $n_s$, bulk flow velocity $u_s$, parallel $T_s^\parallel$ and perpendicular temperatures $T_s^\perp$, and parallel $q_s^\parallel$ and perpendicular heat fluxes $q_s^\perp$ are solved along a given magnetic flux tube through the following set of conservation equations:

$$ \frac{\partial n_s}{\partial t} + u_s \nabla_{\parallel} n_s + \frac{n_s}{A} \nabla_{\parallel} (A u_s) = \frac{\delta n_s}{\delta t}, \quad (1) $$

$$ \frac{\partial u_s}{\partial t} + u_s \nabla_{\parallel} u_s + \frac{\nabla_{\parallel} n_s k_B T_s^\parallel}{n_s m_s} + \frac{k_B}{m_s} (T_s^\parallel - T_s^\perp) \frac{\nabla_{\parallel} A}{A} + \frac{GM_{\odot}}{r^2} \cos(\theta) - \frac{F_s}{n_s m_s} = \frac{\delta u_s}{\delta t}, \quad (2) $$

$$ \frac{\partial T_s^\parallel}{\partial t} + u_s \nabla_{\parallel} T_s^\parallel + 2T_s^\parallel \nabla_{\parallel} u_s + \frac{\nabla_{\parallel} q_s^\parallel}{n_s k_B} + \frac{q_s^\parallel - 2q_s^\perp}{n_s k_B} \frac{\nabla_{\parallel} A}{A} = \frac{(Q_s^\parallel - \Lambda_s^\parallel)}{n_s k_B} + \frac{\delta T_s^\parallel}{\delta t}, \quad (3) $$

$$ \frac{\partial T_s^\perp}{\partial t} + u_s \nabla_{\parallel} T_s^\perp + \frac{\nabla_{\parallel} q_s^\perp}{n_s k_B} + \left( \frac{2q_s^\perp}{n_s k_B} + u_s T_s^\perp \right) \frac{\nabla_{\parallel} A}{A} = \frac{(Q_s^\perp - \Lambda_s^\perp)}{n_s k_B} + \frac{\delta T_s^\perp}{\delta t} \quad (4) $$

$$ \frac{\partial q_s^\parallel}{\partial t} + u_s \nabla_{\parallel} q_s^\parallel + 4q_s^\parallel \nabla_{\parallel} u_s + \frac{3n_s k_B^2 T_s^\parallel}{m_s} \nabla_{\parallel} T_s^\parallel + u_s q_s^\parallel \frac{\nabla_{\parallel} A}{A} = \frac{\delta q_s^\parallel}{\delta t}, \quad (5) $$

$$ \frac{\partial q_s^\perp}{\partial t} + u_s \nabla_{\parallel} q_s^\perp + 2q_s^\perp \nabla_{\parallel} u_s + \frac{n_s k_B^2 T_s^\parallel}{m_s} \nabla_{\parallel} T_s^\perp + \left( \frac{n_s k_B^2 T_s^\perp}{m_s} (T_s^\parallel - T_s^\perp) + 2u_s q_s^\perp \right) \frac{\nabla_{\parallel} A}{A} = \frac{\delta q_s^\perp}{\delta t}. \quad (6) $$

Here, the subscript $s$ indicates the species considered, $k_B$ is the Boltzmann constant, $m_s$ is the atomic mass of the species, $G$ is the gravitational constant, $M_\odot$ is the solar mass, $\theta$ is the angle between the radial direction and the local direction of the magnetic field line along which the equations are solved, $\Lambda_s^{(\perp)}$ are the perpendicular and parallel radiative cooling terms, and $Q_s^{(\perp)}$ are the perpendicular and parallel heating terms. The terms on the right-hand side of the equations containing $\delta/\delta t$ represent collisional effects. The parallel and perpendicular heat fluxes represent the transport along the field lines of the parallel and perpendicular energies.

For ions, $F_s$ corresponds to the electrostatic polarization field contribution,

$$ \frac{F_s}{m_s n_s} = - \frac{Z_s}{m_s} \frac{1}{n_e} \left[ \nabla_{\parallel}(n_e k_B T_e^\parallel) + k_B (T_e^\parallel - T_e^\perp) \frac{1}{A} \nabla_{\parallel} A \right] + Z_s \frac{m_e}{m_i} \frac{\delta u_e}{\delta t}, \quad (7) $$

where the subscript $s$ indicates the ion considered and $Z_s$ represents the number of elementary charges $e$ carried by the ion.

Equations (1) and (2) are not solved for the electrons; instead, we impose quasi-neutrality and no field-aligned current
$J=0$, which allows us to have

$$n_e = \sum_{s=\text{ions}} (Z_s n_s), \qquad (8)$$

$$u_e = \frac{\sum_{s=\text{ions}} (Z_s n_s u_s)}{n_e}. \qquad (9)$$

The total temperature and heat flux are related to their parallel and perpendicular components by

$$T = (T_s^{\parallel} + 2T_s^{\perp})/3, \qquad (10)$$

$$q = (q_s^{\parallel} + 2q_s^{\perp})/2. \qquad (11)$$

We begin the study by computing the coupled transport of electrons, protons, and hydrogen atoms along magnetic flux tubes rooted in the chromosphere and extending outward into the solar wind.

## 2.2. The Numerical Scheme

To solve the above coupled differential equations, ISAM uses the flux-corrected algorithm for solving generalized continuity equations (J. P. Boris & D. L. Book 1976; J. P. L. S. Boris et al. 1993). This finite-volume method is based on an explicit second-order conservative Godunov scheme. The scheme is noncentered and stabilized up to fourth order thanks to a numerical antidiffusive stage, but in the case of steep gradients, the scheme is degraded to first-order accuracy by a stabilization procedure characterized by a diffusive numerical flux. ISAM solves the conservative equations including transport and local source terms. The time dependence is solved using a fourth-order Runge-Kutta algorithm. The conservation Equations (1)–(6) are integrated between an inner boundary located in the chromosphere at ~600 km above the photosphere and all the way to $20 R_\odot$.

## 2.3. The Model Setup

### 2.3.1. The Open Magnetic Field Model

The magnetic flux tube geometry along which ISAM solves the plasma transport equations is given by the widely used form inferred by R. A. Kopp & T. E. Holzer (1976) from the analysis of eclipse observations,

$$A = A_0 \left( \frac{r}{R_\odot} \right)^2 \left( \frac{f_{\text{max}} e^{(r-r_g)/\sigma_g} + f_g}{e^{(r-r_g)/\sigma_g} + 1} \right), \qquad (12)$$

where $A(r)$ is the section of the flux tube with $A_0 = A(R_\odot)$. $r_g$ and $\sigma_g$ are parameters describing, respectively, the height where the expansion rate is the strongest and the steepness of the expansion. $f_{\text{max}}$ is the asymptotic value of $(A/A_0)/(R_\odot/r)^2$ describing the deviation from a radial expansion. $f_g$ is chosen to have $A(R_\odot) = A_0$, so $f_g = 1 + (1 - f_{\text{max}})e^{(R_\odot - r_g)/\sigma_g}$. In this study, only the maximum expansion factor $f_{\text{max}}$ is changed, and the other parameters are kept fixed with $\sigma_g = 0.5 R_\odot$ and $r_g = 1.3 R_\odot$. The values for $f_{\text{max}}$ will be given by a global magnetic model of the solar atmosphere in the context of a magnetic connectivity analysis of Solar Orbiter data.

### 2.3.2. The Atmospheric Heating Rates

The heating processes are still highly debated and could have multiple origins in the different layers of the solar

<figure><img src="image_3.png" /><figcaption>Figure 1. Heating function used for the fast wind solution from the model. Heating on the protons and electrons is indicated by using, respectively, dotted and dashed lines. Total heating (black) is the result of transition region heating (red), coronal heating (blue), and solar wind heating (green).</figcaption></figure>

atmosphere. Our approach is to redefine a mathematical form of the energy deposition and to compare our solar wind simulations with a broad range of remote-sensing and in situ observations. The heating terms $Q_s^{(\perp)}$ in Equations (3) and (4) are a sum of three different heating functions acting on the transition region heating $Q_{\text{TR},s}^{(\perp)}$, the corona $Q_{\text{cor},s}^{(\perp)}$, and the solar wind $Q_{\text{sw},s}^{(\perp)}$.

The heating repartition between directions and species is given in Figure 1. For the transition region and coronal heating terms $Q_{\text{TR},s}^{(\perp)}$ and $Q_{\text{cor},s}^{(\perp)}$, we use the same phenomenological form as used in R. F. Pinto & A. P. Rouillard (2017), which was initially introduced by G. L. Withbroe (1988), J. F. McKenzie et al. (1995), and S. Rifai Habbal et al. (1995):

$$Q = -\nabla \cdot F_h, \qquad (13)$$

where

$$F_h = F_{B_0} \frac{A_0}{A} e^{-\frac{r-R_\odot}{H_f}}. \qquad (14)$$

This functional form supposes a total energy flux $F_{B_0}$ proportional to $B_0$ and a heating scale height $H_f$ that controls the altitude range of energy deposition. As in Y.-M. Wang et al. (2009), this heating scale height is made inversely dependent on the maximum expansion factor $f_{\text{max}}$. Although the mechanisms heating the solar atmosphere are not certain, the dissipation of Alfvén waves channeled by the interaction of upward-propagating and reflected counterpropagating waves reproduce a number of coronal properties (S. R. Cranmer et al. 2007; A. Verdini & M. Velli 2007; A. Verdini et al. 2012). In essence, the value of $f_{\text{max}}$ controls the radial gradients in Alfvén speed along a flux tube, which in turn should dictate the level of wave reflection and associated energy dissipation along that tube. A high value of $f_{\text{max}}$ is expected to lower the scale height of energy deposition by concentrating wave dissipation in the lower corona.

The following functional forms of the terms in $Q_{\text{cor},s}^{(\perp)}$ and $Q_{\text{TR},s}^{(\perp)}$ were found to reproduce well the typical height of the transition region and some general features of the solar corona:

$$F_{B_0}^{\text{cor}} = 0.86 \times 10^5 B_0 \text{ erg cm}^{-2} \text{ s}^{-1}, \qquad (15)$$

$$H_f^{\text{cor}} = 4.35 \times f_{\text{ss}}^{-0.75} R_\odot, \qquad (16)$$

$$F_{B_0}^{\text{TR}} = 0.43 \times 10^5 B_0 \text{ erg cm}^{-2} \text{ s}^{-1}. \qquad (17)$$
$H_f^{TR}$ is kept unchanged in all runs, and we chose $H_f^{TR} = 0.01R_\odot$. The proportional relationship between the energy flux and $B_0$ was demonstrated in Y. M. Wang (2016) using OMNI data (J. H. King & N. E. Papitashvili 2005), where the energy flux was derived from the mass and energy conservation equations.

When turbulence reaches dissipative scales, multiple kinetic mechanisms will heat the plasma and distribute the energy between the different species (W. I. Axford & J. F. McKenzie 1992; J. L. Kohl et al. 1998; Y. G. Maneva et al. 2015). It is well established that ion and electron heating is sustained far out into the interplanetary medium (E. Marsch 2006) as evidenced by the elevated electron temperatures measured beyond 1 au (K. Issautier et al. 1998) and the presence of $T_\perp > T_\parallel$ of ions beyond 0.3 au.

The radial distribution of this interplanetary wind heating $Q_{sw,s}^{(\perp)}$ is chosen here to follow a power law starting at the height $r_s^{sw}$, beyond which the coronal heating is replaced by the wind heating. The power law begins at $Q_{cor,s}^{(\perp)}$ at $r_s^{sw}$ with a fixed slope defined by the derivative of $Q_{cor,s}^{(\perp)}$ at that radial distance,

$$Q_{sw} = Q_{cor,s}^{(\perp)}(r_s^{sw}) \left( \frac{r(r_s^{sw}: 20 R_\odot)}{r_s^{sw}} \right)^{(\Delta \log Q_{cor,s}^{(\perp)}/\Delta \log r)|_{r_s^{sw}}}, \quad (18)$$

where $\Delta(*)$ indicates the difference of (*) between two successive grid points.

In this formulation, the amount of heating in the solar wind is only controlled by changing $r_s^{sw}$. This choice was made to keep a smooth transition between the coronal and solar wind heating while still benefiting from the exponential expression of the coronal heating term.

Heating is applied only to electrons in the transition region. This energy is transmitted to other particles through the frequent Coulomb collisions occurring in that region. Coronal and solar wind heating are applied 60% to protons and 40% to electrons to get, as we shall see, reasonable electron temperatures compared to spectroscopic observations. A total of 80% of the proton heating is applied in the perpendicular direction and 20% in the parallel direction. The same ratio between parallel and perpendicular heating is applied to electrons. This heating repartition gave optimal results from many iteration tests and is kept for all simulations presented in this paper.

## 2.4. Radiative Energy Loss

The radiative cooling term, $\Lambda$, is based on the optically thin function of R. G. Athay (1986) with a correction factor to account for optical thickness in the upper chromosphere for temperatures lower than $\sim 2 \times 10^4$ K,

$$\Lambda = n_e(n_p + n_H)\lambda(T), \quad (19)$$

where

$$\lambda(T) = 10^{-21} 10^{[\log_{10}(T/T_M)]^2} \chi(T) \quad (20)$$

and

$$\chi(T) = \begin{cases} 1 & \text{for } T > T_1, \\ \frac{T-T_0}{T_1-T_0} & \text{for } T_0 < T < T_1, \\ 0 & \text{for } T < T_0. \end{cases}$$

$T_0, T_1$, and $T_M$ are fitting parameters adjusted to simulate radiative losses and medium opacity. In the chromosphere, our

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Scatter plots of the photospheric magnetic field strength $B_0$ measured by Solar Orbiter and the expansion factor $f_{ss}$ computed from a PFSS model (J.-B. Dakeyo et al. 2024). The data shown cover from 2020 August 1 to 2022 March 17, corresponding to the end of the solar minima. The values of the associated solar wind speeds at Solar Orbiter are displayed by the color code defined by the color bar. Panel (a): $B_0 - f_{ss}$ relations for connectivity points rooted at HLs (>45°). The black line corresponds to the relation $B_0 = f(f_{ss})$ used for our HL simulation series (Equation 21). Panel (b): the $B_0 - f_{ss}$ relations for connectivity points rooted at LLs (<45°). The black line corresponds to the relation $B_0 = f(f_{ss})$ used for our LL simulation series (Equation 22). For all panels, the black crosses correspond to the $B_0$ associated with our $f_{max}$ sample.</figcaption></figure>

formulation of radiative losses ensures that the temperature cannot exceed a threshold value defined here by $T_0 = 6500$ K. $T_1$ and $T_M$ define the temperature range of the transition region. We choose $T_1 = 2 \times 10^4$ K and $T_M = 2 \times 10^5$ K.

## 3. Constraining Free Parameters with Solar Orbiter Data

The first solar encounters of the PSP and Solar Orbiter missions have offered an exceptional context for testing ISAM against a multitude of data sets. We provide further tests here by exploiting the results from a recent magnetic connectivity study of Solar Orbiter data carried out by J.-B. Dakeyo et al. (2024). In their work, proton measurements from the Proton-Alpha Sensor (PAS; C. J. Owen et al. 2020; P. Louarn et al. 2021) were used in combination with interplanetary magnetic field observations from the magnetometer MAG (T. S. Horbury et al. 2020) and a magnetostatic model of the corona and interplanetary magnetic field to connect wind properties to their coronal sources.

The free parameters of the model, related to the geometry of the solar magnetic field (Equation (12)) and the heating terms described in Section 2, can be further constrained on the basis of this Solar Orbiter study. This is illustrated in Figure 2, which presents the relation between the photospheric field strength ($B_0$) and the expansion factor ($f_{max}$) of magnetic field lines connected to each solar wind parcel measured by PAS for magnetic field lines rooted at the photosphere at high (Figure 2(a)) and low (Figure 2(b)) latitudes. The wind velocity ($U$) measured by PAS is also indicated by the color scale. These measurements cover from 2020 August 1 to 2022 Mach 17, which corresponds to the rising phase of the solar minima.

The expansion factor computed by J.-B. Dakeyo et al. (2024) and appearing in Figure 2 is calculated with a potential field source surface (PFSS) extrapolation of solar magnetograms (M. D. Altschuler & G. Newkirk 1969). In PFSS models, the
magnetic field is assumed radial at the photosphere and at a theoretical spherical surface of heliocentric radius $r_{ss}$ called the source surface. Typically, this radius is taken to be $r_{ss} = 2.5 R_{\odot}$ to mimic the radial rays of solar eclipses (Y. M. Wang & N. R. Sheeley 1991). Since the field remains radial beyond that surface, the expansion factor $f_{ss}$ of the magnetic field between the photosphere and the source surface is exactly equal to $f_{max}$. We will therefore use $f_{ss}$ for the rest of this study to indicate the asymptotic value of the expansion factor.

Figure 2(a) reveals a positive correlation between the magnetic footpoint strength $B_0$ and the expansion of the field line $f_{ss}$, especially for lower $f_{ss}$, for which the correlation is clearer with a smaller scattering. This dependence of $f_{ss}$ on $B_0$ means that the total energy flux term defined as $F_{B_0}^{\text{Tot}} = F_{B_0}^{\text{cor}} + F_{B_0}^{\text{TR}}$ and the coronal heating scale height are not independent parameters. $B_0$ is less and less correlated to $f_{ss}$ as $f_{ss}$ increases, and Figure 2(b) reveals that solar winds sharing similar $f_{ss}$ can originate on magnetic flux tubes with very different $B_0$.

The correlation between the footpoint field strength and the expansion factor shown in Figure 2 is expected from Equation (12) since $B(r_{ss})$ undergoes little variation in comparison with $B_0$. The reason $B_0(f_{ss})$ in Figure 2(a) levels off much more sharply than in Figure 2(b) is also expected since source regions located at high latitudes (HLs) are not associated with strong magnetic fields, unlike the active region belts.

Most of the wind parcels originating from HL (polar) coronal holes, shown in Figure 2(a), follow the well-known anticorrelation between wind velocity and expansion factor (Y. M. Wang & N. R. Sheeley 1990) as noted in J.-B. Dakeyo et al. (2024). This anticorrelation is, however, not as strongly verified for wind parcels originating from LLs (Figure 2(b)). This is supported by past magnetic connectivity studies (Y. M. Wang 2016), which compared the coronal magnetic field with solar wind properties and found very little correlation between $B_0$ and the solar wind speed measured near Earth.

The upper part of the $B_0(f_{ss})$ distribution of Figure 2(b) is characterized by the presence of fast winds originating from the LL active region belt. These solar wind streams could undergo a different acceleration regime. J.-B. Dakeyo et al. (2024) showed that these winds become supersonic inside the region of magnetic field superradial expansion, well below the source surface. They relate these unusually high speeds as the result of an early supersonic flow that is further accelerated by a de Laval nozzle effect (H. S. Seifert et al. 1947) in the strongly diverging field (very high $f_{ss}$) in the low corona. We will show that this early supersonic flow is driven by a combination of strong heating at the coronal base (due to high $B_0$) and significant flux tube expansion (de Laval nozzle effect).

PFSS reconstructions should be used with caution, especially (1) during elevated activity levels, (2) in the direct vicinity of active regions, or (3) near streamer structures and their associated current sheets. For this reason, we consider statistical trends from Figure 2 instead of individual points that may be prone to significant backmapping errors. An analysis of the uncertainties in backmapping was carried out in J.-B. Dakeyo et al. (2024), and we refer the reader to their Appendix A for further information. During the preparation of their study, they also compared the estimated source locations directly with ultraviolet imaging of corona holes (private communication) and confirmed that the wind population associated with moderate wind speeds of 400–600 km s⁻¹ and elevated expansion factors are rooted in LL narrow (less

than 5° in width) coronal holes with strong magnetic field strengths. This magnetic connectivity was found to be long-lasting (1–2 days) due to the spread of expanding field lines. One should note that LL coronal holes produce typically slow solar wind, especially during solar maximum periods.

The high-U, high-$f_{ss}$ values are related to special conditions at the source of the solar wind. The present paper investigates whether very strong base heating related to high surface magnetic fields (input energy flux) could be the source of extra acceleration. This was of course not necessarily expected since strong base heating induces higher mass flux evaporated into the tube. What we find is that for the strongest base heating (highest values of $B_0$), the enhanced mass flux in the tube contributes to saturate the maximum speed in the range of 500–550 km s⁻¹, which roughly corresponds to the fastest winds measured by PSP and Solar Orbiter close to the Sun.

The correlations shown in Figure 2 can be used to reduce the number of free parameters of our model by providing relations between base heating $B_0$ and field geometry $f_{ss}$. To this end, we fitted two different parts of the distribution to capture the variety of observed solar wind conditions associated with $B_0–f_{ss}$ relations. The two fits follow the lowest and highest parts of the $B_0(f_{ss})$ distribution, corresponding respectively to an HL simulation series and an LL simulation series.

The parametric equations describing these fits are

$$B_0^{\text{HL}} = 70 \times \tanh\left(\frac{f_{ss}}{20}\right) \text{ G}, \quad (21)$$

$$B_0^{\text{LL}} = 3.2 \times f_{ss}^{1.1} \text{ G}. \quad (22)$$

These relations, combined with our heating prescription described in Section 2.3.2, make $f_{ss}$ the free parameter that defines not only the field geometry and its effect on the divergence terms in the conservation equations but also the heating scale height $H_f$ and the energy flux $F_{B_0}$ through $B_0$.

Each simulation series to follow in the next section is made up of nine simulations corresponding to $B_0–f_{ss}$ pairs indicated by the black crosses in Figure 2. The nine simulations correspond to nine values of the expansion factors $f_{ss} = 2, 3, 5, 8, 10, 13, 20, 25, 30$. The three simulation series share the same three first $B_0–f_{max}$ values situated on the superposed portions of the black curves.

## 4. Wind Simulations for HL Coronal Holes

We begin by modeling the solar wind emerging from a polar coronal hole by considering the HL simulation series defined in the previous section. The lower boundary conditions and all input parameters are shown in Table 1. The simulation results are compared to remote sensing and in situ observations in Figure 3.

The profiles extend from the highly collisional photosphere at $1 \times 10^{-3} R_{\odot}$ (~700 km) to $20 R_{\odot}$, in the collisionless solar wind. The transition region visible as the sharp transition in plasma moments located near $3 \times 10^{-3} R_{\odot}$ (~2000 km) agrees closely with the temperature and density jumps given by the non-local thermodynamic equilibrium (LTE) radiative model of E. H. Avrett & R. Loeser (2008), shown as black profiles in the different panels.

Increasing $f_{ss}$ reduces $H_f$, the main effect of which is to concentrate coronal heating below the sonic point, in the still-collisional part of the corona. Heating in this region increases the base temperature, the conductive heat flux toward the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Simulation results for the HL series as a function of height above the photosphere. The color map indicates the $f_{ss}$ value associated with each simulation. Panel (a): proton density compared with observations; blue plus signs correspond to SOHO-SUMER off-limb measurements from E. Landi (2008), and blue diamonds correspond to UVCS measurements from S. R. Cranmer (2020). Panel (b): proton temperature compared with observations; gray patches correspond to UVCS measurements from S. R. Cranmer (2020), and black vertical lines with caps correspond to measurements of PSP (J.-B. Dakeyo et al. 2022). Panel (c): electron temperature compared with observations; blue plus signs correspond to SOHO-SUMER off-limb measurements (E. Landi 2008), gray patches correspond to UVCS measurements from S. R. Cranmer (2020), and black vertical lines with caps correspond to measurements of PSP (J.-B. Dakeyo et al. 2022). Panel (d): proton velocity compared with observations; gray patches correspond to UVCS measurements from S. R. Cranmer (2020), and blue dots correspond to blob velocity in the Heliospheric Current Sheet (HCS) measured by SOHO/LASCO and STEREO/COR (L. Abbo et al. 2016). Panel (e): oxygen CSR $n_{O^{7+}}/n_{O^{6+}}$ calculated using the method described in Appendix A. For panels (a)–(d): black profiles in the chromosphere correspond to the non-LTE radiative semiempirical model from E. H. Avrett & R. Loeser (2008).</figcaption></figure>

<table><caption>**Table 1**
Input Parameters for the Different Simulations within the HL Series Calculated Using Equations (15), (16), (17), and (21)</caption><thead><tr><th>Input Parameters</th><th>HL<sub>2</sub></th><th>HL<sub>3</sub></th><th>HL<sub>5</sub></th><th>HL<sub>8</sub></th><th>HL<sub>10</sub></th><th>HL<sub>13</sub></th><th>HL<sub>20</sub></th><th>HL<sub>25</sub></th><th>HL<sub>30</sub></th></tr></thead><tbody><tr><td>f<sub>max</sub></td><td>2</td><td>3</td><td>5</td><td>8</td><td>10</td><td>13</td><td>20</td><td>25</td><td>30</td></tr><tr><td>B<sub>0</sub> (G)</td><td>2.9</td><td>4.3</td><td>7</td><td>10.9</td><td>13.2</td><td>16.4</td><td>21.8</td><td>24.3</td><td>25.9</td></tr><tr><td>H<sub>f</sub><sup>TR</sup> (10<sup>-2</sup> R<sub>☉</sub>)</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>H<sub>f</sub><sup>cor</sup> (R<sub>☉</sub>)</td><td>2.6</td><td>1.9</td><td>1.3</td><td>0.9</td><td>0.77</td><td>0.63</td><td>0.46</td><td>0.39</td><td>0.33</td></tr><tr><td>F<sub>B0</sub><sup>TR</sup> (10<sup>5</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</td><td>1.25</td><td>1.85</td><td>3</td><td>4.69</td><td>5.7</td><td>7.1</td><td>9.4</td><td>10.5</td><td>11.1</td></tr><tr><td>F<sub>B0</sub><sup>cor</sup> (10<sup>5</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</td><td>2.5</td><td>3.7</td><td>6</td><td>9.38</td><td>11.4</td><td>14.2</td><td>18.8</td><td>21</td><td>22.2</td></tr><tr><td>r<sub>p</sub><sup>sw</sup> (R<sub>☉</sub>)</td><td>4.80</td><td>3.6</td><td>2.5</td><td>1.68</td><td>1.5</td><td>1.2</td><td>0.8</td><td>0.7</td><td>0.6</td></tr><tr><td>r<sub>e</sub><sup>sw</sup> (R<sub>☉</sub>)</td><td>6.8</td><td>4.3</td><td>2.6</td><td>1.63</td><td>1.3</td><td>1.0</td><td>0.6</td><td>0.5</td><td>0.4</td></tr></tbody></table>

Note. $r_e^{sw}$ and $r_p^{sw}$ have been chosen to fit PSP data at 20 R<sub>☉</sub>.

transition region, thereby raising its pressure and chromospheric evaporation. This increases the mass flux and therefore the density in the low corona. These effects are visible in

Figure 3 panels (a)–(d), where *n*, *T*<sub>p</sub>, *T*<sub>e</sub>, and *U* all increase in the lower corona with increasing values of *f*<sub>ss</sub>. The increased coronal pressure pushes the transition region sunward, and the
best correspondence between our model and the E. H. Avrett & R. Loeser (2008) profiles is obtained for the highest values of $f_{ss}$.

Concentrating energy deposition in the lower atmosphere decreases the energy available for solar wind acceleration in the upper corona. The highest terminal wind speeds are obtained for the smallest expansion factor values, $f_{ss} = 2, 3$, and $5$, and the associated curves agree with the temperatures and speeds derived from spectroscopic observations made with SOHO/UVCS of polar coronal holes (S. R. Cranmer 2020).

Increasing $f_{ss}$ results in a progressive decrease in wind speed and pushes the bulk of the wind acceleration higher up in the solar atmosphere. The lowest terminal wind speeds are for $f_{ss} = 20, 25$, and $30$, and their acceleration closely follows that of transient structures released in the heliospheric current sheet observed as “blobs” by coronagraphs on board SOHO and the Solar TErrestrial RElations Observatory (STEREO; L. Abbo et al. 2016; shown as blue dots in Figure 3(d)). We do not expect wind velocities to be lower than those of blobs considered “passive” structures entrained in the slow wind (N. R. Sheeley et al. 1997; V. Réville et al. 2020).

Collisions in the lower corona maintain $T_e \approx T_p$ until the medium becomes collisionless and $T_e$ and $T_p$ begin to differ near $3 \times 10^{-1} R_⊙$ (~2000 km) for the fastest winds but at higher altitudes for the slower winds. The simulations shown in Figure 3 also retrieve the well-known correlation between $T_p$ and $U$ and the anticorrelation between $T_e$ and $U$ in the escaping solar wind (J.-B. Dakeyo et al. 2022). Comparing $T_e$ and $T_p$ at the outer boundary with particle measurements by PSP (black vertical lines with caps shown in panels (b) and (c)) provided by J.-B. Dakeyo et al. (2022), we see that the simulations slightly overestimate the range of temperatures.

We also note from Figures 3(b) and (d) that for the fast wind to have speeds of at least ~600 km s⁻¹, simulated proton temperatures must exceed those inferred in coronal holes from spectroscopy as well as in the nascent solar wind measured in situ by PSP. This points to the well-known result that thermal pressure gradients are not sufficient to produce the highest wind speeds measured in situ, and additional sources of momentum must be considered (S. R. Cranmer et al. 2007; V. Pierrard & M. Peters de Bonhome 2024). This missing momentum source will be addressed in Section 7.

## 5. Source Temperatures at HLs

Figure 4 presents in situ measurements by the ACE of the oxygen CSRs $n_{O^{7+}}/n_{O^{6+}}$ for the period 1998 February 6 to 21 August 2011. The distribution exhibits the well-known anticorrelation between source temperature ($n_{O^{7+}}/n_{O^{6+}}$) and wind speed introduced earlier, which implies a seemingly hotter source for slower solar winds. The slow wind also tends to have a broader range of possible source temperatures if we consider the plot is in logarithmic scale. As already discussed, this is consistent with the multicomponent nature of the slow wind indicating different possible conditions at the source and/or different formation mechanisms. The highest CSRs come from slow winds originating in the vicinity of coronal streamers where enhanced charge states could result from coronal loop material ejected into the solar wind by reconnection with open field lines, which may provide an explanation for their high value (N. R. Sheeley et al. 1997; E. Sanchez-Diaz et al. 2017; V. Réville et al. 2020, 2022; A. P. Rouillard et al. 2021). In the

framework of the quasi-static theory, the modeling of these winds is not considered.

The periods identified in R. D'Amicis & R. Bruno (2015) of slow winds sharing similar Alfvénicity and compressibility to the fast wind are plotted as red dots in Figure 4. The fast wind periods also identified by these authors during this solar maximum period are shown as green dots. These fast and slow winds could be considered a continuum determined by the expansion rate of the coronal magnetic field. For these two wind components, we can still observe an anticorrelation of the CSRs with wind speed that can be investigated using our HL ISAM simulations tested in the previous section against a broad range of observations.

CSRs $n_{O^{7+}}/n_{O^{6+}}$ were calculated from the ISAM simulations in two different ways in the present paper: (1) by using the electron moments modeled by ISAM as input to a solution of the system of continuity equations in Y.-K. Ko et al. (1997; detailed in Appendices A and C) and (2) by solving for the full transport and ionization of heavy ions treated as an additional population in ISAM. The two approaches gave similar results that are presented in Appendix D. A direct solution of the heavy ion transport accounts for differential flows that are found to induce slightly higher ionization levels than without, as discussed in Appendix D.

The CSRs shown in Figure 3(e) and also reported in Figure 4 were therefore computed using the simpler and less computationally intensive method of Y.-K. Ko et al. (1997). As seen in previous studies, the $n_{O^{7+}}/n_{O^{6+}}$ evolves rapidly near the transition region and settles below $1 R_⊙$, remaining invariant further out. Increasing $B_0$ and its associated input energy flux, $n_{O^{7+}}/n_{O^{6+}}$ increases in response to the induced rise in $n$ and $T_e$. The wind velocity $U$ increases between consecutive simulations, but it is not sufficient to decrease $n_{O^{7+}}/n_{O^{6+}}$ via nonequilibrium effects.

These simulations illustrate how a displacement of energy deposition toward the base of the corona decreases the terminal wind speed and increases the coronal base temperature and density, thereby increasing the charge state of heavy ions. This was already shown by Y.-M. Wang et al. (2009) and is further verified here through a more comprehensive and systematic comparison with a broad range of data sets. We note from Figure 4 that the simulated CSRs match the lower values of the distributions but are generally lower than those measured in situ. This discrepancy will be revisited in Sections 6 and 7.

## 6. Simulating Wind Sources from the Active Region Belt

A significantly enhanced input energy flux is expected for wind sources located within the active region belt where the strongest photospheric magnetic fields are encountered. Such conditions correspond to the $B_0(f_{ss})$ values along the upper curve (LL) shown in Figure 2(b). Along this branch, $B_0$ values at $f_{max} = 25, 30$ are over an order of magnitude higher than for the HL runs as reported in Table 2.

Increasing $B_0$, and therefore the related input energy flux, forces the wind to transition into a new expansion regime. For the subset of simulations related to the highest values of $B_0(f_{ss})$, the wind speeds in the low corona are much higher than for similar $f_{ss}$ values in the HL simulations (Figure 3). Supersonic speeds develop already in the region of strong superradial field expansion in the low corona. R. A. Kopp & T. E. Holzer (1976) and J.-B. Dakeyo et al. (2024) studied the effect of such flux tube geometry on the wind velocity profiles using polytropic
<table><caption>**Table 2**

Input Parameters for the Different Simulations within the LL Series Calculated Using Equations (15), (16), (17), and (22)</caption><thead><tr><th>Input Parameters</th><th>LL<sub>2</sub></th><th>LL<sub>3</sub></th><th>LL<sub>5</sub></th><th>LL<sub>8</sub></th><th>LL<sub>10</sub></th><th>LL<sub>13</sub></th><th>LL<sub>20</sub></th><th>LL<sub>25</sub></th><th>LL<sub>30</sub></th></tr></thead><tbody><tr><td>f<sub>max</sub></td><td>2</td><td>3</td><td>5</td><td>8</td><td>10</td><td>13</td><td>20</td><td>25</td><td>30</td></tr><tr><td>B<sub>0</sub> (G)</td><td>2.9</td><td>4.3</td><td>7</td><td>13.8</td><td>17.6</td><td>23.5</td><td>37.8</td><td>48.3</td><td>59.0</td></tr><tr><td>H<sub>f</sub><sup>TR</sup> (10<sup>-2</sup> R<sub>☉</sub>)</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>H<sub>f</sub><sup>cor</sup> (R<sub>☉</sub>)</td><td>2.6</td><td>1.9</td><td>1.3</td><td>0.9</td><td>0.77</td><td>0.63</td><td>0.46</td><td>0.39</td><td>0.33</td></tr><tr><td>F<sub>B<sub>0</sub></sub><sup>TR</sup> (10<sup>5</sup> erg cm<sup>-2</sup>s<sup>-1</sup>)</td><td>1.25</td><td>1.85</td><td>3</td><td>5.54</td><td>7.6</td><td>10.1</td><td>16.25</td><td>20.52</td><td>25.35</td></tr><tr><td>F<sub>B<sub>0</sub></sub><sup>cor</sup> (10<sup>5</sup> erg cm<sup>-2</sup>s<sup>-1</sup>)</td><td>2.5</td><td>3.7</td><td>6</td><td>11.8</td><td>15.2</td><td>20.2</td><td>32.5</td><td>41.5</td><td>50.7</td></tr><tr><td>r<sub>p</sub><sup>sw</sup> (R<sub>☉</sub>)</td><td>4.80</td><td>3.6</td><td>2.5</td><td>1.68</td><td>1.5</td><td>1.2</td><td>0.8</td><td>0.7</td><td>0.6</td></tr><tr><td>r<sub>e</sub><sup>sw</sup> (R<sub>☉</sub>)</td><td>6.8</td><td>4.3</td><td>2.6</td><td>1.63</td><td>1.3</td><td>1.0</td><td>0.6</td><td>0.5</td><td>0.4</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** In situ measurements by ACE of the ion CSRs $n_{O^{7+}}/n_{O^{6+}}$ (panel (b)) and $n_{C^{6+}}/n_{C^{5+}}$ (panel (a)) as a function of the proton wind speed compared with the solutions of Equation (A1). The light gray dots represent the solar minimum corresponding to the period 2006 September 4–2009 September 24. The dark gray dots represent the solar maximum corresponding to the period 1999 March 7–2003 December 22. Black crosses correspond to the solutions of Equation (A1) associated with the HL simulation series. The data corresponding to the time intervals of the fast wind, Alfvènic slow wind, and non-Alfvènic slow wind identified in the study of R. D’Amicis & R. Bruno (2015) are shown as green and red squares and the green stars correspond to a high speed stream measured at solar minimum, respectively.</figcaption></figure>

models and a combination of isothermal and polytropic modeling, respectively. They showed that when the sonic point is brought below or within the region of superradial expansion, the solar wind experiences a significant acceleration.

For such a high regime of $B_0$ and $f_{ss}$, most of the energy is deposited in the low corona, and the wind accelerates only gradually in the upper corona to reach moderately high speeds. Nevertheless, this regime explains the occurrence of the faster wind speeds (>450 km s⁻¹) observed for high $f_{max}$ in Figure 2 and discussed in detail in J.-B. Dakeyo et al. (2024).

The strong base heating induces some significant chromospheric evaporation into the corona. However, the increased plasma acceleration and the high $f_{ss}$ induce a significant expansion of the plasma acting to decrease the density and limiting the rise in temperature in the low corona. Low densities and increased temperatures lead to a drop in the collision frequencies allowing a strong divergence of electron and ion temperatures. The electrons rapidly decouple from the ions, which reduces the energy transmitted by collisions to the

ions, resulting in a decrease in the ionic temperature between $3 \times 10^{-3}$ (~2000 km) and $4 \times 10^{-2}$ (~27,000 km).

According to Equation (A1), the increase in wind speed and the decrease in density limit the rise in the ionization level. This can be seen in the profiles corresponding to the second expansion regime ($f_{ss} = 25, 30$) in Figure 5, where the high velocities (up to ~200 km s⁻¹) and low densities in the collisional region of the corona result in oxygen CSRs no longer following the increasing trend, falling below the charge ratios associated with the $f_{ss} = 13, 20$ profiles.

The carbon and oxygen CSRs from our LL simulation series were calculated from Equation (A1) (as for the HL simulations series) and are compared with the in situ ACE measurements in Figure 6.

The two $B_0(f_{ss})$ relations (Figure 2) considered for the HL (Figure 3) and LL (Figure 5) simulations provide the lower and upper limits of the base heating conditions and possible CSRs. These simulations provide an interpretation for the anticorrelation between wind velocity and CSRs. However, as with other attempts to model CSRs with thermal models (S. R. Cranmer et al. 2007; R. Oran et al. 2015), the CSR values are below those measured in situ. Even the LL series winds with the highest electron temperatures at the base provide CSR values that are below the center of the distribution.

Past modeling studies obtained even lower $n_{O^{7+}}/n_{O^{6+}}$ values below observations by a factor of 10, demonstrating further that thermal models are unable to retrieve the observed charge states (S. R. Cranmer et al. 2007; S. R. Cranmer 2014; R. Oran et al. 2015). Additional processes should be considered, and we now consider the combined effects of nonthermal particles on the solar wind properties and CSRs.

## 7. Wind Solutions Including Suprathermal Electrons

Nonthermal particle populations are ubiquitous in the solar wind as detected in numerous measurements of particle velocity distribution functions (M. D. Montgomery et al. 1968; W. C. Feldman et al. 1975; W. G. Pilipp et al. 1987; M. Maksimovic et al. 1997). Suprathermal electrons tend to follow a kappa or generalized Lorentzian velocity distribution in the solar wind (V. M. Vasyliunas 1968). The high-energy tail then follows a power law in particle velocity,

$$f_{\kappa}(v_e, \kappa) = n_e \left(\frac{m_e}{2\pi k_B T_e}\right)^{\frac{3}{2}} \frac{\Gamma(\kappa + 1)}{\kappa^{\frac{3}{2}}\Gamma(\kappa - 0.5)} \left[1 + \frac{v_e^2}{\kappa\theta_e^2}\right], \quad (23)$$

where $v_e$ is the microscopic velocity of electrons, $\Gamma$ is the gamma function, $\theta_e = \sqrt{\frac{2k_B T_e}{m_e}}$ is the thermal velocity of electrons, and $\kappa$ is the parameter characterizing the amplitude of the suprathermal
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Simulation results for wind sources rooted at LLs as a function of height above the photosphere. Comparisons with observations and models are the same as in Figure 3.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** In situ measurements by ACE of the ion CSRs $n_{O^{7+}}/n_{O^{6+}}$ (panel (b)) and $n_{C^{6+}}/n_{C^{5+}}$ (panel (a)) as a function of the proton wind speed compared with the solutions of Equation (A1). The light gray dots represent the solar minimum corresponding to the period 2006 September 4–2009 September 24. The dark gray dots represent the solar maximum corresponding to the period 1999 April 7–2003 December 22. Black crosses correspond to the solutions of Equation (A1) associated with the LL simulation series.</figcaption></figure>

population with $f_κ(v_e, κ = ∞) = f_{Maxwellian}(v_e)$. It is important that $κ > 3/2$ so that the integral of the distribution function converges and the moments are well defined.

Past modeling studies suggest the existence of suprathermal electrons already in the corona based on considerations of coronal turbulence (D. A. Roberts & J. A. Miller 1998), heating processes involving obliquely propagating finite-amplitude electromagnetic waves (A. F. Viñas et al. 2000), nanoflares (H. Che & M. L. Goldstein 2014), and resonant interactions with whistler waves (C. Vocks & G. Mann 2003).

Suprathermal electrons become collisionless in the mid-corona ($\sim 2 \times 10^{-1} R_⊙$), and a subset of electrons have sufficient energy to overcome the sunward-directed electrical and gravitational forces and escape into the solar wind. These escaping electrons contribute to the acceleration of the solar wind via the ambipolar electric field (V. Pierrard & M. Peters de Bonhome 2024). Our thermal model based on bi-Maxwellian particle distributions could not produce terminal wind speeds of $\sim 700 \text{ km s}^{-1}$ from the measured coronal temperatures (see Sections 4 and 5). Additional effects must occur, and the ambipolar electric field generated by suprathermal electrons could be the missing force to accelerate the solar wind. In the low corona, high-energy electrons should also increase the CSRs of heavy ions (Y.-K. Ko et al. 1996;
S. R. Cranmer 2014; R. Oran et al. 2015). This section jointly investigates the effects of nonthermal electrons on the solar wind energy budget and the ionization levels of heavy ions.

Exospheric models of the forming solar wind are the simplest approach to model the evolution and effects of the electron distribution function. They assume the absence of Coulomb interactions between particles after a certain altitude, referred as the exobase (V. Pierrard et al. 2021). Each electron is therefore uniquely subject to external forces, including the electric, gravitational, and Lorentz forces. These simplifications combined with a stationarity condition and a predefined distribution function at the exobase permit an analytical solution of the Vlasov equation to get the radial profile of the velocity distribution function. Inclusion of the ambipolar electric force in ISAM was carried out by using the exospheric model of V. Pierrard et al. (2023), and the modeled ambipolar electric field was included as an extra force in the momentum (Equation (2)).

We included the ambipolar electric field in three characteristic simulations taken from our two sets of simulations: a fast ($f_{ss} = 5$) and slow ($f_{ss} = 20$) wind simulation taken from the HL set (representative of sources inside the polar coronal hole) and a slow wind ($f_{ss} = 25$) simulation taken from the LL set (representative of sources located near active regions). We will refer to them respectively as HL$_5$, HL$_{20}$, and LL$_{25}$. We used these three simulations to first compute the location of the exobase identified as the inflection point in the wind acceleration profile and found ~0.7 $R_⊙$, 2 $R_⊙$, and 2.3 $R_⊙$, respectively, for HL$_5$, HL$_{20}$, and LL$_{25}$. The ISAM profiles for $n$, $U$, $T_p$, and $T_e$ taken at the exobase were then given as input for the exospheric model. The latter was run for different values of $\kappa$ providing radial profiles of the ambipolar electric field. To avoid the discontinuity caused by the abrupt inclusion of the electric field at the exobase in ISAM, we added a smoothly decreasing function to the electric field in the sunward direction.

$\kappa = 6$ provided the best fit between modeled and observed plasma moments as shown in Figure 7. As expected, the additional momentum source provided by the applied ambipolar electric field induces wind speeds in excess of ~600 km s⁻¹ for observed proton temperatures of ~2 × 10⁶ K. $\kappa = 6$ is consistent with the values 2 < $\kappa$ < 6 inferred by G. Gloeckler et al. (1992) and M. Maksimovic et al. (1997) from electrons measured in situ.

We then calculate the CSRs for these three winds by taking into account the enhancement of ionization associated with a kappa distribution with $\kappa = 6$. To this end, we follow an approach similar to that of Y.-K. Ko et al. (1996). To use the ionization/recombination rates precalculated for Maxwellian electron distributions, the kappa distribution is first approximated by a sum of Maxwellians, each characterized by a different density and temperature. The resulting ionization and recombination rates are then the sum of the ionization and recombination rates determined from the Maxwellians. The method is further described in Appendix C.

The new CSR values for our HL$_5$, HL$_{20}$, and LL$_{25}$ simulations are reported in Figure 8 and compared with CSRs calculated using a Maxwellian eVDF from their associated simulations without the inclusion of an ambipolar electric field.

The addition of a suprathermal population has two effects visible in Figure 8: the ambipolar electric field increases the terminal velocity of the winds, and the high-energy tail of the distribution enhances the ionization and thus the CSRs. For the three simulations considered HL$_5$, HL$_{20}$, and LL$_{25}$, we obtain

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Result of the simulation including the ambipolar electric field obtained with the exospheric model (V. Pierrard et al. 2023) and applied in ISAM for $\kappa = 6$. Solid, dashed, and dotted red lines correspond respectively to profiles from HL$_5$, HL$_{20}$, and LL$_{25}$ simulations. The rest of the notation is the same as in Figure 3.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Comparison of in situ measurements by ACE of the ion CSRs at solar minima with the solutions of Equation (A1). The light gray dots represent the solar minimum corresponding to the period 2006 September 4–2009 September 24. The dark gray dots represent the solar maximum corresponding to the period 1999 April 7–2003 December 22. Black squares correspond to the solutions of Equation (A1) from the HL$_5$, HL$_{20}$, and LL$_{25}$ simulations assuming a Maxwellian distribution function. Red diamonds correspond to the solutions of Equation (A1) assuming a kappa distribution function. The black lines linking the squares and diamonds link the simulations including suprathermal effects to their associated simulation without suprathermal effects. Panel (a): $n_{C^{6+}}/n_{C^{5+}}$ CSR distribution. Panel (b): $n_{O^{7+}}/n_{O^{6+}}$ CSR distribution.</figcaption></figure>

an enhancement of CSRs by a factor of ~3–10 for $C^{6+}/C^{5+}$ and a factor of ~10–30 for $O^{7+}/O^{6+}$. This results in a better fit of the oxygen CSR with solar maximum measurements, while
the carbon CSR fits better with solar minimum measurements. Adding the effect of suprathermals makes it possible to recover CSR values that are consistent with the measurements while retaining the characteristic anticorrelation between wind speed and CSRs.

## 8. Discussion and Outlook

In this study, a plasma transport model extending from the chromosphere to the solar wind was used to model the thermodynamic coupling between the different layers of the solar atmosphere in the estimated source region of the solar wind measured by Solar Orbiter. Neutral particles were included to model ionization processes and equilibria near the transition region.

The spatial distribution of coronal and solar wind heating was based on heating functions (G. L. Withbroe 1988; R. F. Pinto & A. P. Rouillard 2017) that have been shown to approximately match the energy deposition by a turbulent cascade of Alfvén waves. In this respect, the spatial distribution of the energy deposition was made dependent on the large-scale magnetic field properties. The expansion rate of the magnetic field lines, along which the solar wind expands, controls the gradient in Alfvén speed, which in turn controls the reflection of Alfvén waves and the distribution of dissipation with height. The footpoint field strength was taken to be proportional to the Alfvén wave energy flux injected at the base of the tube. This corresponds to perpendicular velocity perturbations at the inner boundary. The expansion factor and the surface field strength, taken from the connectivity study of J.-B. Dakeyo et al. (2024), were found to be correlated over ranges of low to medium expansion factor values ($f_{ss} \lesssim 20$).

With these constraints, our first goal was to verify that the model can reproduce recent in situ measurements of solar wind properties by PSP and SolO in conjunction with a wide range of coronal source properties given by remote sensing observations. A parametric study was carried out, and we present a set of simulations that well reproduced wind properties originating from HLs and LLs. This work confirmed previous results showing that in models coupling the different layers of the low solar atmosphere, the height of the energy deposition strongly influences the wind properties through a number of mechanisms such as chromospheric evaporation and thermal expansion above and below the sonic point. This work retrieves the well-known anticorrelation between wind speed and expansion factors, especially for solar wind originating from HLs. For the solar wind originating from LLs, i.e., footpoints in the vicinity of a strong photospheric magnetic field, which is also associated with large coronal heating, the solar wind undergoes a significant expansion, pushing the sonic point below the region of strong magnetic field expansion. We see that this regime leads to a strong acceleration of the solar wind in the low corona as already discussed by J.-B. Dakeyo et al. (2024) using their isopoly wind model. This expansion regime could also explain the elevated plasma speeds (~100–200 km s⁻¹) inferred near wind sources from remote sensing observations (A. Bemporad 2017) or the persistent outflows observed by coronal spectroscopy on the edge of active regions with velocities exceeding 100 km s⁻¹ as inferred by L. K. Harra et al. (2008). Although the latter outflows have been interpreted as source locations of the solar wind, and our simulations provide some theoretical support, it should be

noted that they could also result from plasma transport along large-scale magnetic loops.

The aim of our study was to gain more insight into the mechanisms that control the charge state of heavy ions or, equivalently, the electron temperature at the sources of the different types of solar winds. For the HL source regions, when energy deposition is concentrated at lower/upper heights in response to strong/weak magnetic flux tube expansions, the charge state increases/decreases. This naturally produces an anticorrelation between ion charge state and solar wind speed and remains clear if the surface field strength is not too strong. This interpretation was already proposed in Y.-M. Wang et al. (2009) and explained in terms of Alfvén wave dissipation by S. R. Cranmer (2020). We show in addition that the electron temperature at the LL sources increases significantly due to the significant heating that occurs in these regions. However, the increase in charge state is limited by the significant plasma expansion.

Overall, the simulated ion charge states are slightly lower than those measured in the solar wind (Figures 4 and 6). Since the coronal source temperatures can also be inferred from spectroscopic observations and are in agreement with the simulated ones (Figures 3 and 5), additional effects not included in our model must lead to higher charge states in the low corona. Our simulations also show that, at the observed coronal temperatures, a purely thermally driven solar wind cannot reach speeds above 600 km s⁻¹, in contrast to observations of the fast solar wind (Figure 3). Both drawbacks could be overcome by inferring the presence of suprathermal electrons in the corona. Using an exospheric model, we show that these nonthermal particles help to increase the ionization of heavy ions through collisions, as well as the acceleration of thermal protons through the development of a polarizing electric field. We find that a nonthermal population represented by a kappa function with a kappa parameter of around 6 leads to good agreement with observations. This value of kappa is slightly lower than the one inferred by S. R. Cranmer (2014).

PSP detected surprisingly little wind near the Sun with speeds in excess of 500 km s⁻¹ (see, e.g., Y. J. Rivera et al. 2024; E. Samara et al. 2024), and additional effects not accounted for in the present study contribute to the acceleration of the fast solar wind to beyond 600 km s⁻¹ in the inner heliosphere (Y. J. Rivera et al. 2024). We note that additional sources of momentum could be imparted to protons by other mechanisms, such as reconnection outflows above transient events such as bright points known to occur at the base of the solar corona (L. Griton et al. 2020; B. Gannouni et al. 2023; C. Hou et al. 2024). These reconnection events have been associated with the occurrence of jets and microstreams in the solar wind (B. Gannouni et al. 2023) as well as potential contributors to additional ionization of heavy ions (C. Hou et al. 2024). A future study could attempt to run the ISAM model along time-evolving magnetic field lines (i.e., after reconnection events) as ISAM is able to model noninertial effects. This is done regularly to model ionospheric convection with the ionospheric version of the model employed here (A. Marchaudon & P.-L. Blelly 2015). Another source of proton acceleration could also come from momentum exchange with alpha particles, which are known to be differentially heated in the source of the fast solar wind (D. Stansby et al. 2019). As alpha particles are an important species, they also modify the electron/temperature scale heights and possibly the
ionization rates of heavy ions along the flux tubes. These aspects are the subject of an ongoing study that will be submitted for publication in the coming months.

## Acknowledgments

The work at IRAP carried out by P.L., A.P.R., M.L., S.T., J.-B.D, P.-L.B. and V.R. was funded by the ERC SLOW\_ SOURCE project (SLOW\_SOURCE-DLV-819189) and CNES through the APR program. The magnetic connectivity study employed the Magnetic Connectivity Tool developed by the Solar-Terrestrial Relations and Modeling Service (STORMS) and the in situ data was made available by the Centre de Données de la Physique des Plasmas (CDPP). This research was also supported by the International Space Science Institute (ISSI) in Bern through ISSI International Team project \#463 (Exploring The Solar Wind In Regions Closer Than Ever Observed Before) led by L. Harra.

## Appendix A
Charge State Calculation Method

As ions are transported through the solar atmosphere, they undergo ionization and recombination where collisions with electrons are significant. In a stationary hypothesis and in the absence of differential flows between different ionization states within an element, we can combine the conservation equations associated with each ionization state to get a matrix differential equation describing the radial evolution of the ionic fraction $y_i = n_i/(\Sigma_{i=0}^Z n_i)$, where the subscript i indicates the ionization state of the element considered (Y.-K. Ko et al. 1997):

$$u \frac{\partial y_i}{\partial r} = n_e (y_{i-1} C_{i-1}(T_e) - y_i (C_i(T_e) + R_{i-1})(T_e) + y_{i+1} R_i(T_e)), \quad (A1)$$

$$\Sigma_{i=0}^{Z} y_i = 1. \quad (A2)$$

Here $n_e$ is the electron density, $T_e$ is the electron temperature, $u$ is the single velocity for all the ionization states, $C_i$ is the ionization rate (including collisional ionization and autoionization), and $R_i$ is the recombination rate (including radiative recombination and dielectronic recombination). The reaction rates are taken from the CHIANTI 10.0.2 Atomic Database (K. P. Dere et al. 1997; G. Del Zanna et al. 2021). This system of equations is solved using the method described in Appendix B.

This method has the advantage of being fast compared to the Runge-Kutta method but could not be used when the lower boundary condition was situated before the sharp gradients of the transition region. We nevertheless chose this method by setting our boundary condition right after the transition region, where collisions are still important and an ionization equilibrium can still be considered. To test this approach and the hypothesis of negligible differential flows, we compared the results given by this method for the CSR $n_{O^{7+}}/n_{O^{6+}}$ with the solution given by a full resolution of oxygen from the bottom of the chromosphere to the solar wind on a fixed background of protons and electrons using ISAM. This point is discussed in Appendix D. The results of the comparison gave no significant difference on $n_{O^{7+}}/n_{O^{6+}}$; indeed, the full coupling of oxygen has shown no significant differential flows in the region where the charge state evolves. Looking at Equation (A1), we see that when collisions become too rare (attained either when the

density is very low, the temperature is very high, or both), the right-hand side vanishes and the ionic fraction remains invariant along the flow trajectory. This is true only when no differential flows are considered, as a relative velocity between ionization states can alter their relative densities, leading to an evolution of the ionic fraction even in the collisionless regime. The vanishing of the right-hand side can also happen when the velocity of the flow is so high that a plasma parcel cannot reach ionization equilibrium while in the collisional part of the corona. Conversely, a low flow velocity and/or a very high collision rate result in ionization equilibrium. It is important to note that the velocity $u$ in Equation (A1) is the velocity of the ion considered. In the case of this parametric study, we ran our model for protons and electrons, not heavy ions. Hence, we could not provide the velocity profile of minor ions; instead, we provided the velocity of protons. This is equivalent to $u_i = u_p$, where $u_p$ is the velocity of protons. This is a reasonable assumption as the region where charge state evolves is still collisional, forcing $u_i \approx u_p$ by friction. This hypothesis has been used in the majority of the literature solving Equation (A1) using $T_e$, $n_e$, and $u_i$ given by solar wind models. The recombination and ionization rates given by CHIANTI are calculated assuming a Maxwellian distribution function, but the presence of non-Maxwellian features such as a suprathermal electron population would have the effect of enhancing these rates and thus the CSRs. In this case, the method presented here has to be adapted. This point is discussed in Section 7 and Appendix C.

## Appendix B
Resolution Method for the Ionic Fraction Evolution

Equation (A1) is of the following form:

$$\frac{dy_i}{dr} = A_{ij} y_j, \quad (B1)$$

where the matrix $A_{ij}$ is defined as

$$A_{ii} = -\frac{n_e}{u}(C_i + R_{i-1}), \quad (B2)$$

$$A_{i,i-1} = \frac{n_e}{u} R_i, \quad (B3)$$

$$A_{i,i+1} = \frac{n_e}{u} C_{i-1}. \quad (B4)$$

The general solution of this equation is

$$y_i(r) = y_i(r_0) \exp(B_{ji}), \quad (B5)$$

where $r_0$ is the height of the boundary condition $y_i(r_0)$ and $B_{ji} = \int_{r_0}^{r} A_{ji}(r') dr'$. The integral $B_{ji}$ is computed numerically using the trapezoidal rule. We then numerically diagonalize $B_{ji} = V_{ji} D_{ii} (V^{-1})_{ji}$, where $V_{ji}$ and $D_{ii}$ are, respectively, the matrix of eigenvectors and the diagonal matrix of eigenvalues. With $D_{ii}$ being diagonal, we can write $\exp(\int_{r_0}^{r} A_{ji}(r') dr') = V_{ji} \exp(D_{ii}) (V^{-1})_{ji}$, where

$$\exp(D_{ii}) = \text{diag}(e^{\lambda_1}, e^{\lambda_2}, \dots, e^{\lambda_n}) \quad (B6)$$

with $\lambda_1, \lambda_2, \dots, \lambda_n$, the eigenvalues of $A_{ij}$. The solution is then computed by solving Equation (B5) using Equation (B6) and taking $r_0$ to be right after the transition region.
## Appendix C
Ionization/Recombination Rates for a Kappa Distribution

The sums of N Maxwellian distribution functions are each characterized by a different density and temperature:

$$ 
\begin{aligned}
f_M^\Sigma(v_e, T_e^1, T_e^2, \dots, T_e^N, n_e^1, n_e^2, \dots, n_e^N) \\
= \sum_{j=1}^N n_e^j \left( \frac{m_e}{2\pi k_B T_e^j} \right)^{\frac{3}{2}} e^{\left( \frac{m_e v_e^2}{2k_B T_e^j} \right)^{\frac{1}{2}}}.
\end{aligned} 
\quad (C1) $$

We choose the following error function:

$$ \epsilon_{\text{err}} = \sum_{v_e} (|\ln f_{\kappa}| - |\ln f_M^\Sigma|)^2. \quad (C2) $$

Fitting $f_κ$ relies on the minimization of the error function $\epsilon_{\text{err}}$ by finding the optimal combination of parameters ($T_e^1, T_e^2, \dots, T_e^N, n_e^1, n_e^2, \dots, n_e^N$) for each Maxwellian distribution. The best fit of $f_κ$ by $f_M^\Sigma$ is determined using the minimization algorithm of the `fminsearch` function in MATLAB. We proceed to this minimization for each grid point of our profiles. The result of the fit is given in Figure 9 for $N=8$, $κ=6$ at the altitude $r=0.1 R_⊙$. The total ionization and recombination rates are then given by

$$ C = \frac{1}{n_e} \sum_{j=1}^{N} n_e^j C_j(T_e^j), \quad (C3) $$

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Fit of a kappa function with $κ = 6$ by eight partial Maxwellians. The black curve corresponds to the kappa function. The red curve corresponds to the sum of the eight partial Maxwellians.</figcaption></figure>

$$ R = \frac{1}{n_e} \sum_{j=1}^{N} n_e^j R_j(T_e^j). \quad (C4) $$

The ionic fraction is then calculated by solving Equation (A1) using the new ionization and recombination rates. We tried different $κ$ values and their associated increase in CSR $n_{O^{7+}}/n_{O^{6+}}$, and we have determined that $κ = 6$ is the most adapted choice, which is consistent with the profiles obtained when the ambipolar electric field induced by a kappa eVDF with $κ = 6$ is included (Section 7).

## Appendix D
Solving for Oxygen Ions

A self-consistent simulation of oxygen and all its ionization states was performed to demonstrate that the method described in Appendix A was an appropriate means of modeling the ionic fraction evolution. The modeling of oxygen with ISAM alleviates the hypothesis of Equation (A1) by allowing for differential velocities between ions. Because of the numerical difficulty of a self-consistent simulation of e, H, O, and all their ions, we simply simulate a wind for e, H, p and fix it when convergence is reached. On this fixed background, we then solve for oxygen and all its ions, O, O$^+$, O$^{2+}$, ..., O$^{8+}$. We justify this procedure by the small impact that the oxygen fluid has on the dynamics of hydrogen because of its very low density in the solar wind. This allows us to reduce the simulation time considerably. CSRs can be modified by the presence of differential flow speeds for two reasons: (1) a difference in velocity between two ions is associated with a difference in density between these two ions by conservation of mass flow, and (2) a difference in velocity is associated with a difference in time spent in the collisional region where ionization and recombination occur (R. Esser & R. J. Edgar 2001).

Figure 10(b) illustrates that the discrepancy between the modeled ratio of oxygen $n_{O^{7+}}/n_{O^{6+}}$ using ISAM and the solution of Equation (A1) is low. At approximately 1 $R_⊙$, the discrepancy is more pronounced and can be attributed to an increase in the differential flow between O$^{6+}$ and O$^{6+}$ ions shown in Figure 10(a). Despite the generation of differential flows, we conclude that the difference between the two methods is sufficiently low to justify the use of the CSRs modeled by solving Equation (A1).
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Results of a self-consistent modelization of oxygen using ISAM. Panel (a): $u_{O^{7+}} - u_{O^{6+}}$ calculated from our model (dashed line). Panel (b): Comparison of $n_{O^{7+}}/n_{O^{6+}}$ from our model (black line) and from solving Equation (A1) (red crosses).</figcaption></figure>

## ORCID iDs

Paul Lomazzi https://orcid.org/0009-0002-9442-5191

Alexis P. Rouillard https://orcid.org/0000-0003-4039-5767

Nicolas Poirier https://orcid.org/0000-0002-1814-4673

Jean-Baptiste Dakeyo https://orcid.org/0000-0002-1628-0276

Viviane Pierrard https://orcid.org/0000-0001-5014-7682

Victor Réville https://orcid.org/0000-0002-2916-3837

Christian Vocks https://orcid.org/0000-0001-8583-8619

Simon Thomas https://orcid.org/0000-0001-7119-9545

## References

Abbo, L., Ofman, L., Antiochos, S. K., et al. 2016, SSRv, 201, 55

Altschuler, M. D., & Newkirk, G. 1969, SoPh, 9, 131

Athay, R. G. 1986, ApJ, 308, 975

Avrett, E. H., & Loeser, R. 2008, ApJS, 175, 229

Axford, W. I., & McKenzie, J. F. 1992, in Solar Wind Seven: Proc. 3rd COSPAR Colloq., ed. E. Marsch & R. Schwenn (Oxford: Pergamon), 1

Bemporad, A. 2017, ApJ, 846, 86

Bleilly, P. L., Marchaudon, A., Indurain, M., et al. 2019, P&SS, 169, 35

Bleilly, P. L., & Schunk, W. 1993, AnGeo, 11, 443

Boris, J. P., & Book, D. L. 1976, JCoPh, 20, 397

Boris, J. P. L. S., Alex, R. M., Oran, E. S., & Gardner, J. H. 1993, LCPFCT: A Monotone Algorithm for Solving Continuity Equations, Naval Research Laboratory, Laboratory for Computational Physics and Fluid Dynamics, 6410

Byhring, H. S., Cranmer, S. R., Lie-Svendsen, O., Habbal, S. R., & Esser, R. 2011, ApJ, 732, 119

Che, H., & Goldstein, M. L. 2014, ApJ, 795, L38

Cranmer, S. R. 2014, ApJL, 791, L31

Cranmer, S. R. 2020, ApJL, 900, 105

Cranmer, S. R., van Ballegooijen, A. A., & Edgar, R. J. 2007, ApJS, 171, 520

Dakeyo, J.-B., Maksimovic, M., Démoulin, P., Halekas, J., & Stevens, M. L. 2022, ApJ, 940, 130

Dakeyo, J.-B., Rouillard, A. P., Réville, V., et al. 2024, A&A, 691, A77

D'Amicis, R., & Bruno, R. 2015, ApJ, 805, 84

Del Zanna, G., Dere, K. P., Young, P. R., & Landi, E. 2021, ApJ, 909, 38

Dere, K. P., Landi, E., Mason, H. E., Monsignori Fossi, B. C., & Young, P. R. 1997, A&AS, 125, 149

Domingo, V., Fleck, B., & Poland, A. I. 1995, SoPh, 162, 1

Esser, R., & Edgar, R. J. 2000, ApJ, 532, L71

Esser, R., & Edgar, R. J. 2001, ApJ, 563, 1055

Feldman, W. C., Asbridge, J. R., Bame, S. J., Montgomery, M. D., & Gary, S. P. 1975, JGR, 80, 4181

Gannouni, B., Réville, V., & Rouillard, A. P. 2023, ApJ, 958, 110

Geiss, J., Gloeckler, G., & von Steiger, R. 1995, SSRv, 72, 49

Gloeckler, G., Geiss, J., Balsiger, H., et al. 1992, A&AS, 92, 267

Griton, L., Pinto, R. F., Poirier, N., et al. 2020, ApJ, 893, 64

Hansteen, V. H., & Leer, E. 1995, JGR, 100, 21577

Harra, L. K., Sakao, T., Mandrini, C. H., et al. 2008, ApJL, 676, L147

Horbury, T. S., O'Brien, H., Carrasco Blazquez, I., et al. 2020, A&A, 642, A9

Hou, C., Rouillard, A. P., He, J., et al. 2024, ApJL, 968, L28

Issautier, K., Meyer-Vernet, N., Moncuquet, M., & Hoang, S. 1998, JGRA, 103, 1969

Janse, Å. M., Lie-Svendsen, O., & Leer, E. 2007, A&A, 474, 997

King, J. H., & Papitashvili, N. E. 2005, JGRA, 110, A02104

Kinnison, J., Vaughan, R., Hill, P., et al. 2020, in 2020 IEEE Aerospace Conf. (Piscataway, NJ: IEEE), 1

Ko, Y.-K., Fisk, L. A., Geiss, J., Gloeckler, G., & Guhathakurta, M. 1997, SoPh, 171, 345

Ko, Y.-K., Fisk, L. A., Gloeckler, G., & Geiss, J. 1996, GeoRL, 23, 2785

Kohl, J. L., Noci, G., Antonucci, E., et al. 1998, ApJL, 501, L127

Kopp, R. A., & Holzer, T. E. 1976, SoPh, 49, 43

Landi, E. 2008, ApJ, 685, 1270

Lavarra, M. 2022, PhD thesis, Université de Toulouse

Li, X. 1999, SSRv, 87, 253

Lie-Svendsen, ø., Leer, E., & Hansteen, V. H. 2001, JGRA, 106, 8217

Louarn, P., Fedorov, A., Prech, L., et al. 2021, A&A, 656, A36

Maksimovic, M., Pierrard, V., & Riley, P. 1997, GeoRL, 24, 1151

Maneva, Y. G., Viñas, A. F., Moya, P. S., Wicks, R. T., & Poedts, S. 2015, ApJ, 814, 33

Marchaudon, A., & Bleilly, P.-L. 2015, JGRA, 120, 5728

Marchaudon, A., & Bleilly, P. L. 2020, JGRA, 125, e27672

Marsch, E. 2006, LRSP, 3, 1
McKenzie, J. F., Banaszkiewicz, M., & Axford, W. I. 1995, A&A, 303, L45

Montgomery, M. D., Bame, S. J., & Hundhausen, A. J. 1968, JGR, 73, 4999

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

Oran, R., Landi, E., Holst, B. V. D., et al. 2015, ApJ, 806, 55

Owen, C. J., Bruno, R., Livi, S., et al. 2020, A&A, 642, A16

Pierrard, V., Lazar, M., & Maksimovic, M. 2021, in Suprathermal Populations and Their Effects in Space Plasmas: Kappa vs. Maxwellian, ed. M. Lazar & H. Fichtner (Cham: Springer), 15

Pierrard, V., & Peters de Bonhome, M. 2024, arXiv:2401.13308

Pierrard, V., Péters de Bonhome, M., Halekas, J., et al. 2023, Plasma, 6, 518

Pilipp, W. G., Miggenrieder, H., Mühlhäuser, K. H., et al. 1987, JGR, 92, 1103

Pinto, R. F., & Rouillard, A. P. 2017, ApJ, 838, 89

Poirier, C. 2022, PhD thesis, Univ. Paul Sabatier-Toulouse III.

Réville, V., Velli, M., Rouillard, A. P., et al. 2020, ApJL, 895, L20

Rifai Habbal, S., Esser, R., Guhathakurta, M., & Fisher, R. R. 1995, GeoRL, 22, 1465

Rivera, Y. J., Badman, S. T., Stevens, M. L., et al. 2024, Sci, 385, 962

Roberts, D. A., & Miller, J. A. 1998, GeoRL, 25, 607

Rouillard, A., Sheeley, N., Lavraud, B., & Davies, J. 2010, in European Geosciences Union (EGU) General Assembly (Vienna: EGU), 14543

Rouillard, A. P., Savani, N. P., Davies, J. A., et al. 2009, SoPh, 256, 307

Rouillard, A. P., Viall, N., Pierrard, V., et al. 2021, in Space Physics and Aeronomy Vol. 1, Solar Physics and Solar Wind, ed. N. E. Raouafi (New York: Wiley), 1

Réville, V., Fargette, N., Rouillard, A. P., et al. 2022, A&A, 659, A110

Samara, E., Arge, C. N., Pinto, R. F., et al. 2024, ApJ, 971, 83

Sanchez-Diaz, E., Rouillard, A. P., Davies, J. A., et al. 2017, ApJ, 851, 32

Seifert, H. S., Mills, M. M., & Summerfield, M. 1947, AmJPh, 15, 1

Sheeley, N. R., Wang, Y. M., Hawley, S. H., et al. 1997, ApJ, 484, 472

Stansby, D., Perrone, D., Matteini, L., Horbury, T. S., & Salem, C. S. 2019, A&A, 623, L2

Stone, E. C., Frandsen, A. M., Mewaldt, R. A., et al. 1998, SSRv, 86, 1

Vasyliunas, V. M. 1968, JGR, 73, 2839

Verdini, A., Grappin, R., Pinto, R., & Velli, M. 2012, ApJL, 750, L33

Verdini, A., & Velli, M. 2007, ApJ, 662, 669

Viñas, A. F., Wong, H. K., & Klimas, A. J. 2000, ApJ, 528, 509

Vocks, C., & Mann, G. 2003, ApJ, 593, 1134

Wang, Y. M. 2016, ApJ, 833, 121

Wang, Y.-M., & Ko, Y.-K. 2019, ApJ, 880, 146

Wang, Y.-M., Ko, Y.-K., & Grappin, R. 2009, ApJ, 691, 760

Wang, Y. M., & Sheeley, N. R., Jr. 1990, ApJ, 355, 726

Wang, Y. M., & Sheeley, N. R., Jr. 1991, ApJ, 372, L45

Withbroe, G. L. 1988, ApJ, 1, 325

Zhao, L., Zurbuchen, T. H., & Fisk, L. A. 2009, GeoRL, 36, L14104