3D MHD Time-dependent Charge State Ionization and Recombination Modeling of the Bastille Day Coronal Mass Ejection

Heavy ion signatures of coronal mass ejections (CMEs) indicate that rapid and strong heating takes place during the eruption and early stages of propagation. However, the nature of the heating that produces the highly ionized charge states often observed in situ is not fully constrained. An MHD simulation of the Bastille Day CME serves as a test bed to examine the origin and conditions of the formation of heavy ions evolving within the CME in connection with those observed during its passage at L1. In particular, we investigate the bimodal nature of the Fe charge state distribution, which is a quintessential heavy ion signature of CME substructure, as well as the source of the highly ionized plasma. We find that the main heating experienced by the tracked plasma structures linked to the ion signatures examined is due to field-aligned thermal conduction via shocked plasma at the CME front. Moreover, the bimodal Fe distributions can be generated through significant heating and rapid cooling of prominence material. However, although significant heating was achieved, the highest ionization stages of Fe ions observed in situ were not reproduced. In addition, the carbon and oxygen charge state distributions were not well replicated owing to anomalous heavy ion dropouts observed throughout the ejecta. Overall, the results indicate that additional ionization is needed to match observation. An important driver of ionization could come from suprathermal electrons, such as those produced via Fermi acceleration during reconnection, suggesting that the process is critical to the development and extended heating of extreme CME eruptions, like the Bastille Day CME.


Introduction
The Bastille Day coronal mass ejection (CME) and associated X5.7 flare was one of the most energetic events of Solar Cycle 25 (Manoharan et al. 2001).The CME occurred on 2000 July 14 at 10:03-10:43 UTC, originated from a filament eruption in NOAA Active Region 9077 on the northwest part of the solar disk, and drove a halo CME directed at Earth.The CME and flare have been studied extensively from remote and in situ observations throughout the solar system.There have been several studies investigating the event and its impact on the heliosphere in the context of the injection and acceleration of energetic particles (Young et al. 2021), ion production rate induced by cosmic rays (Mishev & Velinov 2018), its magnetic field topology and evolution (Wang et al. 2001), its posteruption loop arcade (Reeves & Warren 2002), and electron temperature and speed profile (Reginald et al. 2018).Most recently, it has been studied in connection with its impact on Earth's magnetosphere and associated geomagnetic storm (Brandt et al. 2001;Kubota et al. 2017) and its overall impact on space weather and the heliosphere (Linker et al. 2016;Török et al. 2018).
The Bastille Day CME was one of the most extreme and fastest eruptions, reaching a speed of 2500 km s −1 in the corona and arriving at L1 on 2000 July 15 at ∼12: 00 UT, just over 1 day after its eruption, traveling at >1100 km s −1 following four successive CMEs.Török et al. (2018) simulated the Bastille Day eruption from Sun to Earth, enabling a connection of the ejecta in the corona to in situ measurements near Earth.Modeling the eruption and heliospheric propagation is critical to understanding how the CME disrupts the heliosphere and crucial to tracing the flow of energy that powers the eruption as the plasma is heated and accelerated from the Sun.
In particular, the heavy ion composition measured in the heliosphere can provide useful information regarding the thermal energy input and on conditions during the initial stages of the eruption and liftoff connected to the flow of energy.This is because the ionization and recombination governing the charge states in the plasma are primarily governed by collisions between electrons and ions, making the ionization state sensitive to the electron temperature, density, and bulk flow of the escaping plasma.As the plasma expands and becomes too tenuous to sustain ionization processes at larger heliocentric distances, e.g., out to 5 R e in the solar wind and 15 R e in CMEs, the ionization and recombination stop (Landi et al. 2012;Shen et al. 2017;Boe et al. 2018;Rivera et al. 2019b;Gilly & Cranmer 2020).After this distance, the charge states measured in the solar wind will retain properties of the plasma state during the eruption that can be linked to the energy deposition at the main phase of plasma heating and acceleration (Rakowski et al. 2007;Gruesbeck et al. 2011).The so-called freeze-in distance is the radial distance at which ions no longer change, and they can be distinctively imprinted and connected to different stages in the plasma's thermodynamic evolution.Together, observations of the full charge state distributions from ions of several elements can provide sensitive and detailed plasma diagnostics of the eruption.A comprehensive understanding of the formation of the ionization states congruent with other plasma and field properties can provide means of rigorously quantifying energy release and transport throughout the flare and ejecta to constrain the principal dissipation mechanisms at play.
There have been several implementations of time-dependent ionization simulations that examine and constrain CME heating.Lynch et al. (2011) examined the evolution of charge states in a single-fluid MHD numerical simulation of a CME by tracing a collection of plasma parcels in a flux rope from its eruption in the corona to ∼20 R e in a CME driven by flux cancellation and through the magnetic breakout scenario.The study finds that enhanced ionization states in connection with observations were linked to the slower CME that resulted from the plasma being connected to the flare current sheet, where the highest temperatures reached were ∼10 MK, longer than in the fast CME.In addition, the temperature profile of the cavity region showed continuous heating through 14 R e that was attributed to continual breakout reconnection.The study finds spatial variation of the ionization stage of plasma in both cases, with the highest O 7+ /O 6+ and C 6+ /C 5+ at the back/center of the flux rope as a direct consequence of reconnection heated plasma.That work suggested that the high charge states associated with CMEs are formed during the earliest stages of the eruption from flare heating and the ions do not have a chance to recombine.
Time-dependent ionization is also an important diagnostic for remote sensing observations.It has been applied to UV and EUV observations to investigate the energy budgets of CMEs (Akmal et al. 2001;Lee et al. 2009;Murphy et al. 2011;Kocher et al. 2018;Wilson et al. 2022) in a manner analogous to the present study.It also constrains the parameters of the plasma in the region where the solar wind originates (Landi et al. 2014;Oran et al. 2015;Shen et al. 2017).Perhaps the study most relevant to this investigation is the analysis of the changing ionization states of Fe from AIA observations of a coronal shock wave driven by a CME (Ma et al. 2011).
In this work, we examine the energetics of the Bastille Day CME through the radial evolution of charge states evolving within the 3D CME structure.In particular, the work aims to provide insight to the heating required to reproduce the observed heavy ion quantities measured by ACE/SWICS during the passage of the CME.We use a new version of the MAS MHD code that incorporates a time-dependent ionization calculation (Lionello et al. 2019) to run a CME model using the setup parameters that Török et al. (2018) used to simulate the Bastille Day event.A novel aspect of the simulation was that the eruption is driven from a magnetically stable flux rope configuration at the corona producing more realistic and physical initiation.Additionally, the eruption is propagated through interplanetary space, allowing for a direct connection between coronal dynamics and heliospheric observations.For the present study, the nonequilibrium ionization code is implemented during the CME's coronal evolution prior to where ions freeze in to compare the simulated ions to charge states measured during its passage at ACE.Through the examination of the ionization and recombination history of each ion within the ejecta, we can investigate the conditions of their formation and determine their coronal source.We follow the progression of the heating throughout its 3D structure to determine how ions are formed and evolve in association with where heating takes place.The study brings important insight to the origin of unique ion signatures, beyond those of the Bastille Day observations, that are common to CMEs, such as the bimodal Fe distributions.
The paper is organized as follows: Section 2 describes the observations, Section 3 describes the MAS model, Sections 4 and 5 present the analysis and discussion, and Section 6 summarizes the main conclusions.

Observations
As discussed, the charge states measured outside of the freeze-in altitude remain imprinted with heating and cooling experienced prior to that point, enabling a detailed examination of the plasma's thermal history at critical stages of the eruption.To gain insight to the initial conditions and extended heating reflected in the CME ion signature, we examine the charge state measured during the passage of the Bastille Day CME near Earth to compare to the simulated charge states from the modeled eruption.
Data collected from ACE at L1 include heavy ion composition data from the Solar Wind Ion Composition Spectrometer (SWICS; Gloeckler et al. 1992) and proton measurements from the Solar Wind Electron, Proton, and Alpha Monitor (SWEPAM; McComas et al. 1998) on the Advanced Composition Explorer (ACE).SWICS is designed to measure the thermal heavy ions in the solar wind, and its principle of operation is described in detail in Gloeckler et al. (1998).SWICS allows for the identification of a large range of ions of He, C, N, O, Ne, Mg, Si, S, and Fe, in the E/Q range of 0.49-60 keV e −1 .We also include magnetic field measurements from the Magnetic Field Experiment (MAG; Smith et al. 1998).
Figure 1 shows a stack plot of in situ observations captured by ACE.The Bastille Day eruption reached L1 in approximately 1 day following four CMEs.The period plotted contains some material from the previous CME, which also contained elevated C, O, and Fe charge states prior to the arrival of the Bastille Day event.The boundary for the CME associated with the Bastille Day eruption is shown between the red vertical lines, where the solid line indicates the initial encounter as indicated by the shock and the trailing boundary is identified by the dashed red line.The boundary is selected using the list from Richardson & Cane (2010), which has been expanded using the online list. 5Panels (a)-(c) display the heavy ion charge state distributions for C, O, and Fe.The color bars show the relative abundance of each individual charge state of a given element.Panel (d) shows the charge state ratios of O 7+ /O 6+ and C 6+ /C 5+ .Panel (e) shows the elemental abundances of Ne, N, C, S, Si, Fe, and Mg compared to O. Panels (f)-(h) show the proton temperature, bulk speed, and density, respectively, where the black curves show the highest-resolution data available and the red curves are the time-averaged, lowerresolution measurements illustrating the general profiles from SWEPAM.Lastly panel (i) shows the magnetic field components and magnitude.As summarized in Manoharan et al. (2001), the Bastille Day event is characterized by a sheath region at the front associated with shocked plasma, as indicated by the large fluctuations in the magnetic field components and enhancement of the proton temperature.The sheath is followed by a magnetic cloud as shown in the rotation in the magnetic field components, mainly B n , and the enhancement of the magnetic field magnitude (Wang et al. 2001).The magnetic field B r component was directed such that it was ideal for dayside reconnection with Earth's magnetosphere andgenerated a global geomagnetic storm (Lepping et al. 2001).
The first half of the CME is made up of the most highly ionized heavy ions, panels (a)-(c), in particular as demonstrated by the Fe charge state distributions between day of year (DOY) ∼197.6 and 198.42.The Fe distribution is peaked at Fe 16+ across this period, with an 〈Q Fe 〉 (as indicated by the red curve) larger than 16 and visible contribution from Fe ions spanning from 10+ all the way up to 24+ (the maximum charge state resolved by SWICS).These extremely high Fe charge states are not commonly observed in CMEs or the solar wind, providing a unique opportunity to study the energetics involved in producing these observations.The second part of the ejecta, after the passage of the main flux rope, is characterized by a typical CME feature; a bimodal Fe charge state distribution.The bimodal Fe profile hasa peak at Fe 16+ and at a second peak at a lower ionization state, between Fe 8+ to Fe 11+ (Gruesbeck et al. 2012).The double-peaked feature shows some spatial variability across this structure, where the population of lower charge states dominates around DOY 198.42 and gradually shifts from a peak around Fe 8+ to a peak around Fe 10+ .Gradually this low charge population becomes less dominant compared to the population around Fe 16+ that dominates near the rear end of the CME and even extends beyond the end of the CME compared to the rear boundary defined by Richardson and Cane's ICME list (Richardson & Cane 2010).These changes suggest some gradual variation in the radial ionization and recombination evolution of the ejecta back in the corona.The charge states of C and O do not show the same shift to the highest charge states during the first half of the CME that Fe did.Additionally, panel (d) shows the C 6+ /C 5+ and O 7+ /O 6+ ratios, which dip during the hottest sections of the CME as will be discussed later.The second half of the ejecta is characterized by high charge states of C and O, with an intensifying contribution from the fully ionized charge states near the rear of the ejecta.
The elemental composition of the leading part of the CME structure is characterized by a significant enhancement of lowand mid-first ionization potential (FIP) elements (Mg, Si, S, Fe), with aFIP bias reaching above 10 in some places as shown in panel (e).These extremely high values of FIP enhancement are beyond what would be expected for slow wind but are characteristic of values observed in hot CMEs (Zurbuchen et al. 2016).The enhancement in low-and mid-FIP ions gradually declines toward the second half of the CME and relaxes to slow wind values toward the end.The elevated FIP bias (the ratio of each element to oxygen, compared to the photospheric ratio) indicates material that originated from an active region, where the most extreme FIP biases are observed (Feldman & Laming 2000;Widing & Feldman 2001).
Although this ejecta is made up of some of the most ionized material captured by SWICS, it is also observed to contain lower-ionization material than typically seen in CMEs and the slow wind.In particular, the sheath and the boundary between the flux rope and trailing material show ions of C <4+ , O 3,4+ , and Fe 5+ , potentially from prominence material (Lepri & Zurbuchen 2010;Lepri & Rivera 2021).This complex charge state composition provides a unique opportunity to constrain the heating in CME models and will serve as our reference values to compare simulations against throughout the rest of this paper.

MHD Simulation and Coupled NEI Modeling
The Bastille Day CME analyzed in the present study was previously simulated using the coronal solutions from the thermodynamic MHD code "Magnetohydrodynamic Algorithm outside a Sphere" (MAS; Török et al. 2018).In the Török et al. (2018) study, the CME was initiated from a pre-eruptive flux rope configuration in magnetic equilibrium that propagated through the corona and coupled to a heliospheric MAS framework to simulate its interplanetary propagation (Lionello et al. 2013).We note that the initial study did not include the nonequilibrium ionization modeling that is included in the present study.Török et al. (2018) describe the simulation setup, including the thermodynamic MHD relaxation of the global background corona; the construction, insertion, and relaxation of the pre-eruptive flux rope; and the CME's propagation through the heliospheric domain.
As presented in Török et al. (2018), a comparison between the synthetic and observed white-light, EUV, and X-ray emission shows nice agreement in the generally coronal morphology and the eruption.In line with observations, the model/simulated flux rope accelerates quickly, reaching 2500 km s −1 , compressing the front and forming a shock at the leading edge in the low corona (∼1.5 R e ).There is also a termination shock below owing to reconnection outflows.The shock locations are the sites of significant compressional heating, as will be discussed in the present study.The flux rope also shows a cool, dense prominence feature that is ejected as part of the flux rope, forming a typical multipart CME structure often observed in SOHO/LASCO white-light images.Comparison of synthetic in situ observations from its interplanetary propagation to observations at ACE shows that the bulk speed and magnetic field components are qualitatively consistent for a period that has been time-shifted and a magnetic field with a small multiplication factor.
In the present study, the MHD coronal solutions from Török et al. (2018) are coupled with a nonequilibrium ionization (NEI) module to track the formation of individual ions within the evolving CME structures to learn about the energetics of the eruption.NEI effects are important to consider when investigating the ion evolution in dynamic events, like flares and CMEs, because the rapid changes to the temperature, density, and speed can cause ionization and recombination processes to become unbalanced quickly.If NEI effects are important, the ion abundances will not reflect the local properties of the plasma, which has important consequences for emission diagnostics at the Sun and for how heavy ions measured in situ are interpreted (Rivera et al. 2019a(Rivera et al. , 2019b)).Details of the calculations of the NEI code are discussed in Shen et al. (2015).Following Shen et al. (2015), the set of equations governing the relative abundances is as follows: ( ) where F z i is the relative abundance of ions in an i + ionized state of an element with atomic number z. U is the plasma bulk speed, and n e is the electron density.C and R are the ionization and recombination rate coefficients, respectively, for each ion.Ionization and recombination processes include collisional and autoionization rate coefficients and radiative and dielectronic recombination.We note that the MHD equations are for a single fluid that assumes thermal equilibrium, T p = T e .Furthermore, a Maxwell-Boltzmann velocity distribution for the electrons is assumed.Also,charge states are ionized only a single ionization stage at a time, i.e., C 4+ to C 5+ .The ions are assumed to all flow at the same speed, i.e., no differential streaming.Initially, ionization equilibrium is assumed at the start of the simulation.

Analysis
We compare C, O, and Fe charge state distributions during the Bastille Day CME from ACE/SWICS to the simulated freeze-in charge states based on the model of the numerical simulation of Török et al. (2018).We trace the outflow of different structures within the coronal solutions, between the surface and 20 R e , by following the trajectories of a large number of passive tracer particles advected with the flow within the MHD calculation at the code time step.Then, the MHD quantities and charge state values are interpolated to the tracer positions from the 3D dumps of the calculation.The individual trajectories enable the investigation of charge state evolution within selected plasma structures identified through specific conditions.As such, we survey the parameter space to connect the simulated and observed charge state distributions.Specifically, we are interested in understanding the conditions in which the charge states form and where in the CME they originate, and we quantify the heating and cooling experienced, which ultimately form the freeze-in distributions observed in the heliosphere.
To compare simulations to the general structure observed in situ, we break the in situ observations into two main structures.As shown in Figure 1, the CME substructure is strongly delineated between an Fe charge state distribution peaked at Fe 16+ and one with a bimodal profile, suggesting a common thermal plasma history within plasma in these two structures.Therefore, we naturally divide the structure by the common features of the Fe distributions.We compare the full simulated and observed charge state distributions, in particular, we prioritize reconstructing the full Fe charge state distributions.The basis for starting with Fe is that its distribution will have the most sensitive thermal response, as compared to C and O, to plasma conditions given the number of ions (26).Also, through the extended range in the freeze-in altitudes of Fe charge states, the Fe ions would be the most strongly shaped by the evolution within the 20 R e coronal domain of the simulation.
As such, we designate the first part of the CME within the main flux rope structure, DOY ∼197.8-198.2,where the highest ionization stages of Fe are found.The average values of the C, O, and Fe charge state distributions associated with the first half of the CMEare shown in the top row of Figure 2.This period is characterized by some of the highest average Fe charge states observed in situ, 〈Q Fe 〉 = 16.8 (Rivera et al. 2022).However, this is not reflected in the O 7+ /O 6+ = 0.86 and C 6+ /C 5+ = 0.31, which show mundane slow speed wind values (Richardson & Cane 2004).This discrepancy will be further discussed in Section 5.The bottom row of Figure 2 shows distributions from the second half of the CME, DOY ∼198.8-199.2,that form the double-peaked Fe distributions.The plasma shows an overall 〈Q Fe 〉 = 11.8 and higher O 7+ /O 6+ = 1.24 and C 6+ /C 5+ = 4.3, in contrast to the first half the CME.We discuss the potential reason for this in Section 5.

Charge State Evolution
To match the charge states within the two main structures in the observations, indicated in Figure 2, we examine the simulated tracks capturing the thermodynamic and associated ion evolution.We search for tracks with specific freeze-in ions that escape to the edge of the simulation box at 20 R e , to identify material that would make it into interplanetary space and would be conceivably measured by in situ instruments.Because we are not following the CME structures in the heliosphere, we did not specify that the particles were exactly directed toward Earth given that we find many tracks within the same structure with solutions having a similar ion evolution.Instead, we aim to understand the conditions for the ion evolution that best represented the observations.
To capture the charge states within the first half of the CME, DOY ∼197.8-198.2,we searched for the highest freeze-in ions produced by the simulation.Figure 3 shows an example of a representative tracer particle that tracks the formation of material that produces an Fe distribution that freezes in with a 〈Q Fe 〉 ∼ 13.We note that the jagged appearances of the temperature and density profiles are due to the temporal resolution used.The MHD calculation uses much shorter time steps; however, the simulation snapshots presented in the plot are written out at larger intervals, causing some jumps in the variables.The written-out timescale will not affect the NEI modeling solutions given that the full resolution is used for that calculation.This material reflects the most highly ionized outflowing material produced in the simulation.From top to bottom, the figure shows the heliocentric evolution of the C, O, and Fe relative abundances; proton temperature; proton density; and bulk velocity of the plasma.The bottom row shows the final charge state distributions of C, O, and Fe at the edge of the simulation box, with the black curve representing the observed charge state distribution (top row from Figure 2).The material originates as part of the coronal material embedded in the flux rope as indicated by its initial conditions of approximately T ∼ 3 MK and N ∼ 10 8 cm −3 , which can produce Fe ions up to Fe 16+ but not to higher ionization states.The material is ejected and quickly accelerated to over 2000 km s −1 and heated to over 3 MK.The freeze-in ions reflect ionization states higher than observed for typical ambient solar wind, 〈Q Fe 〉 ∼ 13, C 6+ /C 5+ = 11.23, and O 7+ /O 6+ = 1.772 (Richardson & Cane 2004).However, the highest ionization states of Fe in the outflowing material still did not reach values observed in situ.In general, Fe >16+ is not formed in any significant way throughout the ejected material to match the Fe distributions measured within the magnetic cloud by ACE/ SWICS.
To match with the second half of the CME, DOY ∼198.8-199.2,we find that prominence material that is rapidly heated up reproduces the pronounced bimodal feature observed in situ well.The prominence material is identified by setting the initial conditions of the tracer particle as cooler (T 10 5 K) and denser (N 10 9 cm −3 ) than the surrounding corona, typical of prominence/filamentary material at the Sun (Parenti 2014).Figure 4 presents the radial evolution of a representative tracer particle showing cool, dense filament material rapidly accelerated to over 2000 km s −1 and heated to above 1 MK, below 2 R e .The heating is also reflected by the rapid ionization of C, O, and Fe.The temperature profile quickly cools as the plasma expands out to 20 R e , causing

CME Heating
Figure 5 shows three stages of the CME evolution just after liftoff that includes the two tracks associated with Figure 3, hot flux rope material (black points), and Figure 4, prominence material (gray points).At each time step (time = 2, 4, 6), the magnetic field line associated with the point along the tracks at that time step is included to indicate the location of the plasma.The top row shows isosurfaces at two shock heating values, 0.01 and 0.03 MK s −1 , along with two temperatures, 10 and 8 MK, and the fractional abundance of Fe 16+ .The bottom row shows the same configuration except with isosurfaces at three 〈Q Fe 〉 values, 13.5, 14.5, and 15.5.We note that the funnel shape in time = 2 and 4 in the bottom row is an artifact of the simulation related to a heat conduction front on open field lines prior to CME initiation.However, its appearance does not affect the material traced by the two tracks since the front is cleared prior to the tracer particle material interacting with it.
As discussed in Török et al. (2018), nearly 28% of magnetic free energy is converted into kinetic energy, causing rapid acceleration just after the onset of the eruption.The rapid liftoff generated a shock in the low corona with strong compression, as shown through the shock heating isosurfaces that surround the shock in the top row of Figure 5.The shock heating initially generated temperatures of 8-10 MK, the highest reached in the simulation, at a small region of the leading edge that quickly cooled farther out.However, the hot pocket of plasma was not part of the original flux rope material but rather surrounding coronal plasma that was cleared out by the eruption, which was not followed by the tracer particles.As seen in the subsequent steps (4 and 6), strong shock heating persists at the front and CME itself, as well as near the Sun (Reeves et al. 2019).Plasma at the termination shock is also subject to high ohmic heating owing to the current sheet that forms at the reconnection site.Moreover, in the bottom row and as also indicated by the Fe 16+ isosurface in the top left panel, the front and termination shock also coincide with regions where highly ionized Fe is produced owing to the strong shock heating experienced.Initially, a large fraction of the ejecta contains 〈Q Fe 〉 at 13.5, 14.5, and 15.5.However, as the plasma begins to expand and cool, ions recombine, acting to lower the overall 〈Q Fe 〉 of the outflowing plasma.While the reconnection site maintains the highest-〈Q Fe 〉 plasma.
To examine the distribution of the Fe ions in more detail, we show the distribution of the 〈Q Fe 〉 across a plane intercepting the polarity inversion line at the active region in Figure 6.The wedge illustrates regions of high versus low 〈Q Fe 〉 zoomed in to 6-16+ across the flux rope evolution.The panels also include the addition of two field lines, in light blue, that capture the main flux rope structure.The first step (left panel) shows a clear partitioning of the highest Fe charges states at the compressed front (via shock heating) and at the back of the flux rope associated with ohmic heating via the current sheet that forms.In the following time steps (middle and right panels), we can see the expansion and heating of the flux rope as parts within the structure become elevated in 〈Q Fe 〉 as they become more orange/yellow.It is also seen that significant heating continues to maintain high-〈Q Fe 〉 plasma localized to the reconnection site at the Sun, behind the flux rope.
Moreover, the heating terms at the two tracer particles are presented in Figure 7, showing the shock, conduction, radiation, and ohmic heating in MK s −1 , where prominence material is in black and the hot flux rope material is in gray, corresponding to Figure 5.It is found that heating via conduction, along with shock cooling, is important in both cases.The largest contribution to the heating within the prominence and hot flux rope material occurs through strong conduction within a few solar radii from the eruption site between neighboring material.Conduction is field aligned; therefore, hotter material at the shock front, where strong shock heating is taking place, transfers heat to the location of the tracer particle.Conversely, the negative shock term indicates strong cooling at both structures due to the rapid expansion.The expansion acts to counteract the heating taking place.Ohmic and radiation terms remain relatively small for both cases.

Freeze-in Distances
Figure 8 shows the freeze-in distances of selected ions, C 6+ , O 7+ , and Fe 16+ , within the evolution of the prominence structure in the simulation box extending to 20 R e .The freezein ion distance in the plasma is the location where ions no longer ionize or recombine, becoming fixed after that altitude.The plots follow the tracks of tracer particles with initial conditions like prominence material that originate within the flux rope, with a similar evolution to the example shown in Figure 4.The tracks are colored with the relative abundances of each ion, where the track ends at the point where freeze-in occurs.The gray lines are the projected paths on the r-f plane showing the heliocentric freeze-in distances for each path.In line with previous NEI modeling of CMEs, we generally find that the freeze-in distances are organized as C < O < Fe (Rivera et al. 2019b).We also show that the freeze-in distances are farther from the Sun (5-20 R e ) than is typically expected for solar wind conditions that span 1-5 R e (Landi et al. 2012;Boe et al. 2018), in line with previous work (Rivera et al. 2019b).The large freeze-in distances are a result of the high CME speed, which can transport plasma into the heliosphere quickly.

Discussion
Using a 3D MHD model of the Bastille Day CME event, we trace the C, O, and Fe charge state distribution evolution within its multipart structure using a coupled NEI code.We find that the simulated prominence and surrounding ejected material are significantly heated and accelerated just after liftoff, in line with observations.The prominence and flux rope structure reach a velocity well above 2000 km s −1 , and they are heated to several megakelvins.

Formation of Fe Charge State Distributions
In particular, the thermodynamic history of the main prominence core material produces an Fe distribution that matches the bimodal Fe signature observed within a significant portion of the trailing in situ CME structure.As shown in Figure 4, the initial conditions of the tracer particle indicate typical prominence properties at the Sun, below 1 × 10 5 K and above 10 10 cm −3 .The prominence is accelerated to 2000 km s −1 and reaches T = 3 MK below 2 R e and then rapidly cools, expands, and decelerates.This temperature, density, and speed profile produced significant ionization in the plasma early on and allowed for significant recombination to take place, most notably between 3 and 7 R e in the Fe charge states.Through this thermodynamic evolution, the doublepeaked (bimodal) Fe charge state distribution forms, producing a simulated 〈Q Fe 〉 = 10.366 as compared to the observed 〈Q Fe 〉 = 11.8.However, it failed to reproduce the general profile of the C and O ions, as shown with the black curves in the bottom row of Figure 4.
In both simulated cases, hot flux rope and prominence, we observe that the Fe is ionized up to Fe 16+ but does not ionize significantly after this point.This is due to Fe 16+ belonging to the Ne-like isoelectronic sequence (containing 10 electrons, with two filled electron shells, similar to Ne) that has an ionization potential much higher than the lower ions and therefore requires a significantly higher temperature to produce Fe 17+ ions and above (with a relative abundance peak around a  formation temperature of 6 MK) in the plasma (Del Zanna et al. 2021).Its higher potential is also reflected in the significantly lower total ionization rate within the relevant temperature range of the simulation, between Fe 15+ and Fe 16+ , shown in the top panel of Figure 9, produced using the fiasco Python package (Barnes et al. 2023).As such, the temperature necessary to produce Fe >17+ in a significant manner is not reached within the outflowing structures in the simulation.Subsequently, as the plasma cools, we observe that the Fe ions begin to recombine into lower ionization states.The total recombination rate of Fe 16+ to Fe 15+ decreases with decreasing electron temperature within the temperature range prior to ion freeze-in (1-4 MK), reaching values that are lower than those of neighboring ions, as shown in the bottom panel of Figure 9.The lower recombination rate results in the Fe ions in a lower ionization stage, immediately to the left of the Fe 16+ , recombining faster than Fe ions >16 + , while Fe 16+ lags behind and the plasma shifts to a lower 〈Q Fe 〉.
In the hot flux rope case, the plasma maintains a relatively high temperature, above 1 MK out to 8-9 R e , that maintains high enough ionization rates to combat significant recombination in the plasma prior to freeze-in.The result is the preservation of relatively high ionization states at freeze-in, 〈Q Fe 〉 = 13.149.Conversely, the prominence cools faster, having less time to ionize.Instead, the plasma rapidly recombines, as the recombination rate remains comparatively high.This process results in a larger shift from high to low ionized states with a peak remaining at Fe 16+ that lags behind.The second peak to the left of Fe 16+ is produced through the migration of ions to lower ionization states wherein the bimodal feature is formed.We find that the non-Fe 16+ peak is largely dictated by the freeze-in distance of the material, where plasma that freezes in closer to the Sun forms the non-Fe 16+ peak at higher charge states, while those that continue to recombine longer shift the non-Fe 16+ peak farther to the left of the distribution.This result is in contrast to other studies that suggest that the bimodal feature in CMEs is formed through the mixture of hot and cold CME material measured together, although this may occur as well (Gruesbeck et al. 2011(Gruesbeck et al. , 2012;;Song et al. 2016).As such, future heavy ion measurements below 1 au from the Heavy Ion Sensor on Solar Orbiter (Owen et al. 2020) will supply high temporal resolution measurements that should provide further insight into the spatial structure.

Anomalous Charge States of C and O
Moreover, we find that the C and O charge state distributions were not well reproduced in the hot flux rope or prominence case.In the front side of the CME, the observations reflect less ionized C and O distributions than expected in comparison to the high 〈Q Fe 〉 found in that part of the ejecta (16.8), while the opposite is found to be the case for the trailing part of the CME.The simulation did not capture this.As shown in Figure 1(d), the first part of the CME shows highly uncorrelated O 7+ /O 6+ and C 6+ /C 5+ behavior for the period below DOY 198.5, indicating the appearance of the so-called Outlier wind (Kocher et al. 2017;Zhao et al. 2017;Rivera et al. 2021).The Outlier wind is characterized by anomalously depleted C 6+ ions and other bare ions (i.e., O 8+ , N 7+ ) that are fully stripped of electrons.The Outlier wind is found mainly in CMEs (within 72% of all CMEs) and slower speed wind.The associated characteristics include an enhancement of He/H abundances and proton density compared to the surrounding material.They are generally found during heliospheric current sheet crossings connected to changes in the magnetic field direction and electron strahl character.In Outlier wind, C 6+ /C 5+ is artificially reduced and would be interpreted as lower ionized plasma compared to ions of elements that are least affected, i.e., Fe, given that Fe 26+ is not abundant in the solar wind or CMEs.This would also affect O.The fact that the MHD + NEI model did not reproduce the Outlier effects suggests that the physics driving the Outlier wind is not present in the model.The unusual C and O charge state distributions could be a nonthermal effect (Raymond et al. 2022) since the ionization and recombination rates do not seem capable of explaining the Outliers given that atomic rates relevant to C 5+ and C 6+ are not very different from each other.It is more likely that the anomaly in the bare ions is related to the cyclotron frequency that they all share with He 2+ .

Simulated Hot Flux Rope Charge States
Overall, we find that the highest ionization states observed in situ were not well reproduced by the simulation.In particular, the highest freeze-in 〈Q Fe 〉 values reached approximately 13 in the simulation, while observed values show an average across the first section of the CME of 〈Q Fe 〉 = 16.8.Specifically, from Figure 1 it is observed that the flux rope in situ coincides with the Fe distribution peaked at Fe 16+ , as summarized in Figure 2. Therefore, from the spatial distribution of 〈Q Fe 〉 in Figure 6, additional heating is required within plasma surrounded by the light-blue magnetic field line that constitutes the flux rope in connection with the in situ material.
There are several potential reasons for this result.One that has already been alluded to is that the simulated temperature is too low to reproduce such highly ionized states.This would suggest that the heating experienced by the CME overall is insufficient or spread out too much in space or time.However, there could be other important factors at play.
One possibility could be the presence of nonthermal electrons produced at the flare site that generate additional ionization not captured in the model.Recent studies have indicated that a significant nonthermal electron population is produced by flares and could also be present in noneruptive sites like active regions (Del Zanna et al. 2022;Fleishman et al. 2022).It is well established that energetic electrons will arise from particle acceleration at the current sheet owing to contracting magnetic islands incited by magnetic reconnection such as the case of flares (Drake et al. 2006;Arnold et al. 2021).Nonthermal electrons would be particularly important for the Bastille Day CME, which produced a very energetic X5.7 flare.Higher electron densities at the high-energy tail of the distribution would generally be reflected in the plasma as an enhancement of more highly ionized charge states that would be produced through a Maxwell-Boltzmann profile.In particular, electron distributions represented by low kappa numbers, i.e., κ = 2, 3, 5 near the theoretical limit, would significantly change the ionization and recombination rates compared to a Maxwell-Boltzmann distribution (Dzifčáková & Dudík 2013).As shown in Figure 2 of Dzifčáková & Dudík (2013), Fe XVII, or Fe 16+ , shows that ionization rates would be notably larger, particularly for κ = 2 in the lower temperature range, within 1-4 MK, that is reached prior to freeze-in within the simulated hot flux rope and prominence plasma (see Figures 3 and 4, respectively).The larger ionization rates for Fe 16+ and more highly ionized Fe ions that would follow a similar ionization rate trend indicate that for such an electron distribution (low κ) and within the temperatures reached in the simulation, significantly more Fe >17+ would already be expected.This points to particle acceleration mechanisms, like Fermi acceleration, that form suprathermal tails in the electron distributions as being a critical aspect of extreme CME events like the Bastille Day eruption.
Another possibility would be additional ionization from photoionization through an enhancement of EUV and X-ray flare-related emission.The soft X-ray flux produced by the X5.7-class flare associated with the eruption, as taken from GOES, peaks at 5.7 × 10 −4 W m −2 in the 1-8 Å band, or 2.29 × 10 8 photons cm −2 s −1 near Earth.The photoionization cross section of Fe 16+ peaks at 3 × 10 −19 cm 2 (Reilman & Manson 1979).This indicates that the ionization timescale for the Fe 16 ions is on the order of a year or more.This suggests that photoionization might influence the gas below Fe 16+ but will not produce ions above Fe 16+ to change the results.
Another aspect that is not considered in the NEI modeling, but is likely important, is multiple ionization, which would affect Fe but not C or O. Electrons above 10 keV can remove n = 1 electrons from Fe.That leaves a hole in the K shell that will be filled when an electron from a higher n-level drops into it.That can happen with the emission of a Kα or Kβ photon, so that the result is an ordinary single ionization plus a photon.More often, the ion will get rid of the energy by ejecting another electron (Auger ionization), or even more than one.Thus, a single collision produces multiple stages of ionization.This is not usually too important, but either a κ electron distribution or localized pockets of very hot gas would increase their likelihood.It would be something that would discriminate between C and O, on the one hand, and Fe, on the other.

CME Heating
We find that the main driver of the heating in the ouflowing plasma is due to field-aligned conduction via plasma at the shock front heated through shock compression.As shown in Figure 5, locations of high compression show strong shock heating.The locations of the two tracer particles show high conduction heating during their early stages of propagation, in line with where shock heating is strong, as indicated in Figure 7.As discussed, the simulated temperature of the outflowing flux rope material did not reach values that would be able to produce the highest charge states in the observations connected to the passage of the flux rope.Therefore, although the kinematics, general morphology, and propagated structure of the CME match well with remote and in situ observations from Török et al. (2018), the extreme thermodynamic nature of the eruption was not reproduced.As outlined in the previous section, there could be several nonthermal effects in addition to the extreme heating necessary to produce Fe >16+ .The processes become important for modeling transient events, as compared to ambient solar wind conditions, and should be explored to build a more physically accurate and comprehensive model of eruptions.

Conclusions
In this study, we investigate the charge state evolution using an integrated nonequilibrium ionization module coupled to the MAS MHD simulation of the Bastille Day CME.We connect the heavy ion features measured in situ to the CME's thermal history in the corona through a comparison of the simulated and observed C, O, and Fe heavy ion signatures of the event during its passage at ACE.In situ observations of the Bastille Day showed that it contained highly ionized Fe, with ions ranging between 9 and 24+, with an average Fe charge state of >16+ within the flux rope.The trailing part of the CME showed that a large portion contained a bimodal Fe distribution, with a peak in Fe ranging between 8 and 11+ and the second peak at 16+.
From the simulation, it is found that the heating experienced within the prominence material and flux rope plasma traced in the simulation, which produced the closest matches to the in situ ions, is generated through field-aligned conduction via the plasma at the CME front.Given the large acceleration of the CME, reaching above 2000 km s −1 within a few solar radii, the leading edge experienced significant compression that led to strong and persistent shock heating there.
Results show that the bimodal feature of the Fe charge state distribution observed in the CME, and many CMEs in general, can be generated through rapid heating and cooling of prominence material.However, we note that the possibility of hot/cool plasma mixing and being measured together cannot be entirely ruled out for some cases.High temporal resolution measurements from the Heavy Ion Sensor on Solar Orbiter will provide better means of distinguishing between CME substructure in less evolved ejecta closer to the Sun.
Overall, the C and O charge state distributions were not well replicated.The CME contained abnormally low C 6+ /C 5+ compared to O 7+ /O 6+ connected to the passage of the flux rope in the front of the CME, where the highest Fe ions were measured.These ion dropouts are often observed in transients and during current sheet crossings that alter the C and O charge state distributions and elemental composition (Kocher et al. 2017;Zhao et al. 2017;Rivera et al. 2021).The dropouts abnormally deplete fully ionized ions, C 6+ and O 8+ , in the plasma, which is likely the reason for the C and O distributions not being well constrained.The anomalous depletion of the ions is thought to be produced by means of nonthermal effects that are not able to be replicated through the MHD model (Raymond et al. 2022).
Lastly, we found that the highest Fe charge states observed in situ were not reproduced by the simulated plasma.This result suggests that there may be some physics missing in the simulation, in particular, additional heating terms to reach the necessary temperature to form some of the highest ionized charge states, Fe >16+ , observed in situ.Other studies have also concluded that conduction heating alone is not enough to explain the thermal energy changes in erupting prominences (Landi et al. 2010;Murphy et al. 2011).One idea is the dissipation of small-scale magnetic stresses that are not captured in the initial conditions of the simulation because the magnetic boundary condition at the surface is too smooth.As discussed, the development of nonthermal electrons is likely a significant source of additional, unconsidered ionization to the plasma.This strongly suggests that mechanisms like Fermi acceleration driven by contracting magnetic islands at the current sheet are a critical source of particle acceleration during some of the most extreme eruptions, like the Bastille Day CME.

Figure 1 .
Figure 1.Stack plot showing the in situ observations for the Bastille Day CME, within the red vertical solid and dashed lines.Panels (a)-(c) show the relative charge state distributions for carbon, oxygen, and iron, with a red curve in panel (c) showing the 〈Q Fe 〉.Panel (d) shows the O 7+ /O 6+ and C 6+ /C 5+ ion ratios; panel (e) shows the elemental composition of Mg/O, Fe/O, Si/O, S/O, C/O, N/O, and Ne/O compared to its photospheric ratio; and panels (f)-(h) show the proton temperature, speed, and density, respectively.Lastly, panel (i) shows the magnetic field components and magnitude.

Figure 2 .
Figure 2. Mean charge state distributions of C, O, and Fe from ACE/SWICS for DOY.

Figure 3 .
Figure 3. Hot flux rope material tracer evolution matched to the first part of the CME.Top to bottom: C, O, and Fe fractional charge state distributions; proton and electron temperature and density; and bulk speed.The bottom row shows the final, freeze-in, simulated charge state distributions at 20 R e of C, O, and Fe (blue bars), along with the observed ions from ACE/SWICS (black curve).

Figure 4 .
Figure 4. Heated prominence material tracer evolution matched to the second part of the CME.Panels are the same as in Figure 3.

Figure 5 .
Figure 5. Three time stamps in the simulation showing (top row) two shock heating isosurfaces, along with Fe 16+ and two temperature isosurfaces.The bottom row shows three 〈Q Fe 〉 isosurfaces at the same time steps.The spheres show the location of the plasma tracked in the hot flux rope material (gray) and prominence structure (black) from Figures 3 and 4, respectively, along with the connected magnetic field line at the associated time step.

Figure 6 .
Figure 6.Three time stamps in the simulation showing the two plasma tracks in the hot flux rope material (gray) and prominence structure (black) from Figures 3 and 4, respectively.We include two additional field lines different from Figure 5, in light blue, that capture the main flux rope structure.The panels include the distribution of the 〈Q Fe 〉 across a plane intercepting the polarity inversion line at the active region.The wedge illustrates regions of 〈Q Fe 〉 between 6 and 16+ across the flux rope evolution.

Figure 7 .
Figure 7. Heating terms along the two tracer particle trajectories, where the black curve represents the prominence material and the gray curve represents the hot flux rope material in connection with the tracks in Figure 5.

Figure 8 .
Figure 8. Freeze-in distances of C 6+ , O 7+ , and Fe 16+ within the heated prominence tracks.The tracks flow within the simulation box and stop out at the point where each ion reaches its freeze-in value.The colors indicate the fractional abundance at each point.The gray curves are the tracks projected on the r-f plane.

Figure 9 .
Figure9.Total ionization and recombination rates with increasing electron temperature for Fe 12+ to Fe 19+ produced using the fiasco Python package(Barnes et al. 2023).