Tracking X-Ray Source Movement in a Retracting Flux Tube

Solar flares produce sources of localized, enhanced X-ray emission, thought to be due to the acceleration of nonthermal electrons and the transport of energy away from the reconnection site. The 2002 November 28 C1.6 limb flare showed clear X-ray source motion in the Reuven Ramaty High Energy Solar Spectroscopic Imager observations at 3–10 keV propagating from the apex of the flaring arcade, down toward the footpoints, and then rising back into the corona. Previous work attempted to model this motion using simulations driven by heating with an electron beam or thermal conduction front, finding reasonable agreement only if there were large initial densities. This work extends the previous model by considering a flux tube that retracts through a current sheet away from a magnetic reconnection site. The retraction model includes drag to slow motion in the current sheet, which allows us to vary the energy released by the retraction. This retraction causes a dense and superhot plug of material to form at the loop apex, naturally causing a thermal X-ray source to form in the corona. We find that the observed X-ray source motion, however, is most likely thermal and a signature of the evaporation fronts after initially filling the flux tube.


Introduction
Solar flares are one of the most energetic events in the solar system, thought to be triggered by magnetic reconnection in the solar corona that allows for built-up magnetic stress to be released. In the standard flare model, this reconnection occurs in a current sheet, and the reformed magnetic field is allowed to relax, converting the magnetic energy into kinetic and thermal energy (Carmichael 1964;Sturrock 1966;Hirayama 1974;Kopp & Pneuman 1976). This energy, often in the form of electrons accelerated to high energies (10 keV), is transported from the corona to the lower solar atmosphere, heating the plasma and raising the pressure, which causes upflows and downflows in the chromosphere, known as chromospheric evaporation and condensation (Neupert 1968). This phenomenon was observed in X-ray and radio wavelengths in Ohki (1975), and upflowing material was seen by the Bent Crystal Spectrometer on the Solar Maximum Mission in blueshifted components of Ca XIX and Fe XXV lines (Antonucci et al. 1982). For a review of chromospheric evaporation as observed in the extreme ultraviolet, see Milligan (2015).
Chromospheric evaporation can be driven by heating due to the collisional interaction of accelerated particles with ambient plasma. A significant portion of the energy in a flare is used to accelerate charged particles (both electrons and ions) in the solar atmosphere (Emslie et al. 2012;Benz 2017). Nonthermal radio and hard X-ray bremsstrahlung emission can be used to determine the accelerated electron population through an inversion calculation, but there are many uncertainties in the process (Holman et al. 2011;Benz 2017). These accelerated particles then interact collisionally with ambient plasma in the chromosphere to drive evaporation there.
Reconnection itself is likely to accelerate particles or create an environment that accelerates particles (Sweet 1969). The precise mechanism that accelerates particles is not known with certainty, but there are many possibilities. Some consider that the electric field in the current sheet can accelerate particles to the needed energies (Holman 1985;Litvinenko 1996), while others propose that there are shock-particle interactions that lead to the acceleration (Achterberg 1990;Somov & Kosugi 1997). Another option is the interactions between particles and plasmoids created by the tearing-mode instability (Drake et al. 2005;Karlický & Bárta 2011;Guidoni et al. 2016). Recent work has shown that mergers between successive plasmoids can replicate the observed nonthermal electron spectrum in flares (Arnold et al. 2021).
One-dimensional (1D) models have long been used to model the behavior of chromospheric evaporation. Typically they model energy transport from the corona to the chromosphere via a beam of electrons (MacNeice et al. 1984;Fisher et al. 1985a;Allred et al. 2005;Reep et al. 2020), through thermal conduction (Cheng et al. 1983;MacNeice 1986;Guidoni & Longcope 2011;Longcope & Klimchuk 2015), or Alfvénic heating through waves or turbulence (Kerr et al. 2016;Reep & Russell 2016;Ashfield & Longcope 2023). These works above approach modeling an electron beam by considering a generation site in the corona with particles streaming down to the chromosphere in the "thick-target" approach (Brown 1973). This contrasts with approaches that generate nonthermal electrons near or in the chromosphere to generate X-rays (Fletcher & Hudson 2008). Typically these simulations are done with a constant cross section and consider fluxes high enough to drive explosive evaporation (>5 × 10 10 erg cm −2 s −1 at 20 keV cutoff; Fisher et al. 1985aFisher et al. , 1985b. The threshold between gentle and explosive evaporation found by Fisher et al., however, depends not only on the energy flux but also on the low-energy cutoff and spectral index (i.e., the depth of heating in the chromosphere; Reep et al. 2015).
In this work, we focus on a C1.6 flare observed by the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI; Lin et al. 2002) in X-ray spectra and imaging on 2002 November 28 (Sui et al. 2006). This flare has been the subject of multiple papers tracking the motion of the X-ray sources (Sui et al. 2006;O'Flannagain et al. 2013). Those works concluded that this flare was part of a minority of flares that have early and strong nonthermal emission (Sui et al. 2007). This flare occurred close to the limb, with observed emission ranging from 3 to 50 keV (Sui et al. 2006). The footpoints were unocculted, but the observation sees a flare loop face on. High-energy observations, from ≈20 to 50 keV show the footpoint emission from the presumed electron beam. The lower energy from 3 to 10 keV shows a loop-top source, which splits and descends to the footpoints, then returns to the apex of the loop. This lower-energy emission is typically attributed to thermal emission from heated material. The RHESSI contours can be seen in Figure 1(c) of O'Flannagain et al. (2013), which finds that the 3-7 keV components are dominated by thermal emission, but above 7 keV show evidence of nonthermal emission at the early phases of the flare. We have reproduced measurements of the height of the X-ray sources above the footpoint in Figure 1  This flare has been previously modeled by . Those authors considered a variety of in situ coronal and nonthermal electron-heating parameters to best model the observed source movement. One of the runs from that work is reproduced in Figure 1 as the set of dashed lines, originally from Figure 1 of . That work found that while both types of heating could produce source movement, the specifics were best replicated by the nonthermal electron beam, but reasonable agreement was found only under the assumption that there was significant pre-flare heating producing a coronal electron density above 10 10 cm −3 .
In this work, we consider the possibility that the source motion is not solely due to in situ coronal heating or nonthermal heating near the footpoints. Instead we consider a scenario where there is self-consistent coronal heating due to the retraction of the loop alongside nonthermal electron heating of the same flux tube. We vary the partition of energy between these two effects to consider how they would impact the evolution of X-ray sources like those observed by RHESSI.
In Section 2, we discuss the numerical modeling and an example simulation. In Section 3, we explain the process of calculating and fitting the synthetic X-ray emission. In Section 4, we present the parameter space covered by the simulations and contrast them with observations. Finally, Section 5 discusses and presents the conclusions of this work.

Methods
The thin flux tube (TFT) framework is a model considering the behavior of an isolated bundle of flux against a static background. Originally developed for use in the solar interior by Spruit (1981), it has since been adapted for use in the solar corona by Guidoni & Longcope (2010). This framework allows us to keep the simplicity of a 1D model while also allowing this flux tube to move in space. The advantage to this approach is that the flux tube evolves under external forces not present in more traditional 1D models.
This work uses the Post-Reconnection Evolution of a Flux Tube (PREFT) code (Longcope & Klimchuk 2015). PREFT is a Lagrangian based code, advancing the following equations:   Figure 1 of . In both cases, red is the 3-6 keV band, light green is the 6-8 keV band, and navy is the 8-10 keV band. Earlier attempts to mimic the observed motion fail to capture the gradual movement.
In the above, D Dt is the advective derivative, ρ is the fluid density, B is the external magnetic field, p is the gas pressure, v is the fluid velocity, l is the length along the flux tube, b is the direction tangent to the flux tube, and h is the parallel component of the dynamic viscosity. From Equation (2), σ is the ad hoc drag coefficient, and v ⊥ is the component of velocity perpendicular to the flux tube. In Equation (3), x is the position vector composed of x, y, and z; r is the perpendicular direction, obtained by r y b =´ˆwith ŷ being the direction out of the page in the top middle panel of Figure 2; k is the spring constant; and γ is the damping coefficient. In the energy equation, Equation (4), c ṽ is the volumetric heat capacity, T is the fluid temperature, m is the mean mass of the fluid particles, k B is Boltzmann's constant, κ is the thermal conductivity, n e is the electron density, Γ is the radiative loss function, and Q is additional heating supplied to the flux tube. In Equation (5) q e is the charge of the electron, Λ is the Coulomb logarithm, δ is the power-law index of the electron distribution, μ 0 is the initial pitch angle of the electron beam, B x is the incomplete Beta function, F(t) is the amplitude of the injected distribution of electrons, E is the electron energy, E c is the low-energy cutoff of the distribution, N is the electron column along the loop, and N c is the stopping column at the low-energy cutoff. Equations (1) and (4) are the momentum and energy equations that are advanced through time. Mass conservation is handled through the Lagrangian framework with each cell given a mass element and the cell moving instead of the fluid. Equations (2), (3), (5), and (6) are all additional components that are supplied to the momentum and energy equations. Each of these equations will be discussed in more depth in the following paragraphs.
In the Lagrangian framework, the grid moves under the influence of the momentum equation as indicated by the use of the advective derivative in Equation (1). The right-hand side of the equation represents the sources of the motion. The first term is the total tension acting against the flux tube. The second term is the gas pressure along the tube. The third term is the viscous damping of fluid motions. The fourth term is the magnetic pressure exerted on the tube from the gradient of the external field. The final term is for any external ad hoc forces imposed on the flux tube, two of which are used in this work and discussed below.
The external forcing imposed on the flux tube are proxies for two effects. The first is Equation (2), which is a basic model of drag. This model of drag is sourced by the ambient plasma in the path of the flux tube's retraction and is used to siphon off kinetic energy from the flux tube as it retracts. A more in-depth discussion of this force and its use in PREFT can be found in Unverferth & Longcope (2021). The second force is the damped spring as shown in Equation (3). This force represents the pressure of underlying flux that would halt retraction. This existing flux either previously reconnected or existed below the flare arcade prior to the flare. Here we model it as an  Table 1. Top row: the initial magnetic and spatial configuration of the flux tube: position in the middle, the field as a function of length and height along the current sheet on the left and right, respectively. Bottom row: the initial equilibrium set in the flux tube: density on the left, pressure in the middle, and temperature on the right. The vertical dashed line marks the midpoint of the flux tube in each panel. This set of initial conditions apply to runs 1-5 and 10-25. Runs 6-9 use a similar setup with the retracted length and same coronal mass. overdamped spring to reduce oscillatory behavior. We work to minimize the oscillatory behavior for two reasons. The first is to reduce any unwanted impact of the spring force on the parallel dynamics while effectively stopping the vertical movement. Second is that the type of oscillations that this would lead to have periods similar to or longer than the duration of these simulations (Aschwanden et al. 1999;Wang & Solanki 2004;Ofman & Wang 2008). Oscillations at the base of the current sheet in PREFT have periods on the order of seconds to tens of seconds. The observed transverse oscillations in the works cited have periods on the order of minutes.
PREFT advances a temperature based form of the energy equation. The right-hand side of Equation (4) contains the source terms for temperature evolution. The first term is the adiabatic work from the movement of grid points. The second term is the viscous heating that results from the viscous interaction in the momentum equation. The third term is thermal conduction along the flux tube, using a modified form of Spitzer conductivity, which smoothly transitions to the freestreaming limit (Longcope & Klimchuk 2015). The fourth term is the radiative losses of the fluid, a tabulated fit to radiative losses calculated from CHIANTI 10.0 (Dere et al. 1997;Del Zanna et al. 2021). The final term is the heating added to the loop, which is split into two components. The first component of the heating is the equilibrium heating, which maintains the initial condition in the absence of other dynamics.
The second part of the heating added to the flux tube is the electron-beam heating. The form of this heating is shown in Equation (5), and the form of the injected source distribution is shown in Equation (6). We use a heating form adapted from Emslie (1978). This formulation is an approximation to the full Fokker-Planck solution and works well in the range of lowenergy cutoffs and simplified atmosphere we are using (Allred et al. 2015). We do not include a return current in this approach, which has been shown to alter the heating rate from a given electron flux . The beta function used in this equation transitions from the incomplete to complete form depending on the ratio of the current column to the stopping column. In addition, we have chosen to consider the fully ionized case. We do this as PREFT is a single fluid code and does not track the ionization of the fluid. In addition, the minimum initial temperature is 30,000 K, at which point hydrogen is more than 99% ionized. At this temperature, helium is roughly 55% in helium II and the remainder in helium I; for the temperatures of the transition region and above, it is mostly ionized. The electron heating is sensitive to the parameters used to describe the injected electrons. Equation (6) is a typical power-law injection with a lowenergy cutoff. Above this cutoff, the form is as described, while below this cutoff, the amplitude is zero.
The final step to begin a PREFT simulation is setting the initial and boundary conditions. The initial conditions initialize the flux tube in an Rosner, Tucker, Vaiana (RTV)-type equilibrium (Rosner et al. 1978). We construct this state with a defined length and minimum and maximum temperatures. We add on an equal length of chromosphere to either end of the flux tube. This chromosphere is isothermal and gravitationally stratified, serving only as a mass reservoir. These conditions can be seen in the lower row of panels in Figure 2. The pressure is uniform in the corona and then increases following the density increase in the chromosphere. The ends of the chromosphere are fixed in space and closed, sealing the tube from mass and energy loss through the ends.
With an initialized equilibrium flux tube, we place the tube in a Green-Syrovatskiǐ (Green 1965;Syrovatskiǐ 1971) current sheet. The tube is placed in the current sheet and folded over at a specified height and angle to create a flux tube just after reconnecting. From this configuration, the flux tube is ready to evolve according to Equations (1) and (4). The flux tube is shown in the current sheet in the upper middle panel of Figure  2. The left top panel is the magnetic field strength as a function of the length along the current sheet, and the right top panel is the strength as a function of height where 0 is the edge of the current sheet. From this, the field strength peaks in the middle of the flux tube as it nears the midpoint of the current sheet. Below z = 0 is the region where the flux tube will be affected by damped spring force.

Examining the Simulation in Depth
To explain many of the features we reference in later sections (see Sections 3 and 4) of this paper, we consider the evolution of a single simulation in depth. This simulation is energized by a combination of field-line retraction and an electron beam. Simulations mentioned later in this paper have the same energy input but split into different proportions to explore the effects of the partitioning. This run is calibrated to have an even split between energy in the electron beam and that provided by retracting through the current sheet. The flux tube shortens from 43.9 to 36.4 Mm, and the electron beam peaks at 4.18 × 10 9 erg s −1 cm −2 . The beam heating acts for the first 20 s of the simulation with a triangular time profile peaking at 10 s with a low-energy cutoff of 10 keV and a spectral index of 5. In addition the drag parameter, σ is set to 0.725 Mm −1 , chosen to remove 50% of the energy released by the retraction, such that 50% of the energy input is from the beam and retraction each. This simulation is run 21 in Table 1.
The first row of Figure 3 shows the initial conditions of the simulation. The flux tube is initialized into the RTV equilibrium mentioned previously. In this case, the coronal pressure is scaled to have a consistent coronal mass with the simulations that do not retract. The maximum temperature is 3 MK at the apex, and the minimum is 30,000 K at the chromosphere. The radiative loss function is set to zero below a value of 33,000 K, such that the chromosphere does not require constant heating. The field at the apex of the flux tube is 135 G, and the edge of the current sheet is 71.5 G. This flux tube reconnects with an angle of 116°, where 180°is antiparallel reconnection without any guide field. The field strength and excess length control the amount of energy released, while the angle influences the rate of energy conversion and bulk motions during retraction as described below.
The second row of Figure 3 shows the simulation after 0.15 s. The flux tube's retraction has begun, and drag has only had a minimal effect. The flux tube begins to split the initial bend in two as the apex retracts downward. As the bend passes along the tube, the tube is accelerated downward and inward at the local Alfvén speed. As a result, the jets of material inward toward the apex resemble that of a simulation with no drag, with near-sonic flows of material moving inward. These jets of material interact viscously where they meet at the apex. This collision results in a rise in the density as the material is compressed. The collision is also significantly viscous, and as such much of the energy in the flows goes into heating the plasma there. This is visible in the temperature, density, and pressure plots in row 2 as the spike occurs at the apex.
We show the simulation after 1 s in the third row of Figure 3. The flux tube has continued to retract, and drag has slowed this substantially. As drag acts on the retraction, what would be two clear bends without drag is smoothed out to a much gentler curve. The jets are reduced in strength due to the drag acting on the bends. This reduced strength means that the initial spike in pressure from the plug of hot material stalls the accumulation of material in the apex. After stalling the flows, this plug spreads out under its high pressure. As the dense hot plug of material is no longer growing through the jets, it has peaked in temperature at 55.7 MK. This spread is typically contained by the movement of additional material to the edges of the plug. With weaker jets, this effect is not sufficient to contain the initial plug.
The fourth row of Figure 3 shows the simulation at 5 s, when retraction has just finished. The spring force has halted the motion of the flux tube at z = 0. This force minimizes impact on the flux tube dynamics by working only on the perpendicular velocity but does generate some additional viscous heating as the force is applied, resulting in a slight rise of coronal temperatures. The jets now subside, and no more material is added to the apex. The plug now expands faster without any confinement from inflowing material. During retraction, as the flux tube moves from high field strength to low, the flux tube expands. This expansion reduces the density seen in parts of the flux tube that retracted from the current sheet. The density at 5 s (black solid line) has dropped below the initial condition (solid light gray line), excluding the plug of material at the apex, erasing the enhancement visible previously. This occurs as PREFT is in a constant flux framework; as the magnetic field drops, the cross-sectional area of the flux tube expands. For a more in-depth discussion of this behavior, see Unverferth & Longcope (2020).
The fifth row of Figure 3 shows the flux tube at 7 s, with signs of chromospheric evaporation beginning to drive material into the corona. The electron beam is now at an amplitude where it has a visible effect in the chromosphere. Additional heat has reached the site of the evaporation and condensation via the thermal conduction front. Both sources of heating have raised the pressure there, and now evaporation and condensation can be observed in the velocity as the diverging flow at −14 Mm. The peak in density is the compression at the leading edge of the downflowing condensing material. The narrow enhancement of pressure is also evidence of evaporation as this gradient drives out material in both directions.
The sixth row of Figure 3 shows the simulation at 20 s. The electron beam has finished, marking the end of energy input to the simulation as retraction ended 15 s prior. The evaporation is well developed at this point, reaching approximately halfway up the flux tube. The peak of the evaporation upflow is supersonic. The sound speed is shown as the dashed red line in the velocity panel. The front is enhancing the coronal density of the flux tube by a factor of 3 over the value before the front, in addition to the enhancement from the initial length contraction.
The bottom row of Figure 3 shows the simulation at 29 s. This time step contains the final event of note. The evaporation fronts from each leg of the tube have met at the apex of the flux tube. This collision produces another central spike of density and temperature enhancement though it is minor and persists for only a few seconds. The two fronts will pass through the other and continue to the opposite footpoint. The evaporation fronts will dissipate as they travel, spreading out and reducing the peak velocity of the flow. As they bounce, the loop slowly cools through radiative losses, and evaporation adds material to the coronal segment of the loop at a reduced rate compared to the initial fronts. The evaporation will continue until the pressure spike near −15 Mm decays to meet the rising coronal pressure, which does not happen during this simulation. The simulation is terminated at 150 s to match the duration of the observations in Figure 1.

Data Extraction and Synthetic Observables
We take the output from the PREFT simulation and use it to generate synthetic X-ray spectra. There are two primary components of this, both of which are calculated after the simulation is complete. The X-ray emission can be broken down into the thermal emission of the plasma and the nonthermal bremsstrahlung emission arising from the beam of injected electrons.
We calculate the nonthermal emission of the flux tube for the time that the electron beam is active. We numerically evaluate the bremsstrahlung integral using the relativistic cross section from Bethe & Heitler (1934). This approach is similar to the approach used in  to evaluate the nonthermal emission from HYDRAD 3 (Bradshaw & Mason 2003;Bradshaw & Cargill 2013) simulations.
The thermal X-ray emission is calculated for the entire duration of the simulation. At each time step we use CHIANTI Note. In each run, the electron beam lasted 20 s, with a 10 s rise and fall, and also had a power-law index of δ = 5. 10.0.1 (Dere et al. 1997;Del Zanna et al. 2021) to find this emission. We take into account the free-free, free-bound, twophoton continuum, and line emission. These routines were run using the photospheric abundance of Asplund et al. (2009). We use the photospheric abundance as flare emissions more closely resemble photospheric than coronal abundances in most events (Warren 2014) although the abundances have been measured to vary with time (Mondal et al. 2021). We focus on the synthetic emission from 3 to 10 keV, which matches the range covered by the observations. At each time step, we synthesize the emission of the flux tube. When we sum the emission over the whole loop, we see an energy distribution like the top panel in Figure 4. Thermal emission dominates below 10 keV and nonthermal emission above 10 keV. RHESSI observations of the flare show a significant amount of nonthermal emission during the peak intensity in X-rays (O'Flannagain et al. 2013). This simulation shows similar behavior in that the emission below 10 keV is primarily thermal, with nonthermal emission dominating above that energy. Instead of prescribing an area at the start of the simulation, PREFT is run normalized. To fix the size of the flux tube, we scaled the PREFT simulation to have a magnetic flux of 7 × 10 19 Mx. This value was chosen such that the synthetic emission in the GOES bands matched the GOES light curves of the flare.
With synthetic emissions calculated, the simulation is then put through a fitting procedure to localize the X-ray emission to a source location. For each time step, the emission in each of the bands of interest (3-6, 6-8, and 8-10 keV) is integrated over energy to produce a single emission value along the length of the flux tube. An example of this is the black curve in the lower panel of Figure 4, which is the 6-8 keV band at 10 s. This emission is then taken from the irregularly spaced PREFT grid and placed in 200 regularly spaced bins roughly 180 km each. This is higher resolution than the roughly 2″ of resolution for RHESSI. Reducing the number of bins to match RHESSI's resolution did not significantly alter the result of fitting. To match the observational cadence, instead of fitting each time step as shown in Figures 4 and 5, we integrate emission over 8 s with a 4 s cadence, e.g., 0-8 s, 4-12 s, 8-16 s. This combined emission is then masked in two steps. The first is to consider only half of the loop; in this case we consider from −18 to 0 Mm along the flux tube. The second is to consider only emission greater than 30% of the maximum emission, as was done in O' Flannagain et al. (2013). We use both of these to ensure comparability to the results of O'Flannagain et al.
(2013) and . We then calculate the centroid of the emission in the same method as  and derive a width according to the following: With x being the coordinate along the flux tube, I j is the emission of a given energy band in the jth cell, j the bin number summed over, |x| the centroid location, and w the width of the distribution around the centroid. This allows us to accurately calculate a location even where the shape of the emission would deviate from function fitting of a Gaussian or similar function.
The results of this masking can be seen in the lower panel of Figure 4. The initial curve is shown in the solid black line. The With a source location in each band, we consider the motion of these sources for the lifetime of the simulation. The emission for each band can be seen in the left panel of Figure 5 taken at the interval between 4 and 12 s. Note that the coordinate used here has changed, so 0 is now the left footpoint of the flux tube instead of the midpoint. In O'Flannagain et al. (2013) the height coordinate was length along a reconstructed loop shape anchored to the high-energy footpoint. In PREFT, the height of the chromosphere is dynamic, responding to energy deposition throughout the simulation. We place the 0 point at the end point of the flux tube instead of attempting to have a moving 0 point associated with the top of the chromosphere. While all three bands have a peak occurring close to each other, the centroid of the 3-6 band is near that peak. The 6-8 and 8-10 keV bands instead have significant emission outside their peak, and as such, they show a much broader distribution, and the centroid is placed far from their maxima with a much larger width. Placing these together for each time step, we generate the right panel of Figure 5. We can see that the three bands move similarly. After 20 s there is no nonthermal emission, and so we would expect all sources to move uniformly in the corona, tracking with temperature and density in a similar way.
The X-ray sources move following the events of the simulation. For the initial few seconds, the apex is the brightest feature due to the dense plug of >30 MK plasma. This material is hot and dense compared to the surrounding plasma, but that emission is dwarfed by the much denser chromosphere as the beam heating and thermal conduction front reach it. At 12 s the emission has moved from the apex down to near the chromosphere, as seen in panel (a) of Figure 6 as the dashed cyan line. This footpoint source then migrates up the flux tube, following the evaporation front. This can be tracked as the dashed orange line in the middle panel of Figure 6. As the fronts merge and pass through each other from 25 to 30 s, the source stays near the apex, then follows the fronts down the loop. The sources do not follow the front all the way to the footpoint, instead stopping at roughly halfway between the apex and footpoint before rising back to the apex. This can be seen in the middle panel again as the dashed pink line now sits in the middle of the red downflow path from 50 to 80 s on the left side of that plot. This oscillatory motion would likely continue until the tube cooled or drained such that X-rays are no longer observable. While the motion is noticeable, the density fluctuations smooth out over time, leading to the damping of the source motion, as seen in panel (c).
In comparison with the observations in Figure 1 we can see some differences with this simulation. The observed sources drop from the apex to the footpoint, on the order of 30 s. They then rise back to the apex over the following 40 s. In contrast, the simulated source motion shows a much faster motion, with the initial drop and rise taking less than 30 s in total. The observations also show a separation between the bands as they move down the loop. The 8-10 keV band is closer to the footpoint of the loop than the lower energy bands. The simulation produces the opposite results, with the higher energy bands remaining nearer to the apex of the loop. This is expected given that the synthetic X-ray spectrum is dominated by the thermal emission below 10 keV. The electron beam heating the chromosphere is primarily increasing emission in the 3-6 keV band. The higher energy bands have stronger emission coming from the coronal segment where the plasma is hottest. This moves the centroids for those bands higher up the loop than the 3-6 keV band.
At later times, the observations show more overlap between the bands. During the initial evaporation from 7 to 29 s, we can see the 3-6 keV band is split from the two higher bands. For the following 30 s the centroids are cospatial. At 60 s and for the remainder of the simulation, there is a split of the bands with the lower energy being lower on the loop. This later evolution is consistent with observations, which show a transition from 8 to 10 keV being the lowest to the highest along the loop, with the 3-6 keV being the closest to the footpoint.

Results and Parameter Space
To determine the impact of the energy partition on the evolution of the X-ray sources, we ran a selection of simulations, shown in Table 1. The initial run (run 8) was set up based on results from O' Flannagain et al. (2013) and  in terms of its length and electron-beam characteristics. This simulation resulted in similar behavior in the spectral observations, showing the rise and fall of the nonthermal emission at the higher energies in the observation. Based on this, we built a retraction-powered simulation (run 1) that released the same amount of usable energy as the beamheating run used for the beam.
With those two simulations set, we chose value of the drag parameter that left the retraction with 10%, 25%, 50%, or 75% of the energy in run 1. For each of these drag values, we took the beam heating and reduced its amplitude to match the missing fraction of energy from drag. This creates a set of runs with the partitioning of the energy between the beam and retraction. For each of these energy splits, we varied the low-energy cutoff as shown in Table 1, creating an ensemble of 25 simulations. In addition to this collection we did consider changing the power-law index. This variation did not produce any significant change in the behavior of the sources, and as such we used the value of 5 for all our runs. This value lies in the range found by O'Flannagain et al. (2013) and used in . It is worth noting that for the low-energy cutoffs used in these simulations, the thick target model has shortcomings that could be addressed by using a warm target model (Jeffrey et al. 2019). However, the strong conduction front sweeping down the loop from the retraction site will move substantial energy to the lower atmosphere that beams alone may miss.
The runs in Table 1 all share the same total energy input. We vary the proportion of energy going into either retraction or beam heating, but the total input is the same. Another approach would be to consider the beam-input energy to be fixed and add a variable amount of retraction energy to this. While the following analysis chose the former, we did test the latter and found that it produces qualitatively similar results. As we increased the retraction energy, we found worse correlation between the observed and simulated motions. In part we made the choice to keep the total energy constant because comparison between the various runs in Table 1 are more diagnostic of the energy partition when using the same total energy. In addition, there is no accepted value for the electron flux in this event. The energy release by the retraction of the flux tube away from the reconnection site to the flare arcade is similarly hard to pin down, requiring observation of the loop shortening and a value for the coronal magnetic field. Unfortunately, this event was only observed by RHESSI, so additional constraints are unavailable. Starting with the runs in Table 1 we consider which would be the best fit to observations. We processed each run as described Section 3. Then we considered the X-ray spectrum at the peak amplitude of the electron beam. The observations showed that the emission above 10 keV was consistently nonthermal. If a run did not show nonthermal emission as the dominant source of emission above 10 keV, it was considered a bad match to the observations. The runs that match this characteristic of the observations are shown in Table 2.
The remaining runs' source motions are directly compared to the observed source motion. For this we calculate the cross correlation of the observed X-ray source movement with the synthetic X-ray source movement in each band. This gives us a quantitative measurement of how each band matched with the observations but also the time offset between the two. We average together the values for each of the bands to come up with a single score for each run. This procedure leaves us with the values shown in the second column of Table 2. Additionally, we consider how well the cross-correlation results fit with each other. To that end we introduce the spread value for each of these runs. This value is an indication of how close the time shifts from each of the three cross correlations are. We use the change in time over the maximum possible shift and sum over the three pairs. A value of 0 occurs when all three cross correlations share the same time offset. The maximum possible spread between the three will result in a value of 2.
With quantified comparisons between the synthetic and observed source movements, we can determine the best match. Based solely on the correlation value, run 9 provides the best agreement with observations. However, it does this by matching one band very well, and then the other bands match worse. This run also has a large spread value compared to the other runs in consideration. The 3-6 and 8-10 keV bands had a time offset of 28 and 24 s earlier, but the 6-8 keV band showed the best correlation to its observed counterpart with a 44 s offset later. The best fitting for this run requires that the emission in different bands be shifted more than a minute apart from each other to reach the best correlations, leading us to rank it lower than the correlation alone would suggest.
The three remaining runs are very close in both their correlation and spread values. All three runs show at least half the energy release in the simulation coming from the electron beam, and all of them have a delay of around 40 s to match the observations. Run 21 has a marginally higher correlation than run 17 at the expense of having one band mismatched by 4 s compared to run 17, which has all bands at the same time offset.
Overall, we can see several trends from the runs in Table 2. Of note is all but one of these runs show that the best fit is a shift of 36-40 s. This shift matches the observed dip and rise to the simulated evaporation fronts passing through the apex and back to the footpoints. The comparison of the observations from Figure 1 and run 21 is shown in Figure 7. The left panel overlays the two sets on top of the other, and they initially look out of phase. The right panel shows the 6-8 keV band of each with the observations in solid orange and the simulation in dashed black. The time-shifted simulation matches the observed movement. However, the simulated source does not proceed all the way down to the footpoint in later motion unlike the observed movement. It is worth noting that our initial height definition is somewhat arbitrary and that our source locations could be adjusted left or right relative to the observed source height in Figure 7. This would provide an overall better match but would not address the deficiencies with the motion in the simulation not reaching as far down the loop as the observations. This lack of depth results from two factors. This first is that the evaporation front is a less defined signal against the background of the filled flux tube, reducing the effect it has on the centroid locations. The second is that the front itself spreads out over time, reducing the enhancement from the front as well. With the widths of both observed and synthetic sources taken into consideration at the shifted time of the simulation, the matching timescales in the shifted data lend support to the time shift.
The shifted simulation implies two facts about the evolution of the X-ray sources. The first is that these sources do not necessarily carry any information on the retraction of the flux tube. In these simulations, the flux tube starts with an excess length 23% more than the final length. If the sources started near or at the apex, instead of an extended source down the leg, then the retraction would potentially be diagnostic from that motion. The second point is that the observed source motions are consistent with the density fronts in an evaporatively filled flux tube. Both of these points imply that the "initial" or early X-ray emission is coming from the tube that has been energized some 30-40 s earlier. Evaporation has filled the loop with plasma, and it is now in the radiative cooling phase.

Discussion
The foregoing analysis has explored a novel explanation for the observed movement of X-ray sources in flare loops. We considered an ensemble of 1D simulations to evaluate the effect of partitioning energy between retraction and electron-beam heating. This was done through the simulation of multiple flux tubes using the TFT model and varying the strength of the drag and beam heating. Following this, we created synthetic X-ray observations of those flux tubes through CHIANTI. We summed through height and masked the synthetic X-rays to match the approach used in O' Flannagain et al. (2013) and  before calculating a centroid to track the location of the source. To consider the best partition of energy compared to the observations, we cross-correlated the synthetic and observed source locations. Systematically, we found that the best fit of the observations matched with the behavior of flux tube already filled by evaporation and that at least half of the energy delivered to the flux tube results from beam heating, while the remainder heats the corona through retraction power jets as described above.
The expected behavior was that the fall and rise are the result of the beam properties changing, shifting the heating from near the apex of the loop down toward the footpoints. The duration of the sources dropping from the apex to the footpoints and rising again ranges from 60 to 70 s depending on the channel used. At 40 s in the simulation, the sound transit time from the chromosphere to the apex is 38 s, and as the simulation advances through to 110 s that transit time has increased to 50 s. This increase in transit time mimics the movement of the observed sources slowing down suggesting that this motion could be the result of sound waves in a filled flux tube. In these simulations, the longest time for the rise and fall attributable to the electron beam is around 30 s. We found no way to substantially slow the fall time to match observations. Given this, we are left to conclude that this motion is from already evaporated plasma in the flux tube. In addition to ascribing the motion of the X-ray sources to bulk plasma motion, we can conclude that the addition of retraction to the dynamics did improve the simulations. Even with high drag reducing the effective energy of the retraction, the dense high temperature plug of material at the apex served to keep that area emitting substantial X-rays until retraction had finished. Without this effect, the emission would fall to the footpoints almost immediately, as in . Furthermore, while in this work we only considered one change in length, using an initially longer flux tube to extend the initial period of apex-dominated emission is possible. This would not improve the comparison to observations here as the source motion would show evidence of retraction not seen here.
Including retraction in flare-loop models allows for the creation of features that would otherwise be difficult to replicate. Observations show the presence of superhot hard X-ray sources above the flare loop tops (Veronig et al. 2006;Caspi & Lin 2010;Caspi et al. 2015). There are multiple effects in two-and three-dimensional models that can produce this, such as termination shocks (Shen et al. 2018;Kong et al. 2022), the successive reconnection provided by magnetic islands fragmenting (Karlický & Bárta 2011), and fluid instabilities (Fang et al. 2016). These effects can be difficult to model in 1D, and typical electron beam models have a hard time producing the local density and temperature enhancements in the corona to reproduce these observations (see Figure 6 in Reep et al. 2019 and Figure 1 in Kerr et al. 2020). The plug of hot dense plasma created in PREFT is capable of reproducing such a hot source during retraction.
This approach does have some shortcomings. The analysis considers the flare to be composed of a single loop. This is a gross oversimplification of the process but given the geometry, a necessary one. We attempted to mitigate this by considering multiple heating events in two different ways. We found that doubling the duration of the electron beam from 20 to 40 s improved the fit in the 90% and 100% beam cases. In this sense, the faster retraction occurred, the worse it was to consider a longer beam duration. This was also the result when decoupling retraction from the electron beam. Postponing the beam start time to 10 s led to worse correlations between the observations and simulations. We also considered a single loop that had multiple short-duration heating events, and the results Figure 7. Comparison of the observational emission to the best synthetic emission. Left: the observational data taken from Figure 1 in solid lines and the synthetic emission derived from run 20 in dashed lines. In both cases, the red shows 3-6 keV, the green 6-8 keV, and the blue 8-10 keV emission. Right: 6-8 keV observed emission from the left in solid orange. The simulated data in dashed black lines has been shifted vertically 40 s according to the maximum correlation between the two. The timescale of the evaporation fits well with the observed movement, but the motion is too small. showed no improvement. Finally, we considered blending four of the best flux tubes together with an offset of 20 s. This, too, performed similar to the single-flux-tube approach, and we did not consider further multi-loop attempts. None of these attempts produced any meaningful improvement in agreement with observations despite adding significant complexity to the approach used in each. These efforts led us to conclude that whatever the mechanism responsible for causing the electron beam, at least in this case it is better if the timescale for the beam is consistent with the duration of retraction. This suggests that there is a link between the flux tube retraction and the electron beam acceleration mechanism.
In the future, we could improve on this work through a few efforts. Improving the capability of the beam heating in PREFT to allow for a time-varying power-law index and considering ionization effects may improve accuracy in where the beam energy is deposited as a function of time. In general, the use of multiple loops would also help to provide for a more comprehensive response to the electron beams, including the use of different beams for different flare loops at different times. This approach would allow us to overcome the lack of nonthermal emissions later in the observations, as the current approach results invery early and short-duration nonthermal emission. However, this approach depends on favorable geometry and line of sight in the observations. In addition to this, allowing for turbulent suppression of thermal conduction (Bian et al. 2016;Emslie & Bian 2018;Allred et al. 2022) could potentially keep the hot, dense source at the apex for longer than the current model. We could also consider Alfvénic-wave heating mechanisms (Reep et al. 2018;Ashfield & Longcope 2023), given that retraction in PREFT has Alfvénic signatures.
In light of these results and shortcomings, we believe that more work on the X-ray appearance of flare loops is required. If the movement of these sources is due to beam properties, then understanding that relationship will give insight into the mechanism responsible for the electron acceleration. Care should be taken, however, as it appears the movement of the X-ray peaks is not solely due to the properties of the electron beam, and ignoring these factors is likely to skew results. A promising avenue then would be to continue using multiple 1D models probing the expected responses of the flare loop.