Tracking the Source of Solar Type II Bursts through Comparisons of Simulations and Radio Data

It has long been known that the solar corona is capable of producing a variety of intense radio bursts in narrow frequency bands with different characteristic variations in frequency, time, and intensity, depending on the source (Wild et al. 1963). Among these bursts are type II and type III radio bursts, which were first observed in the 1950s with ground-based radio observatories in frequencies slightly above the 15 MHz ionospheric cutoff (Wild 1950). These observations continued into the space age, with single spacecraft observations like Wind (Bougeret et al. 1995) observing these bursts in the lowest frequencies below 15 MHz. The Wind spectrum in Figure 1 illustrates typical examples of both type II and III bursts, identified by their narrowband structures that drift to lower frequencies, with type II bursts drifting more slowly than type IIIs. Type II and III solar radio bursts have been observed as being associated with the main types of solar energetic particle (SEP) events: gradual and impulsive, respectively (see the review by Reames 2013). Type III bursts are more numerous than type II bursts (Cane et al. 2002), but last for a shorter amount of time and have a lower SEP fluence than the gradual event-associated type II bursts (Reames 1999). SEPs consist of high-energy ions and electrons, generally in the keV–GeV range, and are thought to be produced by solar flares and coronal mass ejections (CMEs) (Tsurutani et al. 2009). The most intense SEP events are caused by CMEs as they accelerate away from the Sun and plow into the extended solar corona (Reames 2013). These SEPs are a form of space weather that poses a threat to both satellites and humans in space. Energetic particles can access the Earth system through the polar regions, potentially causing harm to satellites and passengers in jets on polar routes. The higher the flux of the incoming energetic particles, the further down in latitude their effects can spread. In near-Earth space, energetic particles can cause satellite damage and harm astronauts that are not protected by Earth's atmosphere or other means (Council 2008). Interplanetary CMEs (ICMEs) and their shock and sheath regions are also the leading drivers of geomagnetic storms (Ontiveros & Gonzalez-Esparza 2010; Kilpua et al. 2019). Understanding the physics behind these events is crucial to predicting their impacts ahead of time, to protect human lives and property, and untangling the relationship between SEPs and type II and III radio bursts is instrumental in this endeavor.

Zoom In Zoom Out Reset image size

Electron beams associated with type III bursts were tracked in the 1970s by groups such as Frank & Gurnett (1972) and Alvarez et al. (1972), who found Langmuir waves associated with the electron beams were causing the type III emission (Gurnett & Anderson 1976, 1977). In addition, Lin et al. (1981, 1986) showed that these plasma oscillations occurred alongside a bump-on-tail electron energy distribution, indicating that a plasma emission mechanism was behind the type III radio emission. This same mechanism is thought to be behind the emission of type II bursts as well (Thejappa et al. 2007).

The prevailing physical theory of the plasma emission mechanism was first put forward in 1958 by Ginzburg & Zhelezniakov (1958), and it has since been refined by them and many others over the years, as reviewed by Melrose (2008) and the references therein. The basic idea is that the radio bursts are created in a multistep process that starts with an electron beam or a subpopulation of electrons that are much faster than their neighbors. In the case of type III bursts, relatively small packets of jettisoned particles from flaring active regions on the surface of the Sun make up the beam, with typical electron energies of 10–100 keV (Reames 2013). For type II bursts, the beam is thought to be formed by particles somehow energized by a CME/interplanetary shock (Cane & Stone 1984; Cane et al. 1987), with typical electron energies on the order of ∼1 MeV (Cliver & Ling 2007). Once the electron beam forms, it is unstable to wave–particle interactions that can generate Langmuir waves with frequencies around the local plasma frequency, via the standard bump-on-tail instability. The waves that are the right frequency to resonate with the beamed electrons convert energy from the electron beam to the plasma waves, relaxing the distribution function to a smoother profile. Through a combination of three-wave decay, ion scattering, and coalescence, a fraction of the power in the Langmuir waves can be converted into radio waves at the corresponding local plasma frequency, or its harmonic. These radio waves may then escape from their source region and propagate into space.

The local plasma frequency is only dependent on the density of the electrons n_e , being proportional to $\sqrt{{n}_{e}}$, and so the frequency of the burst trails down through frequency space as the source energetic electrons soar outwards from the Sun. An example of a Wind/WAVES radio spectrum (Bougeret et al. 1995) containing both types of emission is shown in Figure 1. Both types of bursts are identified by their narrowband structures that drift to lower frequencies, and they may be differentiated by their frequency drift rate, with type II bursts drifting more slowly than type IIIs. The full formula for the local plasma frequency ω_pe is shown in Equation (1), where the elementary charge of an electron e (in the electrostatic unit of charge, esu) and the mass of an electron m_e (grams) are constants, and n_e is the local electron number density per cubic centimeter.

$$\begin{eqnarray}&&{\omega }_{{pe}}=\sqrt{\displaystyle \frac{4\pi {n}_{e}{e}^{2}}{{m}_{e}}}=8.98\,\mathrm{kHz}\sqrt{{n}_{e}/{\mathrm{cm}}^{3}}.\end{eqnarray} \tag{ 1 }$$

Type II emissions can be classified by their emitted wavelengths, ranging from metric (m), to decametric–hectometric (DH), all the way to kilometric (km) wavelengths. Type II emissions can start at frequencies around 150 MHz, which corresponds to the plasma frequency a few solar radii from the Sun, and can drift past 1 au, where the local plasma frequency is approximately 25 kHz. Gopalswamy et al. (2005) found evidence that the kinetic energy of the CME seems to determine the lifetime of type II radio bursts, with the most energetic events having bursts that extend from the metric frequencies down to the DH and kilometric ranges. Less energetic CMEs on the other hand, have no type II emission, or an emission that stops in the metric frequencies, before the DH wavelength ranges. The wide-ranging surveys of events by Gopalswamy et al. (2008) found that basically all the DH type II bursts are associated with SEP events at Earth if one accounts for magnetic connectivity. The intensity of the associated SEP fluxes also increases with CME speed and width (i.e., kinetic energy). Another recent broad-ranging study by Winter & Ledbetter (2015) further solidifies this relationship, as it is reported that every single SEP event seen by Wind/WAVES from 2010 to 2013 with a peak >10 MeV flux above 15 protons cm⁻² s⁻¹ sr⁻¹ is associated with a DH type II burst, most of which occur alongside a type III burst.

With this general relationship having been cemented, the primary goal of this work is to better understand which plasma parameters are most important for type II burst generation within the context of CME structure and shock geometry. To do this, a methodology is developed to take CME simulation data and create synthetic type II bursts for hypothesized burst locations across the entire simulation time. Existing methodologies, such as those from Reiner et al. (2007), for using type II radio data to profile the velocity and acceleration of a CME are extended in order to create a "stencil," which is an estimated spectral time profile employing a Gaussian intensity profile fitted over the observed burst intensity. It is this Gaussian profile fit over the observed radio burst that the synthetic spectra are scored against. These two methodologies are then used as the cornerstones of a framework in which the following optimization problem is posed: "What progression of simulated burst locations provides a synthetic spectrum that best matches a single spacecraft-observed spectrum?"

Within this framework, one may track different features of the simulated CME to test different hypotheses of what plasma parameters are causal in creating the type II bursts around CMEs. These burst locations in turn mark where the acceleration of energetic particles is happening as the CME expands throughout the heliosphere. These burst locations can then be fed into energetic particle codes like Energetic Particle Radiation Environment Module (EPREM; Schwadron et al. 2010, 2014) that track the magnetic connectivity with Earth to forecast the energetic particle flux at Earth, although that is not done in this work. One promising long-term goal of applying this framework is finding a common set of thresholds of plasma parameters that creates synthetic type II spectra that match the observed spectra across many recreated simulations of major CMEs. Such findings would be valuable information for constraining physical theories of particle acceleration at CMEs, and would also be immediately applicable to real time space weather operations like those at NOAA's Space Weather Prediction Center and others. A set of key thresholds that marks the site of type II burst generation and particle acceleration would help enable more reliable space weather prediction, especially with the use of faster than real time simulations of incoming CMEs on the horizon.

In this work, these methodologies and the optimization framework are developed with a focus on a single case study in the form of the well-documented 2005 May 13 CME. A general outline of the remaining sections of the paper is as follows. In Section 2, previous work in tracking CMEs and type II bursts is reviewed, with major steps in methodology and subsequent understanding being identified. Section 3 describes the simulation of the CME by Manchester et al. (2014), performed with the Alfvén Wave Solar Model (AWSoM) developed by van der Holst et al. (2014), that is used throughout the rest of the paper. Section 4 develops the methodology to create an event-specific "stencil" of a type II burst from the observed single spacecraft data. In this case, the methodology assumes CME velocity and solar wind density profiles to fit the evolution of the type II burst over time. This approach yields, for each moment in time, an idealized Gaussian spectral profile around a central frequency, which is used to create the "stencil" that the synthetic spectra are scored on according to their similarity with the observed Wind/WAVES radio data. Section 5 develops the methodology to create synthetic spectra from the simulation, using the upstream electron density to determine the observed frequency profile for each moment in time. In Section 6, the methodologies are applied to frame the optimization problem: "What progression of simulated burst locations provides a synthetic spectrum that best matches the single spacecraft-observed spectrum?" Here, a range of CME model features, based on different physical thresholds and geometrical factors, are tested to identify the likely site of burst generation around the CME. In Section 7, possible sources of error and objections to the new methodologies are discussed and addressed, to establish that the methodologies are appropriate for this optimization framework. In Section 8, the results from this investigation are discussed and future work is laid out.

2.1. Theories of Particle Acceleration

2.2. Advances in Space-based Observations

2.3. Advances in Ground-based Radio Observations

2.4. Advances in MHD Simulations

2.5. Using in situ Plasma Data

The CME simulation (described in Manchester et al. 2014) was performed using the AWSoM (van der Holst et al. 2014), which models the time average ambient solar wind over a Carrington Rotation, based on photospheric radial magnetic field measurements from synoptic magnetograms. This level of detail is necessary for accurate simulations of upstream solar wind density, as Manchester et al. (2005) showed that the ambient solar wind structure strongly affects the evolution of the CME structure. The simulated CME is driven by a Gibson & Low (1998) magnetic flux rope that has some common preevent characteristics, such as a dense core surrounding the flux rope and a dense helmet streamer with a cavity. This particular CME has been well studied, with works like that of Bisi et al. (2010) combining data from several suites of instruments to piece together a more complete picture of the event. A SOHO C2 coronagraph showing the "halo" CME for this event alongside its corresponding Wind/WAVES radio spectrum is shown in Figure 2.

Zoom In Zoom Out Reset image size

The AWSoM model, recently developed at the University of Michigan (van der Holst et al. 2014), includes a phenomenological model of low-frequency Alfvén wave turbulence to heat the corona and accelerate the solar wind. At the inner boundary, the Alfvén wave energy is injected with the Poynting flux proportional to the magnetic field strength. Counterpropagating waves are present on open field lines, due to Alfvén speed gradients and field-aligned vorticity reflecting some of the outward-propagating waves backwards. Wave energy enters both ends of closed field lines to produce near-balanced turbulence at the apex of the loops. AWSoM currently has three temperatures—an isotropic electron temperature and separate parallel and perpendicular ion temperatures–although the results here were obtained from a two-temperature simulation, one for electrons and a second for isotropic ions. To distribute the heating of the turbulence dissipation across the different temperatures, the model uses the linear wave theory and nonlinear stochastic heating as presented in Chandran et al. (2011). AWSoM also incorporates both collisional and collisionless heat conduction into its treatment of electron thermodynamics. This all leads to a more realistic simulation that better matches the real data (Sachdeva et al. 2019).

The results from the simulation of the 2005 May 13 CME event are seen in Figure 3. Here, different physical quantities relevant to SEP/acceleration and radio emission are shown 20 minutes after initiation. The panel on the right shows the magnetic flux rope being expelled from the corona, with a peak speed of nearly 2000 km s⁻¹, which closely matches the speed fit by Reiner et al. (2007) for this same event. The left panel shows the full extent of the region of the simulation, with local densities that are at least 3.5× the preeruption values to identify the density-derived shock, colored to show the upstream plasma frequency. The central panel shows a region of strongly shocked plasma, identified by selecting regions of the simulation with entropy values that are at least 4× the preeruption values, with the surface colored to show the shock angle Θ_Bn . Entropy is defined as $S=\tfrac{p}{{n}^{5/3}}$ for plasma number density n and pressure p. The right panel shows only the magnetic field configuration. In all panels, the magnetic field lines are colored to show the plasma frequency.

Zoom In Zoom Out Reset image size

The data from the MHD simulations contain a wealth of information that would be impossible to get at similar cadence and resolution with in situ or remotely sensed data. Each data point contains cell-averaged quantities of magnetic field, plasma density, momentum, and internal energy. These parameters may then be used to calculate second-order parameters like entropy, de Hoffmann–Teller frame velocity, the angle between the shock normal and upstream magnetic field Θ_Bn , and more. Some of these derived plasma quantities have been found to be highly correlated with Langmuir waves upstream of shocks in data from in situ observations (Pulupa et al. 2010). Model data extractions (isosurfaces and field lines), using these plasma parameters, along with geometrical subselections will be used in the following sections to test different hypotheses of burst locations by comparing their synthetic spectra with the observed Wind/WAVES spectrum.

Synthetic spectra do not have physical units such as V² or spectral brightness, like spacecraft data, as they are not made by directly simulating the physics of burst generation. Briefly put, a synthetic spectrum is created by taking a suspected source region from an MHD simulation, calculating the upstream plasma frequency for each point, and plotting the resulting distribution in a histogram for each time step, with the frequency bins being the same as those of the Wind/WAVES data. As such, it is not feasible to directly compare synthetic spectra with spacecraft-observed spectra, as they have dimensionless units, consisting only of raw counts of simulation points with a certain plasma frequency, compared to the power received at an antenna in a given frequency. Thus, it is necessary to create a custom "stencil," i.e., a distribution over frequency and time derived from the observed data, so that one may analytically assign a similarity score between the synthetic spectrum and the observed spectrum, using the stencil as a proxy. This methodology gives one the flexibility to be agnostic regarding the actual physics generating the type II burst, assuming only that the emission frequency is the upstream plasma frequency. This forms the basis of a framework in which one can find and test hypothetical burst locations by quantifying how well the formed synthetic spectra match the stencil.

Note that this framework would not be possible without the assumption of the emission frequency being an upstream plasma frequency. It has been suggested by Bastian (2007) that some events have a type II-like emission that comes not from the plasma emission mechanism, but from incoherent synchrotron radiation from near-relativistic electrons around the CME or its sheath. If this were the case, one could potentially back out the distribution of energetic elections, but could probably not constrain their location, since the emitting frequency would not be tied to a location like it is for the plasma emission mechanism, where the upstream plasma frequency is affected by the distance from the Sun.

Since every CME is unique, having a different geometry and speed profile, each event requires a custom stencil to be created that the synthetic spectra can be compared to. In the current literature, there are examples of using basic parameterizations of type III bursts to fit the underlying kinematics of the electron packet driving the burst. Dulk et al. (1987) simply fit a constant speed of the outgoing energetic electrons, while Reiner & MacDowall (2015) included a constant deceleration term alongside an electron density and Parker spiral model to accurately capture the radio wave travel times. They argue that the additional degree of freedom allows their model to better fit the frequency drift from the spacecraft-observed spectrum as well as match the arrival time of the Langmuir waves or local emissions. Reiner et al. (2007) characterized the kinematics of CMEs from the type II frequency drift, using SOHO LASCO white light observations to constrain the initial speed of the ICME. One of the events they analyzed was actually the same 2005 May 13 event studied in this work. Here, this methodology is extended to create a physically motivated stencil based on CME kinematics.

The initial kinematics of the 2005 May 13 CME were derived from the projected distances and timestamps of observations of the event with SOHO LASCO C2 (17:22 UT) and C3 (17:42). These times correspond to an initial speed of 1690 km s⁻¹ and an extrapolated liftoff at 16:47. This initial speed, however, is not the true speed; it is only a plane of sky-projected speed. Since the event in question is a halo CME, the true speed is higher than the plane of sky estimate. The event was detected in situ by the Wind spacecraft at 1 au on 02:13:30 UT on 2005 May 15, implying an average speed of 1246 km s⁻¹. In addition, Reiner et al. (2007) found the measured radial shock speed at 1 au to be 778 km s⁻¹. This large deceleration between liftoff and 1 au is expected for fast CMEs; Woo et al. (1985), Sheeley et al. (1999), and Manchester et al. (2004) have shown that this deceleration increases with the initial speed of the CME, explaining the rather narrow range of speeds of shocks observed at 1 au.

As Reiner et al. (2007) describes, there are many possible simple kinematic solutions satisfying these constraints, with higher initial speeds being matched by quicker deceleration for shorter periods of time, before settling at 778 km s⁻¹ at 1 au. They also note that the structure of the type II emission may be used to find the "true" velocity profile. By assuming a solar wind density model, one can convert a kinematic velocity profile for a CME-driven shock to a predicted radio profile over frequency and time. Reiner et al. (2007) opted for a simple r⁻² density model, determined by a prescribed density at 1 au. They then adjusted the parameters of the kinematic profile until they found a good fit by eye to the center of the apparent type II event as the fundamental or harmonic emission. This kinematic model had v₀ = 2500 km s⁻¹, a deceleration of 26.7 m s⁻¹ s⁻¹, and a deceleration time of 17.9 hr. The frequency profile described by this solution will be central to our own scoring stencil below. With this kinematic solution, one can see in Figure 4 that in the first 2 hr the type II emission is primarily at the fundamental plasma frequency, while after 19:30 or so the harmonic emission is dominant, continuing for another 31 hr before the shock arrives at 1 au, coinciding with the in situ observations of the shock by multiple spacecraft like Wind and ACE (Bisi et al. 2010). This arrival of the shock may also be seen in the Wind/WAVES data as an enhancement to the type II burst at around 50 hr on the bottom panel of Figure 4.

Zoom In Zoom Out Reset image size

One may also extend this characterization of the emission, fitting by eye a customized envelope created by manually choosing three time–frequency profiles using this quadratic model to partition the initial fundamental emission from 17:00–19:00. This can be seen in Figure 5(a), where the three lines surround the strongest part of the primary type II emission. The center frequency line from the Reiner et al. (2007) fit is used as the target frequency for each time step, with the lines surrounding it used to impose a Gaussian profile over the observed burst. This Gaussian profile is normalized so that the middle frequency has weight 1. For every moment in time, each half of the identified emission in frequency space is given a different sigma value to determine its frequency width, and is defined with 3σ being the upper/lower end of the by-eye fit of the data. The simulation data used for this work is constrained to the first 2 hr of the event, which coincides with the period of the event in which the fundamental emission is dominating. Further studies would require a straightforward extension to the stencil in order to include the harmonic frequencies.

Zoom In Zoom Out Reset image size

The resulting Gaussian profile over frequency and time is shown in Figure 5(b), and it can be thought of as a stencil that a synthetic spectrum from an MHD simulation may be compared to and scored from, enabling a template-matching technique that may be used to solve the inverse problem of finding the site of the type II burst generation. The stencil is dimensionless, having no physical units, and thus it may be directly compared to synthetic spectra if they have the same time and frequency structure. An example of the Wind/WAVES spectrum overplotted with the described Gaussian stencil is seen in Figure 6, 10 minutes into the type II event at 17:10. One can see in this figure that the type II burst is far from the brightest feature in the noisy data, making manual partitioning necessary for the creation of the stencil. A global element-wise multiplication between the synthetic spectrum and stencil, followed by a global sum of all the channels over frequency and time, will yield a similarity score. This similarity score is linearly proportional to the fractional amount of data points in a simulated CME feature whose upstream plasma frequency in its synthetic spectrum overlaps with the scoring stencil. The generation of these synthetic spectra will be further explained in the following section.

Zoom In Zoom Out Reset image size

The primary goal of this work is to better understand which plasma parameters are most important for type II burst generation within the context of CME structure and shock geometry. Making this connection will lead to greater understanding of how gradual SEP events are created by CMEs. By employing state-of-the-art AWSoM simulations, one can recreate historical major space weather events. Sophisticated simulations can provide synthetic in situ data over the entire simulated domain. Notably, this means that one can derive the plasma parameters when the CME is within 20 solar radii of the Sun, when most of the SEP acceleration happens. By taking different data extractions of the simulation, based on physical thresholds for different physical processes or based on CME geometry, synthetic spectra can be made and compared with the actual radio data. This radio data needs to come from space, since the lowest parts of the type II and III bursts emit below the Earth's ionospheric cutoff. By matching the synthetic spectra with the spacecraft-observed data, one may attempt to determine which physical thresholds must be exceeded in order for burst generation to occur. The confirmation of similar physical thresholds applying to synthetic spectra that generally match observations across multiple events would pin down the specific physics triggering the burst generation, as well as the location of the particle acceleration. This work is focused on a single radio-loud event from 2005 May 13 at 16:47 as a test case for this new framework.

To this end, the first 2 hr of an MHD simulated recreation of the radio-loud CME is compared with the Wind/WAVES radio spectrum to identify the likely regions in which the type II bursts are generated. The data extractions from the simulation are based on combinations of physical thresholds, as well as the CME/shock geometry, including solar longitude and latitude. From these data, synthetic radio spectra are made and compared to the stencil described in the previous section.

Subsets are extracted from the 3D data files of the AWSoM simulation, typically isosurfaces that meet some physical criteria. These data files contain some basic plasma parameters that can be used to calculate the secondary plasma parameters at the included simulation points, including the upstream plasma density used to calculate the upstream plasma frequency. The four main data types used here are: isosurfaces where points with at least a 3.5 MK temperature are used to identify the CME current sheet; shocked plasma with a density compression factor of 3.5 or more; and points where the entropy changes by a factor of 4 or more, indicating the strongest (highest Mach number) shock; and data where the Fast Magnetosonic Mach number is 3 or more. Each of these data features is produced with a time cadence of 2 minutes for the first 40 minutes of the simulation, later transitioning to every 5 minutes, producing a total of 35 data files over 120 minutes. 2 hr of simulation time brings the CME out to 20 solar radii, corresponding to a local upstream plasma frequency around 300 kHz. The one exception to this simulation sampling cadence is the fast magnetosonic Mach # ≥ 3 subselection, which, due to a lack of resources, only combines the six times of [10, 20, 30, 40, 60, 90] minutes after the start of the simulation. This results is a coarser data set that may still be used to inquire into the usefulness of this parameter, defined as the ratio of the shock speed to the fast magnetosonic speed, which has been shown to be associated with energetic particle acceleration at values >3 (Rouillard et al. 2016). With these broad data sets in hand, finer extractions can then be applied, based on plasma parameters considered to be important for SEP acceleration, such as the upstream velocity magnitude and the de Hoffmann–Teller frame velocity (Pulupa et al. 2010), as well as subsets based on the geometry and spatial location of the data. Synthetic spectra are created by calculating the upstream plasma frequency for each point, and then plotting the resulting distribution in a histogram. This is done for each time step, using the same frequency bins as the Wind/WAVES data. For every time step of data saved, some preliminary processing must be done. Depending on the simulation domain, different coordinate transformations must be performed.

The SC component of the MHD simulation goes out to 20 R_⊙ and is done in a rotating reference frame; the Space Weather Modeling Framework (SWMF) calls it Heliographic Rotation, while SunPy calls it Heliographic Carrington. The Inner Heliosphere (IH) simulation domain starts at 20 R_⊙ and can go out to several au. The IH domain is simulated using an inertial reference frame; SWMF calls it Heliographic Inertial, while SunPy calls it Heliocentric Inertial. Data from both of these domains are transformed into the Heliocentric Earth Equatorial (HEEQ) frame, where the X-axis is the Sun–Earth line, and the Z-axis is the solar north pole. These coordinate systems are described fully in Thompson (2006).

The HEEQ data for each time step is then sorted by x coordinate, and a routine is applied to cull any low-frequency points that appear to be behind higher frequency points from the Earth's perspective. This is done efficiently by employing fast algorithms to conduct global nearest neighbor calculations. Starting with the largest x coordinate (closest to Earth), the plasma frequency is checked, then compared to any sunward neighbors within a given angle to the sky, chosen here to be 0.001 solar radii from the vantage of Wind at L1, corresponding to 1'' in the plane of sky from Earth. Any neighboring data points within this obscuration angle with a lower plasma frequency value than the earthward neighbor are culled. This reflects the physical fact of the absorption of these lower frequency waves when passing through higher density areas. Wind/WAVES was parked at L1 on the Sun–Earth line during this time period, so this line-of-sight culling is best done in this HEEQ frame, to match the observed spectral data. Secondary scattering effects, which can enable lower frequency emissions to take non-line-of-sight paths around higher density areas that would normally absorb them, are not considered in this work. The remaining points can then be further reduced using customized physical cutoffs to investigate the resulting distribution in upstream plasma frequency. A special correction is also applied to the entropy subselection to remove artifacts in the data that are much closer to the Sun than the rest of the extracted data. These points have upstream magnetic field magnitudes on the order of 1 Gauss, as opposed to the rest of the data, with values of less than 0.1 Gauss, and constitute around 1% of the total data in the subselection. This correction is applied to all figures in the paper using the entropy subselection.

To form the base of a synthetic spectrum, a histogram is created for each time step with the same frequency bins as the Wind data. This histogram constitutes a single column of the synthetic spectrum. The entire synthetic spectrum is a 2D histogram, with each column sum normalized to 1. This synthetic spectrum is then interpolated to the same time cadence as the underlying data (1 minute for Wind/WAVES), while still enforcing the column sum normalization. An example of this interpolation process is shown in Figure 7. Then, by employing an element-wise multiplication between the stencil and the interpolated synthetic spectrum, and globally summing the resulting array, one receives a score for how well the synthetic spectrum from the selected simulation feature matches the spacecraft-observed data. Since the column sum of the synthetic spectrum is 1.0, and the stencil maximum is normalized to 1.0, a maximum score of 1.0/time step is achievable. The current system most rewards synthetic spectra that match the center frequency, with the Gaussian profile only marginally rewarding close-by frequencies, but a constant 1.0 weight can also be applied throughout the entire stencil envelope for a less picky algorithm. The synthetic spectra for the strong shock region are shown in Figure 7. The left panel shows the synthetic spectrum as an uninterpolated 2D histogram, while the right panel shows the synthetic spectrum interpolated to the Wind/WAVE cadences, plotted over the stencil envelope. The similarity score is seen in the title of this right panel; here, 29.26 out of a possible of 107. The maximum possible score is 107 because there are only 2 hr, or 120 minutes, of simulation data, but the simulation starts at 16:47, and the type II stencil starts at 17:00, leaving a maximum overlap of 107 minutes. The slight distortion in the synthetic spectrum around 1 MHz is due to the transition from the Wind RAD1 receiver below 1 MHz, which has a frequency resolution of 4 kHz, to the RAD2 receiver above 1 MHz, which has a frequency resolution of 50 kHz. When plotting these sets of frequency channels with different resolutions, an apparent distortion at their interface near 1 MHz can be seen, as a greater portion of points per frequency channel are seen above 1 MHz with the increase in bin width frequency.

Zoom In Zoom Out Reset image size

By analyzing in detail the plasma parameters from an MHD simulation of a radio-loud event, one may attempt to solve the inverse problem of determining the location of the burst generation site. One may extract different subsets from the simulated CME to create time-varying synthetic radio spectra using the upstream plasma density, and compare these to the actual Wind/WAVES radio data in order to determine which hypothesized source region best matches the observed spacecraft data. This bypasses the need for simulating the exact physics that are creating the bursts, and tells us what areas of the CME are most likely to be the site of particle acceleration. These locations may then be inspected more closely with the additional data from the simulation to see which plasma parameters are important for type II burst generation to occur.

It is known from the statistical analysis of type II bursts in Aguilar-Rodriguez et al. (2005) that these bursts have a small bandwidth-to-frequency ratio (BFR) of between 0.1 and 0.4 in the frequency range 0.03–14 MHz of the RAD1 and RAD2 antennas. However, the MHD recreation of this 2005 May 15 event shows that if the entire density-enhanced shock were emitting a type II burst, the frequency spread would be far larger. This effect can be seen in the left panel of Figure 3 and the top right panel of Figure 8, where the density-enhanced region of the CME has upstream plasma frequencies spanning over an order of magnitude. In this section, thresholds of plasma parameters will be used to identify and track smaller parts of the CME through space and time in order to create synthetic spectra that can better match the spacecraft-observed radio spectrum.

Zoom In Zoom Out Reset image size

6.1. CME Data Extractions

6.2. Geometrical Data Selections

One may further refine these data extractions to test the similarity of the synthetic spectra derived from different geometrically defined areas in the CME. For each moment in time, one can consider only keeping data points above some set fraction of the furthest out point, with some radial distance maxr, and discarding the rest. This would, in effect, only leave the portion of the shock furthest out from the Sun, which would typically have lower upstream plasma densities, and therefore lower plasma frequencies. One can see from Figure 9 that increasing the inner radial cutoff to 0.7*maxr (upper right), 0.8*maxr (lower left), or 0.9*maxr (lower right) in this manner can lead to a synthetic spectrum with a higher similarity score, with more of its remaining data within the envelope of the stencil. The fact that this is true for the entropy-based subselection, and not the density or magnetosonic Mach # subselections, is significant. It signals that the outer edge of the entropy-derived shock is a potential site of radio burst emission. Here, we take a moment to emphasize that this framework is a diagnostic tool, and is not expected to provide a guarantee of particle acceleration around simulated regions that provide high-scoring synthetic spectra.

Zoom In Zoom Out Reset image size

Each of the regions across the shock surface can then be subdivided even further to better isolate the potential locations of radio emission. By looking at subselections based on the solar longitude and latitude used to create the synthetic spectra, one can test hypotheses of stationary (with respect to the outward-moving shock structure) burst production sites. The synthetic spectra from two geometrical simulation features performed using the entropy-enhanced data with no global radial cutoff can be seen in Figure 10. The left panel shows the synthetic spectrum from the portion of the shock within −15° to −12° longitude and −9° to −6° degrees latitude, which clearly lies inside the defined stencil envelope. The right panel shows the synthetic spectrum from the portion of the shock within −33° to −30° degrees longitude and −9° to −6° degrees latitude, which clearly lies mostly outside the defined stencil envelope. This comparison shows how certain areas of the shock better fit the stencil of the type II burst determined from the spacecraft-observed data.

Zoom In Zoom Out Reset image size

In addition to studying synthetic spectra from single sectors of the entropy-based shock, looking at the geometrical distribution of the scores of the event can be enlightening as well. Figure 11 shows a 2D spatial color contour graph, where each cell corresponds to the similarity score for a synthetic spectrum created for a 3 × 3 degree geometrical simulation feature, normalized to the maximum score of 107 over 120 minutes of simulation data. The upper left panel of Figure 11 shows the similarity score distribution for all the nonartifact data in the entropy feature, showing that the area of the CME to the southwest of the Sun–Earth line around 6° longitude and −12° latitude has a robust maximum. This panel has a white circle and a red circle to indicate the good and bad regions, respectively, of the shock used to make the synthetic spectra in Figure 10.

Zoom In Zoom Out Reset image size

The other panels are also 2D similarity score graphs, but with progressively stricter global radial subselections, where only the data points furthest out from the Sun are counted, with a time-dependent radial threshold of 70%, 80%, or 90% of the maximum distance traveled by any point in the entropy-enhanced feature at that time. All of these similarity score distributions with various radial cutoffs share a global maximum around (6°, −12°). One can see that the similarity score distributions across the shock for a global radial cutoff of r > 0.7maxr maintain the maximum normalized score achieved with the no global radial cutoff subset of roughly 0.73, as seen in the top two panels of Figure 11. However, the r > 0.7maxr radial cutoff also has a region in the northwestern portion of the shock where the similarity scores have improved relative to the no global cutoff subset. This implies that the r > 0.7maxr cutoff gets rid of some secondary population points slightly behind the main shock, with slightly higher plasma frequencies, reducing the similarity score for those regions in the uncut subset. This r > 0.7maxr cutoff subset will be used for the remainder of the plots, unless stated otherwise. Also apparent in these panels is the feature of the shock around ( −30°, 0°) that has a constant normalized similarity score of around 0.3 across the different cutoffs, indicating that this area of the shock is at the leading edge. These eastern and southwestern edges of the shock will be more closely investigated in the following sections.

6.3. Derived Plasma Parameters

One can also use the lists of the potential key parameters from Pulupa et al. (2010) that are known to be correlated with upstream Langmuir waves to guide the selection of simulation features. This work focuses on the upstream plasma data, starting with an examination of the upstream velocity magnitude and plasma frequency across the surface of the entropy-enhanced MHD feature. Figure 12 shows the 2D spatial color contour graphs of the upstream plasma frequency and upstream velocity magnitude across the entire surface of the shock at 20, 40, and 90 minutes after the start of the simulation, displaying the mean values of the parameters for each longitude and latitude bin. The plasma frequency was calculated with the upstream simulated electron density, and is used in the creation of the synthetic radio spectra. The upstream velocity magnitude is shown as an example of a variable that is fairly well correlated with in situ observations of Langmuir waves, ranking ninth out of the 25 plasma parameters tested in Pulupa et al. (2010). The general pattern of upstream velocity magnitude across the shock front having higher values on the southwestern side of the shock, and smoothly transitioning to lower values on the eastern edge of the shock, is shared with the upstream magnetic field magnitude and upstream Alfvén speed, other plasma parameters that are fairly well correlated with in situ observations of Langmuir waves. The 2D spatial color contour graphs for the upstream magnetic field magnitude and upstream Alfvén speed are not pictured here to save space.

Zoom In Zoom Out Reset image size

Pulupa et al. (2010) found that shocks with lower upstream plasma densities are more likely to exhibit Langmuir wave activity, indicating that the upstream solar wind structure is important in forming the shape of the shock and in any acceleration of nearby particles. This is another reason that it has taken this long for simulations to get to a point where studies like this are possible; only recently have simulations, like the one conducted in Manchester et al. (2014), accurately simulated the steady-state solar wind structure of the inner heliosphere specific to the period preceding a simulated CME, before launching said CME and propagating it though the heliosphere.

Another plasma parameter classically thought to play an important role in particle acceleration at CMEs (and then Langmuir wave creation, followed by type II burst generation) is Θ_Bn , the angle between the shock normal $\hat{n}$ and the upstream magnetic field vector B. Pulupa et al. (2010) have shown Θ_Bn to be fairly well correlated with energetic electrons, with higher angles near 90° having a higher level of association with the electron beams near shocks. Another correlated parameter is the de Hoffmann–Teller frame velocity, V _HT, which was found to be the plasma parameter most correlated with energetic electrons by Pulupa et al. (2010). V _HT is defined by the following equation:

$$\begin{eqnarray}&&{{\boldsymbol{V}}}_{\mathrm{HT}}=\displaystyle \frac{\hat{{\boldsymbol{n}}}\times \left({{\boldsymbol{V}}}_{u}\times {{\boldsymbol{B}}}_{u}\right)}{{{\boldsymbol{B}}}_{u}\cdot \hat{{\boldsymbol{n}}}}.\end{eqnarray} \tag{ 2 }$$

V _HT is closely related to Θ_Bn , with the denominator, ${{\boldsymbol{B}}}_{u}\cdot \hat{{\boldsymbol{n}}}=\cos ({{\rm{\Theta }}}_{{Bn}})| {\boldsymbol{B}}|$, approaching zero as Θ_Bn nears 90°, causing V _HT to become very large. Here, V _u is the upstream velocity of the plasma in the shock rest frame. An approximation for V _u is made by taking ${{\boldsymbol{V}}}_{u}={{\boldsymbol{V}}}_{\mathrm{sim}}-1500\hat{{\boldsymbol{n}}}$, assuming for simplicity a constant shock speed of 1500 km s⁻¹ in the direction of the shock normal. Here, ${{\boldsymbol{V}}}_{\mathrm{sim}}$ is the simulated upstream plasma velocity in the regular HEEQ frame, whose magnitude is shown in Figure 12. Figure 13 shows the 2D spatial color contour graphs of Θ_Bn and V _HT across the surface of the shock at 20, 40, and 90 minutes after the start of the simulation. There is more structure within these variable contour graphs than the graph of ${{\boldsymbol{V}}}_{\mathrm{sim}}$ in Figure 12, with the higher values around the edges of the shock matching well with the structure of the similarity score contour graphs. The area of maximum similarity around (6°, −12°) coincides with an area of consistently high V _HT over the 2 hr of simulation time. The other feature of note is the area of highly elevated V _HT around the eastern edge of the shock at (−30°, 0°), with V _HT falling off after about an hour into the simulation.

Zoom In Zoom Out Reset image size

One may use these correlated plasma parameters to provide additional restrictions on the MHD features and examine the resulting synthetic spectra in order to ascertain which thresholds produce the best match with the actual Wind/WAVES spectrum. This is done in the following subsection with V _HT, as its structure across the entropy-derived shock seems to have the best match with the features seen in the uncut similarity score distribution in Figure 11.

6.4. Subselections Based on Both Geometry and Plasma Parameters

Starting at the broadest scales and working down, Figure 14 shows the synthetic spectra from a further reduction of the global radial cutoff subsets, after imposing a minimum threshold of V _HT > 2000 km s⁻¹ across the entire shock. This value was chosen as it seems to partition nicely the areas of elevated V _HT, and corresponds to roughly 1.3 times the simplified 1500 km s⁻¹ shock speed used to calculate the V _HT values. Since V _HT is inversely proportional to $\cos ({{\rm{\Theta }}}_{{Bn}})$, adopting a threshold of 1.3 times the shock speed means that $1/\cos ({{\rm{\Theta }}}_{{Bn}})\gt 1.3\Longrightarrow {{\rm{\Theta }}}_{{Bn}}\gt 40^\circ$, which will cut out most of the low Θ_Bn data, but still be inclusive of elevated regions of Θ_Bn , not just including the very highest. The in situ shock analysis in Figure 3 of Pulupa et al. (2010) shows that the cumulative distribution function for shocks that exhibit intense upstream Langmuir wave activity only really starts going past 0.1 after Θ_Bn = 40°. These synthetic spectra show a broad overlap with the stencils, suggesting that the process creating the type II burst may be nonlocalized, with multiple points all across the shock passing some physical thresholds, such as V _HT > 2000 km s⁻¹, to start the burst generation process. However, analyzing the similarity distribution as a function of geometry may reveal that the emission is concentrated in distinct regions of the shock. The following figures analyze these trends in the distribution of similarity scores in the synthetic spectra taking both geometry and a threshold of V _HT > 2000 km s⁻¹ into account.

Zoom In Zoom Out Reset image size

One may zoom in on that area of the shock that does best in stationary burst hypotheses, and see if it holds up under various subselections, based on thresholds of plasma parameters. This is done with a cutoff of V _HT > 2000 km s⁻¹ in Figure 15, where the left panel shows the synthetic spectrum for the 3 × 3 degree sector of the simulation around (−30°, 0°) that matches the area of high V _HT in the first hour of the simulation, and the right panel shows the synthetic spectrum for the 3 × 3 degree sector of the simulation around (6°, −12°) with a high degree of similarity to the Wind data. On the left, one can see that the eastern edge of the shock matches the stencil quite well, before tapering off after about an hour. On the right, one can see that the similarity score around the (6°, −12°) edge is nearly unchanged from the left panel of Figure 10, lowering from 82 to 67, only a 20% dip, indicating that a population of points within this geometrical region of the shock with high V _HT could be where at least part of the burst generation occurred.

Zoom In Zoom Out Reset image size

One can also conduct a geometrical analysis of similarity scores while also including the variable threshold subselection V _HT > 2000 km s⁻¹, the results of which are shown in Figure 16. The top two panels show similarity score 2D spatial color contour graphs calculated using only the first hour of simulation data, and are normalized as such. The top left panel shows the distribution using the cutoffs r > 0.9*maxr and V _HT > 2000 km s⁻¹, while the top right panel shows the distribution using the cutoffs r > 0.7*maxr and V _HT > 2000 km s⁻¹. The bottom two panels show similarity score contour graphs calculated using the same respective cutoffs, but using all 2 hr of simulation data, and are normalized as such. So, when comparing the scores of the top row panels and the bottom row panels, one must take roughly half of the upper scores to compare them on the same scale with the bottom scores.

Zoom In Zoom Out Reset image size

From these contour graphs, one may surmise a few things. First, the eastern edge of the shock around (−30°, 0°), which coincides with a large swath of the high V _HT, is the leading edge of the shock for the first hour of the event, and matches the observed burst envelope quite closely in this period (see Figure 15(a)). By looking at the two left panels, one can see that after 1 hr, the eastern edge of the shock does not add any additional points to its similarity score, indicating that after that an hour, not much type II emission originates from there. This can also be clearly seen in the synthetic spectrum for this area of the shock in the left panel of Figure 15. This area of the shock also has high Θ_Bn and V _HT values for the first hour of the simulation, as seen in the variable 2D spatial color contour graphs in Figure 13(a). This eastern edge of the shock can also been seen as a feature in the similarity score histograms, being most prominent in the r > 0.9maxr subselections. One can also see an obvious maximum around (6°, −12°) in the bottom two panels, and a feature in that area in the top right panel, coinciding with the previous similarity score contour graphs around the southwestern edge of the CME shock. This corresponds to the areas where there is a strip of consistently high Θ_Bn and V _HT on the southern and western edges of the shock.

A picture emerges where the key emitting parts of the shock are likely those with high V _HT, with the eastern edge of the shock in the first 40 minutes or so dominating the emission, then falling off, due to not being shocked strongly enough past a certain point or some change in upstream magnetic field decreasing the local V _HT. The western and southern edges of the shock are also seen to have points with consistently high V _HT, with the maximum similarity centered around (6°, −12°), roughly the southwestern corner of the main shock. This high V _HT signature along much of the western edge of the shock persists through all but the beginning of the 2 hr of simulation time, achieving a high similarity score across a wide area in longitude and latitude before the V _HT-based subselection. After the V _HT-based subselection, the maximum around (6°, −12°) is more pronounced, but the area around it remains relatively high, indicating that there may be some level of type II emission from much of the shock or the shock edge where areas of high V _HT exist.

MHD simulations have many limitations and simplifying assumptions, and this one is no exception. A significant source of error arises in how the flux rope is inserted into the simulation to recreate a solar eruption and CME. The CME is initiated by inserting a Gibson–Low magnetic flux rope (Gibson & Low 1998) into the corona, which is line-tied to the bottom boundary at the location of active region 10759 at 16:03 UT. The specific deformation of the flux rope (a mathematical stretching to draw the flux rope toward its origin, while preserving its angular size) leaves the flux rope in a state of force imbalance, causing an instantaneous eruption. This treatment of the flux rope fails to capture the energy buildup and early acceleration phase of the eruption. There is also the matter of magnetic reconnection, which is achieved in the model through numerical diffusion, but which is likely leading to reconnection rates faster than what would actually happen. More recent simulations by Török et al. (2018) capture the liftoff phase. Nonetheless, the important aspects that the simulation needs to capture for this study are the speed profile of the CME and the "upstream" mass density within 20 solar radii to calculate the plasma frequency. As discussed in Section 6, the initial speed of the simulated CME is just above 2500 km s⁻¹, close to the estimate of Reiner et al. (2007), which used the LASCO observations as an initial plane of sky estimate of the CME speed, before fitting the kinematic profile with the radio data. In some sense, the best test of the simulation's accuracy of the density structure close to the Sun is how well the synthetic spectrum generated by the shock, or whatever the emitting structure is, matches the real data. As is seen in the previous sections, the synthetic spectra from simulations of the entropy-derived shock largely match the observed radio data, lending confidence to the simulation's accuracy near the Sun.

There is also an argument to be made that it is not appropriate to compare a synthetic spectrum, which is really a column-normalized 2D histogram of the plasma frequency of points in a selected feature over the simulation domain, to an actual radio spectrum, which has physical units of spectral flux density. This would be a legitimate argument if this framework were attempting to deduce the exact physical mechanism behind particle acceleration and type II burst generation and propagation all at once in this single work. This is not the case. This work is attempting to provide physical thresholds of plasma parameters that, when surpassed, seem to create spectra similar to those actually observed. If these thresholds are consistent across simulations and observations of multiple events, that will be a significant takeaway that can inform forecasting models of earthbound energetic particles. It would also provide heuristic checks that can be performed on fully fleshed-out theories of particle acceleration and type II burst generation, potentially ruling out entire classes of theories.

There is also the issue that the model being used includes an adaptive mesh refinement scheme. This implies that each point being counted may be of a different size. In this case, the blocks are refined in a spherical wedge, centered on the trajectory of the CME, with the blocks closer to the Sun being more highly refined. We argue that this fact is not an issue for the purpose of type II burst template-matching. Many of the points identified in the initial phase of the eruption come from areas of high entropy that may be more highly refined than regions that are more weakly shocked. In general, the relatively small frequency bandwidth observed from type II bursts suggests that the emission source is localized or occurs in a region of uniform density and thus uniform radial distance from the Sun. It is therefore likely that the model emission sources that are successful in recreating the observed radio spectrum will come from localized areas of similar density and plasma frequency that are likely to be of the same level of refinement. Thus, for the proposed purposes of identifying useful physical thresholds or geometrical criteria that correspond to type II burst generation, the apples to oranges comparison between synthetic spectra and spacecraft-observed spectra is not a fatal flaw.

In this work, the simulated structures of the 2005 May 13 radio-loud ICME were compared with the Wind/WAVES recorded radio spectrum to identify the likely regions in which the type II burst was generated. This work attempts to solve the inverse problem of the location of the burst generation site by analyzing in detail the plasma parameters from an MHD simulation of the real event. Different physical structures (shocks and current sheets) of the MHD simulations were extracted to create time-varying synthetic radio spectra, and these were compared to the actual Wind/WAVES radio data in order to determine which hypothesized source region best matches the spacecraft-observed data.

An event-specific "stencil" was developed and used to assign a similarity score to the synthetic spectra, based on how much of their data aligned with the observed Wind/WAVES data. The stencil is an estimated spectral time profile employing a Gaussian intensity profile fitted over the identified type II fundamental emission. This fundamental emission was identified using a kinematic profile of the CME made by Reiner et al. (2007). The larger the portion of data from a given MHD feature whose upstream plasma frequency lies within this stencil, the higher similarity score it will receive. This scoring stencil is described in depth in Section 4.

This methodology bypasses the need for simulating the exact physics that are creating the bursts, and it outputs which areas of the CME are most likely to be the site of particle acceleration. These locations were then inspected more closely, with the additional subselections to the simulation data, to see which plasma parameters are important for type II burst generation to occur. This methodology combines remotely sensed type II radio data with virtual in situ measurements through simulation data, enabling a detailed look at the plasma parameters at the best fit for the site of the emission.

This technique was employed with several different model extractions to test the feasibility of different classes of particle acceleration mechanisms. The four main features used here were: points with at least a 3.5 MK temperature used to identify the current sheet; points that were shocked to a factor of 3.5 or more in their density compared to before the shock wave; points that had their entropy enhanced by a factor of 4 or more compared to before the event; points that had a fast magnetosonic Mach # ≥ 3. Each of these data sources had 2 hr of simulation data, with a time cadence of 2 minutes at the beginning of the run, transitioning to every 5 minutes later in the run, for a total of 35 snapshots over 120 minutes. The one exception to this simulation sampling cadence is the fast magnetosonic Mach # ≥ 3 subselection, which, due to a lack of resources, only combines the six times of [10, 20, 30, 40, 60, 90] minutes after the start of the simulation. The synthetic spectra showing the upstream plasma frequencies of these four broad MHD features alongside the spacecraft-observed radio data reveal that, of the four broad synthetic spectra, the one from the entropy-enhanced feature best matches the Wind/WAVES observations, both in spectral width and overall shape. This similarity is reflected in the computed similarity scores of each of these synthetic spectra, with the entropy-enhanced feature having the highest similarity score.

With these broad MHD features in hand, finer subselections were also applied, based on the plasma parameters thought to be important for SEP acceleration, starting with subselections based on the geometry of the CME. The framework created here explores the optimization problem with a stationary burst assumption, one where the burst location has a constant longitude and latitude, and continues outward with the outer edge of the entropy-derived shock. An analysis of the similarity scores for the different geometrical simulation features from the entropy-enhanced data over the 2 hr of simulation time indicates that the area of the shock near (6°, −12°) has the greatest overall similarity with the main band of fundamental emission in the stencil identified by the kinematic solution of the CME.

In addition to investigating the similarity of the synthetic spectra generated from different geometrical regions across the entropy-derived shock, further constraints on the data are considered, namely a threshold based on the de Hoffmann–Teller speed ( V _HT). This plasma parameter is chosen both for its high correlation with in situ observations of Langmuir waves (Pulupa et al. 2010), and also for the fact that the distribution of this parameter across the surface of the shock generally matches the elevated features seen in the similarity score distributions around the eastern and southwestern edges of the shock, more so than other parameters like upstream velocity magnitude. The spatial distribution of the similarity scores following an additional constraint of V _HT > 2000 km s⁻¹ was also examined, where it is found that the key areas of high similarity are unchanged, with the matching peaks lending confidence to the suitability of using enhancements to the de Hoffmann–Teller speed as a predicting factor of particle acceleration at CME-driven shocks. Overall, the areas with the highest similarity in their generated synthetic spectra coincide with the areas of high Θ_Bn and V _HT on the southern and western edges of the shock, as well as the eastern edge for the first hour. This shared structure was not present in the other plots of plasma parameters across the entropy-derived shock structure. The fact that the structure of de Hoffmann–Teller speed has the best match with the similarity structure of the synthetic spectra across the shock is not particularly surprising, since Pulupa et al. (2010) found that this same parameter was most highly correlated with in situ observations of Langmuir waves at 1 au.

In the future, additional studies utilizing this analysis framework across other historical events could be used to find common physical thresholds that create synthetic spectra matching the observed spectrum. If these thresholds are consistent across simulations and observations of multiple events, that will be a significant takeaway that can inform forecasting models of earthbound energetic particles. It would also provide heuristic checks that can be performed on fully fleshed-out theories of particle acceleration and type II burst generation, potentially ruling out entire classes of theories.

This analysis framework would also become more powerful with the addition of multiple spacecraft for imaging or localization at these low radio frequencies. Single antenna observations can measure the total flux density of solar radio bursts, but cannot pinpoint where the emission is coming from. Previous efforts, such as that by Krupar et al. (2014), have been made at localizing the low-frequency emission from type III radio bursts using the two STEREO spacecraft, applying goniopolarimetric analysis techniques to deduce a wavevector direction and apparent source size at each spacecraft. Combining the measurements from two or more spacecraft allows the approximate triangulation of the source. However, these direction-finding techniques are not a substitute for the true coherent imaging that is enabled at low frequencies by interferometry, as imaging yields information about the sky brightness pattern up to the resolution of the array. The imaging approach will be taken by SunRISE (Alibay et al. 2018; Hegedus et al. 2019), which will utilize six spacecraft equipped with dual-polarization low-frequency radio antennas to form a radio interferometer in space. If reliable localizations of the source of the type II emission are available from SunRISE, or other means, this MHD analysis framework will allow the probing of the plasma parameters at the observed location of the emission, enabling in-depth studies of the physics of acceleration. At least one future paper will be dedicated to the combination of this MHD analysis framework with a simulated SunRISE data analysis pipeline, to prove the feasibility of such an array, and hopefully there will be more studies with actual data after SunRISE launches in 2023. The combination of observations of multiple type II events with SunRISE, alongside simulated MHD recreations utilizing this framework, will hopefully lead to a satisfying answer to the question of how CMEs generate type II bursts and accelerate particles.

A. Hegedus, W. Manchester, and J. Kasper are supported by NASA funds for SunRISE 80GSFC18C0009. A. Hegedus and J. Kasper are also supported by the NASA Solar System Exploration Research Virtual Institute cooperative agreement number 80ARC017M0006, as part of the Network for Exploration and Space Science (NESS) team out of University of Colorado Boulder 1555631. J. Kasper is additionally supported by the NASA Wind fund, 80NSSC18K0986. W. Manchester is also supported by NSF PRE-EVENTS grant No. 1663800 and by NASA grants NNX16AL12G and 80NSSC17K0686. The CME simulation was performed with the Pleiades system provided by NASA's High-End Computing Program, under award SMD-11-2364, and the NCAR Yellowstone supercomputer. This research made use of NumPy (Harris et al. 2020) and matplotlib, a Python library for publication-quality graphics (Hunter 2007). The authors would like to thank the many individuals and institutions who contributed to making Wind/WAVES possible. The authors would also like to thank the reviewers for their helpful suggestions.