Charge State Calculation for Global Solar Wind Modeling

The charge state composition of the solar wind carries information about the electron temperature, density, and velocity of plasma in the solar corona that cannot always be measured with remote sensing techniques, due to limitations in instrumental sensitivity and field of view as well as line-of-sight integration issues. However, in situ measurements of the wind charge state distribution only provide the end result of the solar wind evolution from the source region to the freeze-in point. By using 3D global modeling it is possible to follow solar wind plasma parcels of different origin along the path of their journey and study the evolution of their charge states as well as the driving physical processes. For this purpose, we implemented nonequilibrium ionization calculations within the Space Weather Modeling Framework’s solar corona and inner heliosphere modules, to the Alfvén Wave Solar Model (AWSoM). The charge state calculations are carried out parallel to the AWSoM calculations, including all the elements and ions whose ionization-recombination rates are included in the CHIANTI database, namely, from H to Zn. In this work, we describe the implementation of the charge state calculation, and compare simulation results to in situ measurements from the Advanced Composition Explorer and Ulysses spacecraft, and study charge state evolution of plasma parcels along different wind trajectories and wind types.


Introduction
The solar wind is a continuous stream of highly ionized particles released by the Sun into the heliosphere, and it is of critical importance for space weather. Once released from the Sun's surface, the solar wind fills the heliosphere and determines the local plasma properties in interplanetary space, greatly affecting the propagation and geoeffectiveness of traveling disturbances such as coronal mass ejections (CMEs); also, the interaction between solar wind streams of different speeds creates regions of shocked material that can have effects on the near-Earth space environment. Thus, understanding the origin and the acceleration of the solar wind is of critical importance regarding the ongoing efforts in predicting the arrival time and geoeffectiveness of solar storms for the purpose of mitigating their adverse effects.
The charge state composition of the solar wind plasma is of particular importance because ionization and recombination processes are very sensitive to the evolution of the electron density, temperature, and bulk velocity; therefore, charge states carry information not only about the solar wind source but also about the physical processes taking place in the low solar corona (SC), until the ionization and recombination processes stop being effective. This results in different solar wind types coming from different sources in the solar innermost atmosphere carrying different ionization signatures (see, e.g., Neugebauer et al. 2016;Cranmer et al. 2017;Fu et al. 2017;Zhao et al. 2017). Solar wind charge states freeze-in within few solar radii from the solar surface due to fast decreasing electron density, where coronal heating and the wind acceleration mechanisms are also occurring (Hundhausen et al. 1968a). Therefore, studying the ionization of the solar plasma can provide important information about energy deposition in the low SC.
There are two ways to compare models of solar wind heating and acceleration with observations: via in situ measurements of solar wind properties and via remote sensing observations of the solar wind source, especially through high-resolution spectra. The most stringent constraints are obtained when both types of data types are used; however, in order to carry out such a comprehensive comparison, it is necessary to use theoretical models encompassing the whole domain from the solar transition region all the way to the heliosphere. Such solar models have been developed to encompass this domain; a few early examples are given by Lionello et al. (2009), Downs et al. (2010, and the Space Weather Modeling Framework (SWMF; Tóth et al. 2012). In these models, semiempirical heating functions were used to heat the corona accelerating the solar wind via Alfvén wave pressure gradient. Comparison is focused on narrowband imaging observations in extremeultraviolet (EUV) and X-rays (e.g., Sachdeva et al. 2019), and in situ measurements of plasma properties. Only recently, synthetic spectra are being used for model validation (Szente et al. (2019), Shi et. al, 2021, submitted).
Observations by the Hinode spacecraft (De Pontieu et al. 2007) and the Solar Dynamics Observatory (McIntosh et al. 2011) have suggested that there is enough energy in the outward propagating magnetic fluctuations in the chromosphere, transition region, and low corona to obtain and maintain the coronal temperature at 1 MK. Following these results, advanced global 3D solar wind models were developed utilizing Alfvén wave turbulence as the engine for coronal heating and solar wind acceleration in a self-consistent way, such as those by van der Holst et al. (2014) and Mikić et al. (2018). In particular, the Alfvén Wave Solar Model (AWSoM; van der Holst et al. 2014) is an extended magnetohydrodynamic model that includes low-frequency, reflection-driven, Alfvén wave turbulence. AWSoM accounts for three different temperatures: isotropic electron temperature and the parallel and perpendicular proton temperatures. This model was later been combined with a threaded-field-line model (AWSoM-R, Sokolov et al. 2021) for heliospheric distances R S < R 1.1 R S for the purpose of providing time accurate simulation results from Sun to Earth faster than real time.
Many spacecraft have been providing in situ measurements of solar wind plasma properties; among recent ones, the Advanced Composition Explorer (ACE; Stone et al. 1998) and Ulysses (Wenzel et al. 1992;Balogh 1994;Marsden 2001) produced solar cycle-long solar wind data sets, which include velocity, magnetic field, ionization, and composition properties of the solar wind. Using input bulk speed (v), electron density (n e ), and electron temperature (T e ) along the trajectory of the solar wind, one can predict the wind charge state composition and compare it with the in situ measurements. Several charge state models have been developed that enable such a comparison following individual solar wind plasma parcels along their trajectory, intrinsically working in 1D (e.g., Gruesbeck et al. 2011 and references therein). The Michigan Ionization Code (MIC; Landi et al. 2012a), for example, combines v, T e , and n e profiles along the wind parcel's path with the ionization and recombination rate coefficients of the CHIANTI database (Dere et al. 1997;Del Zanna et al. 2021): a fourth-order Runge-Kutta method in combination with an adaptive step size are used to solve the ionization equations as a function of time, using ionization equilibrium at the wind source region as the initial condition. Shen et al. (2015) developed another charge state model using an eigenvalue method with adaptive time step.
More recently, charge state models started to be combined with 3D magnetohydrodynamic models as a post-processing tool. In these cases, the wind's v, T e , and n e are obtained along the flow lines of the 3D model and used as input for the 1D charge state calculations, and the results are compared to in situ charge state composition data measured by Ulysses/Solar Wind Ion Composition Experiment (SWICS; Landi et al. 2014;Oran et al. 2015). Results showed that the ionization rates were underpredicted compared to observations and it was suggested that the difference was due to the unaccounted suprathermal electrons. Lionello et al. (2019) integrated a fractional charge state code in the time-dependent 3D Magnetohydrodynamic Algorithm outside a Sphere (MAS) model with Alfvén wave turbulence for both the steady-state global wind and CMEs. Here, again the ionization rates were underpredicted, but the authors ascribed the discrepancy to the excessive wind speed of the 1D model. Artificially lowering this bulk speed provided a more favorable comparison. Coupling of 1D charge state calculations to partial results of a 3D model provides invaluable information regarding a plasma parcel with specific solar wind properties, but does not allow us to reach an understanding of how the solar wind evolves on a global level, and limits comparisons only to the few times where the calculations are made. On the contrary, a systematic 3D determination of the solar wind charge state composition evolution can open a window to the solar wind global evolution and also provide us with a tool that allows to (1) carry out comprehensive comparisons with measurements obtained from multiple spacecraft anywhere in the heliosphere, and (2) predict the environment that current in situ instrumentation such as those from the Solar Orbiter (Müller et al. 2020) and Parker Solar Probe (Fox et al. 2016) will face in their orbits.
Nonequilibrium effects can affect line emission close to the Sun (Landi et al. 2012b;Shi et al. 2019). The combination of modules predicting the EUV and X-ray emission from the 3D global model with global charge state calculation can contribute to improve the quality of the information obtained by comparing EUV and X-ray emission with narrowband and spectroscopic data. This motivated us to develop a newly implemented nonequilibrium ionization calculation as an integral part of the 3D AWSoM model, which already possesses both narrowband imaging and spectral calculation (Szente et al. 2019) capabilities.
The paper is organized as follows: We first describe the implementation of nonequilibrium ionization calculations for ions H-Zn into the AWSoM model in Section 2. Then we discuss the background solar wind obtained with AWSoM in Section 3. In Section 4, we analyze the results and discuss the freeze-in process along a select number of flow lines of various footpoints. In Section 5, we compare the model output with observations from the SWICS (Gloeckler et al. 1992) instruments on board the Ulysses and ACE spacecraft. We summarize our findings in Section 6.

Implementation
The charge state composition of the solar wind is determined by ionization and recombination due to inelastic collisions between free electrons and ions. Landi & Lepri (2015) used MIC to show that photoionization from background UV, EUV, and X-ray solar radiation also contributes to further ionize a few species. In the present work, we only consider collisional processes (radiative and dielectronic recombination, collisional ionization, and excitation-autoionization), and defer the implementation of photoionization to the future as these are less effective processes in forming the frozen-in charge state distributions (Hundhausen et al. 1968b). We also did not consider the first ionization potential (FIP) effect during implementation.
The temperature-dependent ionization and recombination rate coefficients are taken from the CHIANTI database. At every location of the solar wind trajectory, the local collisional ionization and recombination rates depend on the plasma electron density (which regulates how many collisions an ion undergoes) and electron temperature (which determines the efficiency of each collision at ionizing/recombining the colliding ion). The wind speed determines the time that the solar wind parcel spends at any given location: if this time is long enough, the plasma parcel reaches ionization equilibrium; however, Landi et al. (2012c) showed that once the solar wind is released from the source, it almost immediately departs equilibrium. Also, if the local electron density is lower than a certain threshold value, the probability of an ion undergoing collisional ionization and recombination becomes very small, so that this ion remains unperturbed. This threshold density is different for each species, but as the wind density monotonically decreases with distance, more and more species stop ionizing and recombining. Once all of them have stopped, the plasma ionization status freezes in and its charge state distribution does not evolve anymore. This status is attained at different distances for different wind parcels depending on the wind conditions. The resulting frozen in charge state distribution is the end product of the thermodynamic and dynamic history of the solar wind, and thus it provides invaluable information about the plasma evolution through the lower solar atmosphere from the wind source region to the freeze-in point.
In the present implementation of these processes, the SC and inner heliosphere (IH) modules of the SWMF provide the plasma background for the calculations. Both modules include an option to choose the physical model used to calculate the plasma parameters. In the present work, charge state calculations are performed when the model takes into account three types of temperatures: electron temperature, and anisotropic proton temperatures (relative to the local magnetic field). While only the electron temperature is used to calculate the ionization and recombination rates, having a model that decouples the thermodynamics of free electrons from the one of protons causes the resulting electron temperature values to be vastly different from those obtained with a one-temperature or even a two-temperature model (van der Holst et al. 2014).
The charge states are calculated throughout the threedimensional domains of the SC and IH components in a cellby-cell manner, similarly as done by MIC in one dimension, by solving the system of equation that regulates the plasma ionization and recombination processes: where T e is electron temperature, n e is the electron density, R m and C m are the total recombination and ionization rate coefficients, respectively, and y m is the fraction of the element in charge state m, so that The total ionization and recombination rate coefficients ( ) C T m e and ( ) R T m e are read from tables generated using CHIANTI; these values are tabulated on an electron temperature grid element by element.
The charge state distributions of ionization equilibrium for the boundary condition at the solar surface (1 R S ) and the initial condition throughout the whole domain in the beginning of the simulation are obtained by assuming equilibrium at the source: setting the right-hand side of Equation (1) to zero.
The results presented in this paper are calculated using the coronal abundances (Feldman et al. 1992).
The SWMF architecture naturally provided the possibility to implement the above calculations so that any three temperature AWSoM SC and/or IH simulations can be performed with the inclusion of charge state calculations in line without postprocessing and with no modification to the original model, regardless of whether they describe the steady-state SC or dynamic events (CMEs, jets, and blobs). At any time, the charge state distribution is calculated over the entire computational domain, providing a convenient tool to predict the values for any place along the orbit of past and current satellites, such as Ulysses, ACE, and Solar Orbiter. Doing the nonequilibrium ionization calculation in line has the advantage that the radiative cooling can be computed self-consistently in future model modifications.  There is also an asymmetry between the northern and southern poles: the plasma is warmer, denser, and slower above the north pole than above the south pole. These calculations are carried out under the assumption that the free-electron velocity distribution is Maxwellian, so that the CHIANTI rates, which are provided under this assumption, can be used. Such a distribution may not be appropriate for the solar wind (Montgomery et al. 1968;Ralchenko et al. 2007); the presence of non-Maxwellian tails of high energy electrons can indeed enhance the ionization rates and change the overall distribution (Cranmer 2014): when empirically included, they improved the AWSoM results compared to observations, see Oran et al. (2015). Also, we did not include photoionization, which affects certain ions (oxygen) more than others (carbon and iron), as shown by Landi & Lepri (2015). Both these approximations are expected to provide lower predicted ionization levels than observed especially in case of oxygen. Furthermore, our treatment of charge state evolution ignores the contribution provided by charge exchange between solar wind ions and hydrogen and helium atoms outgassing from circumsolar dust, which, as Rivera et al. (2020) showed, can enhance the abundance of He 1+ by orders of magnitudes over the values predicted by MIC. In addition, our treatment assumes that all heavy ions flow with the plasma at the same speed, not experiencing differential acceleration nor any interaction with the turbulence which energizes the background plasma.

Solar Wind Background
The solar wind background is provided by the SWMF's AWSoM model (van der Holst et al. 2022). The SWMF is an open source software containing multiple modules used for physics-based space environment modeling (Tóth et al. 2012). AWSoM is the model that is used in the the SC and IH modules.

Simulation Setup for the SC and IH
The SC component starts from a uniform parallel and perpendicular proton temperature and isotropic electron temperature of 50,000 K at the inner boundary. The proton number density at the boundary is N p = 2 × 10 17 m −3 : such a large, overestimated value prevents chromospheric evaporation in the same way as in Lionello et al. (2009). The initial condition for the solar wind is the isothermal Parker wind solution. The boundary and initial conditions for the charge states are ionization equilibrium as discussed in Section 2. The magnetic field at the inner boundary is prescribed via magnetograms: in our case we used the Global Oscillation Network Group (GONG; Harvey et al. 1996) magnetogram of Carrington rotation (CR) 2063 (between 2007 November 7 and 2007 December 4). Because of the consistent underestimation of magnetic field at 1 au we experience when using GONG magnetograms, we enhanced the radial magnetic field strength at the boundary (50,000 K temperature, 1 R e ) by a factor of 3.7, empirically selected after a trial-and-error procedure. The energy density of the outward propagating Alfvén waves is set via the Poynting flux S A : The solar atmosphere is heated by Alfvén wave turbulence and plasma is accelerated by the Alfvén wave pressure gradients. As previously mentioned, the simulation takes into account the proton temperature anisotropy relative to the local magnetic field, it calculates a separate electron temperature (three-temperature model). Also the model incorporates radiative losses from CHIANTI 8.0 using coronal abundances (Feldman et al. 1992) calculated under ionization equilibrium, Coulomb collisional heat exchange, and collisional and collisionless electron heat flux.
AWSoM uses two different grids for the SC and IH components. For the SC component we use a spherical grid starting at 1 R e and ending at 24 R e . The grid is stretched toward the Sun as to resolve the steep gradients in the transition region and low corona. To accurately resolve the charge state evolution near the solar wind source we need a fine grid close to the solar surface; the final number of grid cells used is about 23 million. For the IH component we use a Cartesian grid, which extends from −750 to 750 R S . The domain is such that the total number of used grid cells is about 40 million. Since at 24 R S the solar wind is already frozen-in, no further refinements were necessary. To obtain a steady-state solution we iterated the AWSoM equations for 200,000 steps in the heliographic rotating frame. Then we performed time accurate simulations started from this point of steady solution and followed several plasma parcel's evolution (see Section 4 for 50,000 s and in one case 80,000 s.

Model Validation: Narrowband Images and In Situ Bulk Wind Properties
To assess the quality of the global model calculation we first compared the plasma results to in situ observations from the WIND and Solar Terrestrial Relations Observatory (STEREO) spacecraft at 1 au (as extracted from NASA/Goddard Space Flight Center (GSFC)'s OMNI data set through OMNIWeb), as well as EUV narrowband images of the inner corona taken from the Solar and Heliospheric Observatory (SOHO) and STEREO spacecraft. During CR 2063 the two STEREO spacecraft locations allowed them to observe the Sun from a line of sight (LOS) ≈20°from the Sun-Earth direction, where SoHO and WIND were located. Figures 1-3   SoHO/EIT, observing from Earth, and STEREO-A/EUVI, observing the plasma that is rotating away from SoHO's view from the west limb. As this CR corresponded to the the very quiet minimum of solar cycle 23, no active regions were present in the field of view, while coronal holes were present at both poles, with the southern one being quite extended. The AWSoM predictions easily capture the presence of both holes, and also well reproduce the enhanced limb emission from the streamers at both limbs, as well as the region with decreased intensity in the southern hemisphere seen by STEREO-B at the central meridian. The color scale of each channel of each spacecraft is the same between observations and model, and it shows that AWSoM overpredicts solar emission at the limb in the quiet Sun, while at both poles the disagreement is reduced.
The comparison at 1 au for all three spacecraft are shown in Figure 4. For each spacecraft, the solar wind bulk speed, proton number density, isotropized proton temperature T = (2T ⊥ + T ∥ )/ 3, and magnetic field strength are shown. In all cases the simulation successfully reproduces the rapid change of plasma properties seen by WIND around 2007 November 20, after which the solar wind speed is overestimated significantly at all three locations. All other quantities are successfully reproduced, overall confirming the quality of the model predictions. Figure 5 shows the main coronal properties (plasma density, electron temperature, and plasma speed) involved in the charge state calculation along the meridional plane as seen from Earth on 2007 November 4. The configuration of the solar atmosphere is typical of solar minimum, with coronal holes at both poles where faster solar wind is accelerated, and a system of streamers in the equatorial region. The streamer belt as well as the current sheet are tilted from the ecliptic plane, especially at the east limb where a large and hot streamer extends southward at almost 40°inclination, while at the west limb the northward tilt is much lower. Also, a third, weaker structure, which the magnetic field configurations indicates to be a pseudo-streamer, is present in the north-east sector of the image, whose temperature however is much lower than the other two structures and barely reaches 1 MK in the entire field of view. Furthermore, the north polar coronal hole has warmer plasma than the southern one, and its temperature exceeds 1.5 MK at much lower altitudes than the southern coronal hole; also, the north pole wind speed exceeds 100 km s −1 at lower heights than in the south polar hole. The interplay of increased temperature and speed in the north coronal hole results in the competing effects of increased ionization due to the higher temperature, and in a lower ionization due to the smaller time span spent in the high density regions: the consequences for the charge state evolution of such different properties between the two coronal holes will be discussed in Section 4.

Nonequilibrium Effects versus Equilibrium Charge States
Figures 6 and 7 show selected charge state ratios from N, Ne, Si, and S in the same plane as Figure 5. Figures 8, 9, and 10 show selected charge state ratios from C, O, and Fe in the same plane as Figure 5 along with the equilibrium ionization results in the same plane. N 5+ /N 6+ , N 6+ /N 7+ , Ne 6+ /Ne 7+ , Ne 7+ /Ne 8+ , C 4+ /C 6+ , C 5+ /C 6+ , O 5+ /O 6+ etc. ratios are expected to decrease in a hotter, denser, slower moving plasma, while the S 11+ /S 10+ , O 7+ /O 6+ , Fe 11+ /Fe 10+ ratios and the average Fe charge are expected to be larger.
The different properties of the two polar coronal holes result in different charge state ratios, with the north pole consistently showing higher ionization in all ratios, indicating that the higher temperature leaves more lasting signatures than the shorter time spent in the inner corona due to the larger acceleration. It is worth noting that oxygen charge state ratios show a stronger evolution in the inner corona at the north pole than in the south pole, where they seem to freeze in at lower heights than in the north pole.
The nonequilibrium solutions for C, O, and Fe are compared to the equilibrium solutions in Figures 8, 9, and 10, while the individual nonequilibrium to equilibrium ion abundance ratios for carbon and oxygen are shown in Figures 11 and 12.
Carbon ionization solutions show that in the fastest regions (polar coronal holes) nonequilibrium effects result in larger fractions of carbon to be C 4+ and C 5+ than expected at equilibrium at distances of only a few tenths of radii, while C 6+ is predicted to be significantly lower than at equilibrium at all distances. In the streamer belt nonequilibrium solutions are closer to the equilibrium, while in the eastern located pseudostreamer the nonequilibrium solutions are the opposite of what is observed in the fast wind: there are less C 4+ and C 5+ and more C 6+ in the nonequilibrium solution. The interplay of these results is shown in Figure 11: C 4+ /C 6+ and C 5+ /C 6+ both show that the plasma is less ionized in the nonequilibrium solution in the fast polar wind, and more ionized in the pseudostreamer. Departure from equilibrium is less in the slow-wind streamer belt region.
In the case of oxygen the departure from equilibrium is overall smaller in all regions. The fast wind originating from the polar coronal holes is predicted to have more O 5+ , O 6+ and less O 7+ , while nonequilibrium is also affecting the oxygen charge state composition at the edge of streamers in a different way than in their centers. These differences are then propagated to the charge state ratios: the values in the nonequilibrium and equilibrium solutions are shown in Figure 12; the nonequilibrium values for the fast wind indicate lower values for both the O 5+ /O 6+ and O 7+ /O 6+ ratios, suggesting that the wind is more concentrated in the dominant, He-like O 6+ stage and did not have the chance of ionizing further. In the slow wind we find ratios closer to 1, which means the departure from equilibrium are less pronounced, but still significant. In the streamer on the east we see increased ratios compared to the equilibrium solution. It is interesting to note that nonequilibrium effects are found also in closed field structures, due to the combination of the presence of flows with lower plasma density.
In the streamer belt, the charge state ratios and the average Fe charge state indicate higher ionization than in the solar wind, as expected from a plasma with temperature exceeding more than 2 MK. Overall, Figures 11 and 12 indicate that speed-induced departures from ionization equilibrium are  . Radial magnetic field component is colored on the 1 R s sphere of the solar surface, and #1-10 flow lines of plasma parcels are indicate the path the parcels travel during the simulation. The field of view is 4 R e . The enumerated field lines' footpoint on the solar surface correspond to the starting point of the plasma parcel's history followed as shown in Figures 14 and 15, with the exception of field line #10, whose corresponding plasma parcel did not depart from the solar surface beyond 1.022 R e , even within 80,000 s of simulated real time. We show the southern (left), northern (middle) poles, and the eastern (right) side of the solar sphere to show the footpoints of the plasma parcels followed. present both in closed and open field structures, though their values is larger in the latter. Departures from equilibrium in closed field structures have been predicted in the past, either as a result of siphon flows (e.g., Spadaro et al. 1990) or nanoflares (Bradshaw et al. 2012). In both cases, the changes due either to speed-induced or nanoflare-induced variations in temperature were faster than the speed with which the plasma could adapt to the new temperature. Figures 11 and 12 indicate that such variations can be widespread in closed field structures, especially when the electron density is lower and thus at larger heights, although in the case of AWSoM departures from equilibrium are entirely due to the effect of speed.
It is important to note that these departures occur well within the range of heights covered by past and current highresolution spectrometers and narrowband imaging instruments, and therefore may be expected to affect the analysis of spectral line intensities of the lines emitted by each of these ions. Usually, ionization equilibrium is assumed throughout the inner SC, but Shi et al. (2019) discussed the effects of windinduced departures from equilibrium on coronal plasma diagnostics, concluding that these effects lead to significant changes in the measured plasma elemental abundances. The key parameters in the corona, whose temperature ranges from 1-3 MK, are the electron density and the speed. In most closed  structures at low heights, the speed v is of the order of 10 km s −1 or less and the electron density n e densities are larger than 10 7 cm −3 . Assuming an isothermal loop around 10 5 km long at a temperature in the 1-3 MK range lying within the field of view that current EUV imaging instruments reach (i.e., within ≈1.3 R e ), these values correspond to values n e × t larger than 10 11 cm −3 s: according to Smith & Hughes (2010) such values imply that the plasma is likely to be close to equilibrium, especially for Fe, which is the main contributor to the observed emission in the coronal channels of EUV imagers. However, in the presence of flows starting from the chromosphere and traveling into the corona, the e-folding time toward approaching equilibrium becomes temperature dependent, so that it is difficult to estimate departures from equilibrium a priori, and how these will propagate into the emission observed by such instrument. Landi et al. (2012b) carried out this calculation for the fast wind using the fast wind model from Cranmer et al. (2007), finding that with the exception of O, the departures from equilibrium of all elements were limited: this would suggest that narrowband imagers should be close to equilibrium. Still individual loops might host larger speeds in the transition region, which could make these effects larger: in order to draw definitive conclusions, a full nonequilibrium calculation of the solar spectrum is necessary. This work will be pursued in a future paper, where nonequilibrium ionization will be coupled to SWMF's SPECTRUM module to calculate nonequilibrium line intensities.
Another example is the difference in the predicted O 5+ nonequilibrium/equilibrium abundance ratio in the core and in the leg of both the west and east streamers present in the model. A nonequilibrium effect leads to larger O 5+ values in the leg than in the center, which translate into larger line intensities in the streamer legs than in the center. Raymond et al. (1997) reported brighter streamer legs than centers observed by the SOHO/Ultraviolet Coronagraph Spectrometer (UVCS) instrument (Kohl et al. 1995) with the O 5+ bright 1031.9 and 1037.6 Å lines: using the standard assumption of ionization equilibrium, they concluded that the oxygen element was depleted in the streamer center relative to the streamer legs. The present result suggests that such differences in O 5+ line brightness may be in part due to departures from equilibrium. Specifically, the predicted higher O 5+ abundance, once factored in, decreases the oxygen elemental abundance required to account for the observed intensity. Assuming that the present streamer is a reasonable representation of the structure studied by Raymond et al. (1997), the present result indicates that the oxygen elemental abundance is even more depleted in the corona than reported in Raymond et al. (1997); also, the difference between the leg and core is less pronounced. Still, Figure 12 reports nonequilibrium effects that are unlikely to account for the entire discrepancy found by Raymond et al. (1997); also, LOS integration effects may even reduce the overall effect on measured intensities. Our results merely indicate that speedinduced nonequilibrium effects can influence the analysis of line Figure 15. Oxygen and carbon ion charge state compared to the final, frozen-in value is shown for the plasma parcels traveling along paths #1-9 as shown in Figure 13. The difference observed in the different freezing-in locations are due to the difference of temperature, density, and velocity distribution along the path, as shown in Figure 14. Carbon ions are plotted with solid, oxygen ions with dashed lines.
intensities even in closed magnetic structures, and that future studies that couple the present implementation with SPEC-TRUM are necessary to thoroughly assess the impact of nonequilibrium on spectroscopic diagnostics.

Individual Wind Flow Lines
In order to study the evolution of the solar wind along individual flow lines, we selected 10 points on the solar surface where the plasma is likely to escape (i.e., make it out to 5 R e ) and followed the plasma's path in time looking at how the charge states in the plasma parcel react to the changes in the ambient solar wind. The flow lines of the parcels are shown in Figure 13 and their characteristics are listed in Table 1, where we also include the distance and speed reached after 50,000 s of simulated real time. We selected these flow lines identifying footpoints where the magnetic field lines are open and we expect the plasma to depart from the solar body, and allow us to sample different physical situations: center and edges of both polar coronal holes, including locations close to the helmet streamer or the open-closed boundary, as well as low-latitude regions close to the current sheet. We looked to the evolution with distance of the ratio of the carbon (C 4+ − C 6+ ) and oxygen (O 5+ − O 7+ ) charge states to their frozen-in value.
The velocity, density, and temperature profiles along the path for flow lines #1-9 are shown in Figure 14. #10 did not travel beyond 1.022 R e , even after 80,000 s, and thus it is not shown there. As seen, #1 and #5 are much slower parcels than the average; they originate from close to the ecliptic plane, while #8 and #9 are the fastest ones, both from the southern coronal hole. #10, the slowest parcel, originates from the southern coronal hole boundary: with the exception of #1 and #10, all parcels traveled beyond 3 R e . #1 reached a (presumably) frozen-in state already at about 1.4 R e . Figure 15 shows the ratio between the charge state values and their frozen-in value along each flow line: results of the fast wind coronal holes centers are qualitatively similar to those obtained in 1D by Landi et al. (2012a Figure 7), using the fast wind from the coronal hole model of Cranmer et al. (2007). Plasma traveling along flow lines corresponding to open magnetic field freeze-in below 3 R e , although the precise height changes from ion to ion and flow line to flow line, with #1 freezing in already at 1.5 R e despite coming from a lowlatitude source region. Note, that #1 does not reach beyond the height 1.54 R S , spending all 50,000 s in higher density plasma than #5, whose speed magnitude is comparable, but seemingly more radially oriented, reaching about 3 R S during the same time period. What is unforeseen based on results shown by Landi et al. (2012a), is that along certain flow lines (#2, #3, and #6 are from the edge of coronal holes, and #5 is low latitude and not related to the coronal hole open field) the ionization history passes the frozen-in value, and then later resumes evolving. On the contrary, plasma properties from coronal hole centers evolve in a smoother manner and so the ionization states relax in a monotonic manner to the final, frozen-in values, unlike in the above-mentioned examples.
Flow line #6 is the line showing oscillation-like behavior, with the charge state ratios changing in a manner reminiscent of numerical instability. However, such a variability is due to a different cause, namely, the plasma parcel is traveling along a helmet streamer boundary as can be seen in Figure 16, and the close proximity of the higher density plasma affects the sensitive charge state calculation. The SWMF uses the Block-Adaptive-Tree Solar wind Roe-type Upwind Scheme Figure 16. Field line corresponding to footpoint #6 is marked with magenta color, and neighboring magnetic field lines are plotted in black, the solar body is colored respective to the radial magnetic field. Parcel #6 travels along the boundary of a helmet streamer, which explains the unusual fluctuations seen in the corresponding image of Figure 15. (Powell et al. 1999), in which we use a second order scheme for this simulation. Due to the finite grid resolution, the second order scheme occasionally samples from the high density streamer right next to the flow line resulting in the perceived oscillations. It is expected that at higher resolution such behavior would be greatly decreased, and confined to an even narrower range of location close to the real streamer boundary.

Comparison with Observations
We compare charge state distributions and charge state ratios for multiple ions to what it was observed along the paths of Ulysses spacecraft during the 2007 February 15 UT 00:00:00 and 2008 January 15 UT 00:00:00 time period, and of the ACE spacecraft during the time span of CR 2063: 2007 November 4 UT 09:59:00 to 2007 December 4 UT 09:59:00. While ACE was observing in the ecliptic plane for the whole duration of CR 2063, Ulysses was undergoing its third polar pass, which however lasted for an entire year so the boundary condition based on the radial magnetic field of CR 2063 only captures a small part of Ulysses' third pass. Still, a comparison between the data collected during the entire Ulysses pass and CR 2063 predictions along the Ulysses path can provide a meaningful assessment of the model performance because this CR was extremely quiet, and during it there were only 2 days affected by CME events (2007 November 19 UT 23:00 and 2007 November 20 UT 12:00) according to the Richardson-Cane CME list (Cane & Richardson 2003;Richardson & Cane 2010) available online at http://www.srl.caltech.edu/ACE/ASC/ DATA/level3/icmetable2.htm. Thus, CR 2063 can be taken as a proxy for the year-long Ulysses pass, which occurred during solar minimum.
Level 2 ACE/SWICS data was provided by the ACE Science Center 1 while the Ulysses/SWICS data was obtained from the ESA-NASA Ulysses final archive at http://ufa.esac. esa.int/ufa/#data. Since the two CME-affected days alone do not cause any significant bias in the distribution, we did not exclude them from the comparison with ACE.

Comparison with Ulysses Measurements during the Third Solar Pass
We compared the AWSoM predictions along Ulysses' entire third polar pass, which covered almost the entire range of latitudes between the south and north poles. This comparison allows us to assess AWSoM predictions outside the ecliptic, and is particularly important for the measurements that Solar Orbiter will take in the future.
However, Ulysses completed its third polar pass between 2007 February 15 and 2008 January 15: this time period exceeds the boundaries of our simulation, which only include CR 2063, which roughly corresponds only to the month of 2007 November. As the magnetic boundary we used for this CR changed significantly during the rest of the year, the comparison between the present results to the Ulysses measurement is not consistent, in the sense that only a small part of the Ulysses data sets was effectively observed during CR 2063. Still, since the Ulysses pass took place during the depth of the very weak minimum of solar cycle 23, the Figure 18. Carbon, oxygen, iron, and magnesium charge state ratios and iron average charge state are compared between ACE data (black) and model output (green) for the same time period of CR 2063. configuration of the Sun was evolving slowly so that a comparison is still qualitatively meaningful, as the solar structures predicted by AWSoM are qualitatively similar to those that took place during the year of the Ulysses pass. Figure 17 shows the comparison between the predicted and measured C 6+ /C 5+ , O 7+ /O 6+ and average charge state of Fe, computed from Fe 6+ to Fe 16+ , for all latitudes sampled by Ulysses during its polar pass. Predicted values are taken from the trajectory that Ulysses followed. Figure 17 shows that the variability and higher charge state distribution of the slow wind is captured at the right nearequatorial latitudes, confirming that the overall structure of the Sun was not very variable, and that AWSoM is capable of accounting for the short term changes between faster and slower wind streams occurring around the streamer belt. As far as the absolute values of the charge state ratios of C and O, they are correctly reproduced: given the relatively low spatial resolution that memory limitation imposes on global models, the small-scale details are lost, but the solar wind ratios are in the correct range, indicating that overall the predictions of the model are accurate. The average charge of iron is overestimated at the poles, more on the north pole than on the south pole. Still, the value of the discrepancy is limited.
For the same time period Oran et al. (2015) predicted lower charge state distribution than observed: their result prompted Landi & Lepri (2015) to evaluate the effect of photoionization, finding that it could be significant but not enough to account for the discrepancy. This was due to the fact that during the minimum of solar cycle 23 the X-ray and photoionizing flux from the SC was too low to enhance the overall ionization; in fact, in the present model has been able to reproduce the observed charge states even without the inclusion of this additional process.
At the north pole, the predicted charge states are overestimated. This may be due to some improvement needed by the model in treating the north pole, as during CR 2063 Ulysses was sampling the latitude range 57°.62-70°.75 so that the predictions should be most accurate (this period is marked by vertical lines in Figure 17, with about a 6 day extra margin to cover the time plasma travels to the approximately 1.6-1.7 au distance). It is Figure 19. Carbon charge state ratio distribution is compared between ACE data (black) and model output (green) for the same time period as CR 2063.
worth noting that magnetic field signatures at high latitudes are very weak, resulting in uncertain measurements, which lead to imperfect boundary conditions for the SC model: these uncertainties might contribute to these discrepancies.

Comparison with ACE Measurement at the Ecliptic
The comparison with ACE measurements provides a more direct assessment of the predictive capability of the model, as it represents data corresponding to the exact time period of boundary condition of the simulation, namely, CR 2063, although such assessment is limited, considering that it focuses AWSoM predictions in the ecliptic plane only. The results are shown in Figure 18, where select measured charge state ratios are compared with AWSoM prediction as a function of time. Synthetic data is sampled once every 1800 s, while ACE data was sampled once about every 7200 s. In the case of the comparison of the Ulysses passing time period, the ACE data is represented by one data point per day. We focused on the C 6+ /C 4+ , C 6+ /C 5+ , O 7+ /O 6+ , Mg 10+ /Mg 9+ , and Fe 11+ /Fe 10+ charge state ratios, and on the average Fe charge. The average Fe charge was computed including all charge states from Fe 6+ to Fe 20+ , as the others had too few counts in the available data set to be considered reliable.
There are two main things to notice from these results. First, all predicted charge states are in the same range as observed by ACE. This is very encouraging as it indicates that the evolution of the solar wind is also correctly captured by AWSoM in the ecliptic plane, where the solar wind coming from coronal hole edges is modeled. Second, the simulation seems to be even able to predict, with a few days' delay, the decrease in the carbon, oxygen, and magnesium charge state ratios occurring at around 2007 November 22, meaning that a large scale structure corresponding to the faster wind stream seen by OMNI in Figure 4 at around 2007 November 20-25 produces compositional signatures that AWSoM can successfully capture. Interestingly, the lower carbon, oxygen, and magnesium charge state ratios do not correspond to lower charge states for iron. AWSoM actually predicts an increase in those values, while the observations seem to suggest a steady value: although AWSoM seems to overestimate the iron ionization status, it is nevertheless capable of capturing the different behavior between these elements. Since the resolution of global models is not large enough to enable a direct comparison between individual wind streams, we compared the distribution of predicted charge state ratios with observed ones. The direct comparison for CR 2063 is shown in Figures 19 and 20 for the C 4+ /C 6+ , C 5+ /C 6+ , O 5+ /O 6+ , and O 7+ /O 6+ charge state ratios. AWSoM predictions are in the same range as the observations, although they tend to be more sharply clustered around some values than the observations suggest. There is one exception: the O 5+ /O 6+ ratio, which is predicted to be one order of magnitude lower than observed. As noted by Landi et al. (2012c), the large ionization rates from O 5+ to O 6+ cause this ratio to freeze at larger heights than the other oxygen ions, and mimic ionization equilibrium until the plasma freeze-in, so that this ratio is expected to have values resembling the SC in the 1-3 MK temperature range, which are predicted to be ≈2-4 × 10 −3 consistent with AWSoM results. On the contrary, ACE observations are more typical of colder plasmas, at altitudes far below the oxygen freeze-in height. This discrepancy, however, may be due to another process recently discovered by Rivera et al. (2020) to increase the He 1+ /He 2+ ratio orders of magnitude from the expected value, namely, recombination induced by wind particles colliding with neutral hydrogen and helium outgassed from circumsolar dust. The large O 6+ -H charge exchange rates and the relative abundance of O 6+ make this process a promising candidate: work is in progress to study whether this process is responsible for this behavior (Y. J. Rivera et al. 2021, in preparation).
The question, however, is why is only O 5+ affected by charge exchange and not O 6+,7+ , or the carbon ionization stages. The reason may be that oxygen is a peculiar element because it is mostly concentrated in the He-like O 6+ stage, and the abundances of O 5+ and O 7+ are much lower. This means that small amounts of O 6+ recombining into O 5+ by charge exchange can enhance the abundance of the latter element by a large factor, without significantly changing the overall abundance of O 6+ . In the same way, charge exchange can significantly decrease O 7+ , without altering O 6+ . However, these trends in O 5+,7+ are directly counteracted by photoionization (also not taken into account by the current implementation). It is possible that our neglect of photoionization is behind our underestimate of the O 7+ /O 6+ ratio because for O 7+ photoionization is more Figure 21. Carbon charge state ratio distributions are compared to ACE data (black) as simulated (green) for a year, using the magnetic background of CR 2063. important than charge exchange due to the very large abundance of O 6+ . On the other hand, photoionization is not enough to counterbalance the charge exchange contribution to O 5+ , again because the large abundance of O 6+ causes the charge exchange to dominate. For carbon, charge exchange and photoionization may be balancing each other because this element is not dominated by a single charge state like oxygen, and the major charge states (C 4+ to C 6+ ) have comparable abundances, and Landi & Lepri (2015) showed that photoionization is less important for carbon than for oxygen. However, a detailed calculation is necessary to confirm or improve this scenario, which we defer to a future paper.
Distributions of carbon and oxygen charge state ratios (C 4+ /C 6+ , C 5+ /C 6+ , and O 7+ /O 6+ ) and the average charge state of oxygen (calculated using charge states from 5+ to 8+) have also been compared to the distribution of the measurements taken by ACE for an entire year. The comparison is shown in Figures 21 and 22, and shows that AWSoM is capable of reproducing the range of values under all conditions during an entire year. Oxygen is under-ionized, with the highest peak of the predicted distribution indicating a ratio a factor 2 or less than observed. This difference can be ascribed to our model having neglected photoionization: Landi & Lepri (2015) showed that in 2007 photoionization was able to increase the predicted O 7+ /O 6+ by around 1.3-1.8 (see their Figure 3) for typical ionizing fluxes observed at the ecliptic by Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED)/Solar EUV Experiment (SEE). It is worth noting that such a photoionizing flux is expected to be lower for the wind observed by Ulysses (see Section 5.1) and the portion of the solar disk affecting high-latitude wind will be dominated by polar coronal holes, whose EUV and X-ray flux is significantly lower than the quiet Sun values affecting the wind observed in the ecliptic plane.

Summary
We described the implementation of nonequilibrium charge state ionization in AWSoM, which allows determining the nonequilibrium ionization distribution everywhere in the SC/IH computational domains of the SWMF, combining CHIANTI ionization and recombination rate coefficients with AWSoM's predicted plasma electron temperature, density, and speed in every cell. We carried out a simulation for CR 2063 during the minimum of solar cycle 23, extending the results to 1 au, and compared them with in situ measurements of charge state composition from the SWICS instrument on both the ACE and Ulysses spacecraft. We directly compared ACE measurements for CR 2063, and the data for the entire Ulysses third polar pass, even though it extended beyond the CR we modeled. Predicted charge state distribution compared favorably with observations, demonstrating the accuracy of both the implementation and of the AWSoM model's predictions. Charge state calculations in the 3D solar wind model are very useful to study the evolution of the solar wind plasma in regions inaccessible to both remote sensing and in situ observations, where the heating and acceleration mechanisms are still active.
During the comparison with observations, we discovered that solar circumstellar dust might play a role in altering the charge state composition of oxygen (Y. J. Rivera et al. 2021, in preparation). The agreement with observations improves on the results obtained with previous versions of the code providing agreement with in situ measurements, while earlier comparison showed too low predicted charge states, indicating improvements in the AWSoM's model performance. Also, this new update to the AWSoM model provides an excellent tool to predict solar wind conditions anywhere in the heliosphere, where retired, future, or current (Parker Solar Probe, Solar Orbiter, etc.) spacecraft are (or were) located. Overall, the reported discrepancies in charge state composition are similar to other nonequilibrium ionization modeling efforts, such as in Oran et al. (2015), where the suprathermal electron population was used, or Shen et al. (2017) where the validation was not directly compared with in situ measurements, but by deriving the synthetic emission from the region of interest and compare with UVCS observations. Future development steps will be the inclusion of photoionization in the equation for ionization and recombination, whose effects were small for the CR analyzed in the present work, but are expected to be more significant during solar maximum. Also, the coupling of the present nonequilibrium charge states with he post-processing tool SPECTRUM will allow the prediction of nonequilibrium line intensities: nonequilibrium effects are expected to be important while studying the emission of coronal holes and other plasmas (Shi et al. 2019). Also Landi et al. (2012b) showed that wind-induced departures from ionization equilibrium can affect radiative losses in the innermost solar atmosphere. As the AWSOM model incorporates radiative loss calculations, and these departures are significant in the transition region, it would be beneficial to conduct a study on how the induced radiative loss changes the energy deposition at a global level.