M K Crosley

R A Osten

J W Broderick

Stéphane Corbel

J Eislöffel

Jean-Mathias Griessmeier

J A Van Leeuwen

A Rowlinson

Philippe Zarka

C Norman

J.-M Grießmeier

Anton Pannekoek

The Search for Signatures Of Transient Mass Loss in Active Stars

Keywords: stars: coronae, stars: flare, methods: observational

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

M dwarfs are the most common type of star (Henry et al. 2006). They have been shown to have a high rate of hosting exoplanets. Johnson et al. (2007) found that the rate of M stars hosting Jovian planets within 2.5 AU is 1.8% ± 1.0% for M dwarfs. Dressing & Charbonneau (2013) claim, with a 95% confidence lower limit, that the occurrence rate of earth sized planets in the habitable zones of cool stars is 0.04 planets per star.

In studying the potential habitability of these exoplanets, phenomena such as flaring and the potentially more impactful Coronal Mass Ejection (CME) have to be taken into consideration. A CME impacting a planetary atmosphere can cause space weather and other deleterious affects (Lammer et al. 2007;Khodachenko et al. 2007). Due to tidal locking, the rotation rate of exoplanets residing in the habitable zone of an M dwarf is expected to be low (Kite et al. 2011). The weak magnetic moment of these exoplanets produce a weakened magnetosphere because of the low rotation rate (Khodachenko et al. 2007). Therefore they are not expected to be as capable at deflecting CMEs as the Earth (Lammer et al. 2007). CMEs have the potential to erode the atmosphere when the pressure exerted by the CME compresses the magnetosphere into the atmosphere (Khodachenko et al. 2007).

Solar scaling relations have been used for these phenomena in order to relate stellar flares that are observed to the properties of CMEs which might be occurring with the flares (Aarnio et al. 2012;Drake et al. 2013;Osten & Wolk 2015). A few observed stellar counterparts to solar flare phenomena include white light flares (Hawley et al. 1995), nonthermal gyrosynchrotron emission (Osten & Bastian 2006), and coherent radio emission (Osten & Bastian 2006). Although the details are still unclear, this multi-wavelength nature of the flare process suggests that both solar and stellar flares share the same basic physical processes. This motivates a connection between stellar flares and stellar CMEs in the same way as between solar flares and solar CMEs.

Although applying solar scaling relations provides an easy approach for modeling, there is reason to be skeptical. Detailed comparisons of solar and M dwarf flares show a difference in the nature of accelerated particles. M dwarf flares produce more accelerated particles than solar flares. Due to a "quiescent" nonthermal microwave component during times when no obvious flares are occurring and a perpetuated hot plasma, M dwarf flares produce more accelerated particles. This suggests that the flare like mechanism of particle acceleration and coronal heating operate at all times in the corona (Güdel et al. 1996). Also, solar particles do not penetrate the photosphere deeply enough to reproduce the observed magnitude of the white light stellar flare signal (Kowalski et al. 2015) causing difficulties in the application of the standard solar flare model to M dwarfs (Allred et al. 2006).

Despite the fact that stellar flares are routinely observed on different types of stars, clear signatures of stellar CMEs have been less forthcoming. In the past, observations of stellar radio bursts at meter wavelengths were limited to those bursts with intensities greater than several hundred mJy (Bastian 1990). Only a few attempts have been made (Jackson et al. 1990; Kundu & Shevgaonkar 1988;Güdel & Benz 1989;Davis et al. 1978) to search for radio emission. The emission from flare stars at these wavelengths were often assumed to be analogous to the solar type II and type IV radio bursts. Jackson et al. (1990) used the now closed Clark Lake radio telescope to search for radio emission from flare stars in the 30.9 to 110.6 MHz range. No events were detected in the 30.9 to 110.6 MHz range within limits of 1 Jy (for 1 hr integration) to 50 Jy (for 1s integration). Kundu & Shevgaonkar (1988) used the 333 MHz capability at the VLA to detect highly polarized emission from YZ CMi. UV Cet was observed at 333 MHz by Kundu & Shevgaonkar (1988) and Güdel & Benz (1989), but had no detections. Boiko et al. (2012) used the ground-based radio telescope UTR-2 to observe flare stars AD Leo and EV Lac. They claimed, with high confidence, that they found radio emission originating from AD Leo that had burst structure similar to solar type III bursts in the decameter range. A stellar analogue of the solar type II burst has not been detected (Leitzinger et al. 2011).

Traditional solar observations of CMEs use coronagraphs to observe Thomson scattering of photospheric photons off coronal electrons (Vourlidas & Howard 2006). This white light emission is faint when compared to the integrated solar disk emission necessitating the use of a coronograph. Current astronomical coronagraphs cannot achieve sufficient star contrast to detect a CME at a distance of 1 -2 R * [START_REF] Mawet | Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave[END_REF]) making it unfeasible for CME observations. Therefore, only integrated stellar disk emission can be used for observations. Flare observations are not affected by this constraint, but the ability to use Thomson scattering to determine CME properties such as its size, mass, speed, and occurrence rate has been removed.

The type II burst, non-thermal radio emission originating from fast mode magnetohydrodynamic (MHD) shocks, is an alternative diagnostic (Gopalswamy 2006). The shocks accelerate non-thermal electrons, which in turn produce radio emission at the fundamental and harmonic of the local plasma frequency via well-known plasma processes (Gopalswamy 2006). Type II bursts appear as a slowly drifting radio burst following an exponential path through frequency with time related to the speed of the shock and the density scale height in the corona (Gopalswamy 2006).

Typical solar emission frequency drift rates of type II bursts have a range based on starting frequency, shock speed, and characteristics of the corona. Starting frequencies of 100, 10, and 1 MHz, result in frequency drift rates of ∼ 10 -1 , 10 -3 , and 10 -5 MHz/s respectively (Gopalswamy 2006). The duration in which the type II burst drift from their starting frequency (∼ 100MHz) down to frequencies below the ionospheric cutoff is many minutes (Gopalswamy 2006).

Typical solar CME speeds range from a few 100 km/s to about 3000 km/s (Aarnio et al. 2011). The probability of having an associated type II burst increases with CME speed. The association rate reaches 100% at speeds near 2000 km/s (Gopalswamy et al. 2008). Although not all CMEs produce type II bursts, a type II burst is a direct indication that a CME has occurred (Gopalswamy 2006). Thus, the type II rate can be use as a lower limit on the CME rate.

Given the coronal characteristics of M dwarfs, these shocks are expected to be observable in the frequency ranges between Earth's atmospheric cut off of ∼10 MHz to ∼1 GHz. The ∼1 GHz limit corresponds to emission from the base of the stellar corona where electron densities of ∼ 10 10 cm -3 occur (Ness et al. 2004). LOFAR is a new facility that is capable of exploring part of this frequency space for detections of transient mass loss.

LOFAR is capable of exploring the frequencies between 10 -190 MHz. LOFAR's millisecond to microsecond time resolution (van Haarlem et al. 2013) is a valuable asset for looking for variable-timescale transient objects such as a type II burst. This paper presents the search for bursting low frequency stellar radio emission and a methodology for interpreting radio observations in concert with flare recordings. In Section 2 properties of the star YZ CMi are discussed and the details of LOFAR's beam-formed mode are presented. The methodology for data reduction is shared in Section 3. The results of the observations are discussed in Section 4.1. Physical properties derived from the results are examined in Section 4.2. Section 4.3 presents the methodology for estimating frequency drift rates of observed type II bursts from the physical properties found in Section 4.2. An examination of the assumptions used in analysis is shared in Section 4.3. A methodology for connecting flares and CMEs to further improve analysis for future experimentation is proposed in Section 5. Section 6 concludes.

Observations

The primary observing target was YZ CMI, an active M dwarf, because of its high optical flare rate and previous flare observations allowing for detailed knowledge of its stellar parameters. It has a radius of ∼0.36 R [START_REF] Mullan | [END_REF], mass of ∼0.34M (Lim et al. 1987), and is at a distance of 5.93 pc (Perryman et al. 1997) from Earth. The magnetic and plasma properties have been determined from previous observations: the magnetic field strength at flare sites between 50-100G (Raassen et al. 2007), a coronal temperature of 10 7 K (Dupree et al. 1993), and a best value of quiescent state coronal density of ∼ 3 × 10 10 cm -3 (Ness et al. 2004). YZ CMi has an optical flare rate of 0.4 flares/hr for flares above U-band energy of 5×10 31 ergs (U-band is centered around 364 nm) (Lacy et al. 1976). Also, previous observations of YZ CMi have detected low frequency (1500 and 300 MHz) microwave events (Kundu & Shevgaonkar 1988).

By observing YZ CMi at radio frequencies, the presence of CMEs can be confirmed and their properties may be constrained via detections of type II radio bursts. LOFAR's tied-array Beam Formed Mode was used for the pilot project LC0 013 for this purpose.

Instead of producing interferometric visibilities, LOFAR's beam-formed modes can combine the LOFAR collecting area into "array beams", either the coherent or incoherent sum of multiple station beams (Stappers et al. 2011). These data are used to produce time-series and dynamic spectra for high-time-resolution studies (intensities as a function of time and frequency). The Coherent Stokes sub-mode produces a coherent sum of multiple stations ("tied-array" beam) by correcting for the geometric and instrumental time and phase delays (van Haarlem et al. 2013). This produces a beam with restricted field-of-view (∼5 deg maximum), but with the full, cumulative sensitivity of the combined stations (van Haarlem et al. 2013). These observations used LOFAR's core observation stations and did not utilize any of the remote or international stations. Between 20 and 24 core stations The five observations presented in this paper took place on February, October, and November of 2013 and started at 03:45:00 UT on their respective days. The 15 hours of observations were split between either HBA or LBA modes in three hour intervals as listed in Table 1.

Observations include simultaneous on-beam (science target) and off-beam (an adjacent blank patch of sky) signals. Off beams had a declination offset of 10 degrees for HBA and 15 degrees for LBA relative to the on-beam. The frequency channels were averaged in frequency to create 244 subbands. HBA observations were averaged over 64 channels while LBA were averaged over 256 channels to create a subband. This creates the same effective bandwidth per subband of 195 kHz for both HBA and LBA observations. See Table 1 for a comprehensive summary of the observations.

The sensitivity of a combined set of coherent tied-array beams is equivalent to that obtained from the sum of the collecting areas of all the stations being combined (Stappers et al. 2011). By using the sensitivity of N equal stations and the System Equivalent Flux Densities (∆S) listed on LOFAR's sensitivity web page1 with a integration time (δt) and bandwidth (δν), a theoretical sensitivity can be determined ∆S =

Ssys √ N (N -1)2δtδν
. These observations did not include a calibration target for absolute calibration so theoretical limits will be imposed as upper bounds for sensitivity.

Data Reduction

In order to reduce the data, an average of frequency channels for both LBA and HBA are taken. Frequency is averaged, as explained above, and the time scale is averaged from the 0.131s integration time used in observations to ∼1s timescales or longer. To account for any frequency dependent signal intensity variations, each frequency subband (244 subbands in total) is scaled by the median value of that subband providing a scaled unit-less intensity.

Figure 1 shows the data after initial averaging in time and frequency. A variety of transient and pseudo-persistent effects plague the observed spectra. One such aspect, not displayed here, is the ionospheric frequency cutoff (∼10 MHz) affecting the lowest observable frequencies of LBA spectra. Any prominent characteristic present in both beams is a sign that the signal source is inherently terrestrial as it is independent of pointing and is considered to be Radio Frequency Interference (RFI). Therefore any shared characteristics between the on-beam and off-beam shall be flagged before proceeding to the next step (and after averaging has occurred). Note. -The bandwidth used for calculations is 195 kHz and the integration time was 60 second.

1 The order of the sensitivities listed matches the frequencies listed in Table 1 for each Observation ID.

To rigorously identify statistically-significant events throughout the data, a frequencydependent noise level is determined for the data. Figure 2 shows graphically the start of this process. The first step is to model how the sensitivity in each frequency channel is changing in time to help remove systematic variations. A 3 rd degree polynomial is used to generate the baseline. The deviations from this baseline are then calculated and plotted in Figure 3. Using RMS variations from this fit, a Gaussian distribution can be modeled to determine outliers and determine frequency-dependent noise level (σ).

After repeating for both on and off-beams, the two are then compared to look for potential signals.

Results and Discussion

Results

These observations have yielded no plausible detections of drifting low-frequency (type II) radio bursts on YZ CMi and thus no indirect CME detections either. A type II burst is the result of a CME creating an MHD shockwave as it travels though the stellar atmosphere. Since not all CMEs create type II bursts, the frequency of type II's will always be less than or equal to the frequency of CMEs. Therefore the frequency of type IIs is a lower bound for the frequency of CMEs. Assuming the rate of shocks is a lower limit to the rate at which CMEs occur, no detections in a total of 15 hours of observation places a limit of ν typeII < 0.0667 shocks/hr ≤ ν CM E for YZ CMi due to the stochastic nature and of the events and limits of observational sensitivity.

Section 5 proposes a methodology for performing simultaneous observations of flares (optical observations of a star) while observing for type II shocks (radio observations of a star) in the manner presented in this paper. The results of this methodology is to advance and improve further experiments by providing a more in-depth list of constraints on CME

properties than what can be achieved by radio observations alone.

Limits

Yashiro et al. (2006) showed that X-ray solar flare energies above 10 31 erg in the GOES bandpass (1-8 Å) have a 100% correspondence between solar flares and CMEs. Adjusting the flare energy limit placed by Yashiro et al. (2006) to be representative of u-band energy (Osten & Wolk 2015), this would imply the limit would be ∼ 2 × 10 31 erg. YZ CMi has an optical flare rate of 0.4 flares/hr for flares above U-band energy of 5×10 31 (Lacy et al. 1976).

Assuming that the solar relationship holds for M dwarfs, this would imply that YZ CMi would have an expected CME rate, while not accounting for observational constraints, of ν CM E ≥ 0.4 CMEs/hr. This is much greater than the upper limit 0.0667 shocks/hr ≤ ν CM E .

An expectation of 1.2 optical flares per 3 hour observation session leads to a predicted 1.2 shocks per 3 hr observation period for YZ CMi. The Poisson statistics state that the probability of getting 0 detections is ∼ 0.25%. This is suggestive of a more definitive result as opposed to an unfortunate low-activity period for YZ CMi during observations.

With the lack of confirmed detections of stellar type II bursts, the sensitivity of LOFAR shall be used to determine the limits on (T b ) and source size (Ω) of the Type II bursts.

The flux density can be estimated from the type II burst brightness temperature and the size of the emitting region. The flux density (S ν ) is proportional to the brightness temperature (T b ) and solid angle (Ω) (S ν ∝ T b ν 2 dΩ). The brightness temperature T b can be estimated from models of the waves involved (Benz & Thejappa 1988) and is about 10 14 K for solar type II bursts. Stellar coherent bursts show evidence for larger T b (Osten & Bastian 2006); therefore the stellar case is taken to be T b ≥ 10 14 K.

To have a better understanding of where expected Type II bursts will produce detectable results, Figure 5 displays flux densities that are a factor of 5 above their respective noise levels at LOFAR for various frequencies. Table 2 is used for the frequency-dependant noise levels used in calculations for Figure 5. In order to have detectable bursts, large angular scale, and high brightness temperature CMEs are required. Very bright bursts (10 16-17 K) may still be undetectable if they do not fill a large enough fraction of the stellar surface area. This is likely a contributor to the lack of detections as even a CME that produced a type II shock in the correct parameter space would have the additional constraint of having high enough brightness temperature and/or large enough angular size to be detected by LOFAR. the CME is visible), Figure 5 shows that at roughly 1.6 times the stellar surface area, the minimum brightness temperature of ∼ 5 × 10 13 K for higher would be required for LOFAR's observable frequencies. This is near the T b ≥ 10 14 K limit set previously for expected brightness temperatures.

Estimation of Drift Rates

The characteristic shape of a Type II burst is a consequence of the emission mechanism for this event. Langmuir waves are produced by electrons which are accelerated by the MHD shock generated from the movement of the CME through the corona and radiate via the local plasma frequency and its harmonics. As the shock moves through the corona, the frequency of the accompanying radiation also varies. The frequency will vary as:

dν dt = ∂ν ∂n ∂n ∂h ∂h ∂s ∂s ∂t ( 1
)
where ν is frequency, n is electron density, h is radial height above the star, s is distance along the path which the shock travels, and t is time. The drift rate is thus composed of four terms describing the changing environment around the shock. As the shock propagates away or towards the star, the local electron density (n) will change. The shock will emit via the local plasma frequency (ν p = ne 2 0 me in SI) where n is electron density, m e is the electron mass, e is the electron charge, and 0 is the permittivity of free space. An outward moving shock will emit at a lower frequency as it travels towards lower densities

ν p ∝ √ n, ∂ν ∂n = ν 2n .
The second term describes how the density changes as a function of height radially above the stellar corona. The start frequency gives information on height in the stellar atmosphere via a barometric model ∂n ∂h = -n H 0 . The final two terms describe the path and speed the CME takes as it travels. The CME does not necessarily travel perpendicularly outward and so ∂h ∂s = cosθ describes how the vertical height changes as a function of the path traveled. Here θ is the angle at which the shock is traveling relative to the radial direction. The CME will travel a sufficiently small distance over the course of our observations such that we may assume that it will have a constant speed v s . Taking this all into consideration leads to a simplified expression for the drift rate of:

dν dt = ν 2n - n H 0 (cosθ) (v s ) = - νv s cosθ 2H 0 ( 2
)
where ν is frequency, v s is velocity of the shock, θ is the angle at which the shock path is oriented relative to the radial direction, H 0 is the density scale height.

The emitted frequency from the type II burst is an observable and therefore can be directly measured. Measuring how the observed frequency changes with time provides a direct measurement for the drift rate. Density scale height can be estimated from coronal base temperatures, leaving v s cosθ as the remaining unknown parameter. Assuming a perpendicular path (cosθ ≈ 1), a detection of a drifting radio burst provides a constraint on the minimum speed of the CME shock and therefore the CME itself.

In order to evaluate drift rates, the known stellar parameters of YZ CMi stated above will be used for calculations. The parameters are: coronal density scale height will be 1.67 × 10 3 T cm with 10 7 K (Dupree et al. 1993), radius R * =0.36R [START_REF] Mullan | [END_REF]), mass 0.35M , magnetic field strength at base of the corona B 0 = 100 G (Raassen et al. 2007), and n 0 = 3 × 10 10 cm -3 (Ness et al. 2004) for number density. The magnetic scale height H B for an exponential field is assumed to be comparable to R for the low-order magnetic field. The Alfvén v A = 2.03 × 10 11 B √ n speed is used as the minimum speed required for a shock to be generated. An Alfvén speed scale height v a ∝ v 0 e R/R 0 can be estimated using an effective scale height equal to 2HnH b 2Hn-H b = 9.2 × 10 10 cm. Taking account of all stated values, Table 3 displays various values of drift rates and initial emission frequencies for various shock formation distances.

Figure 6 shows Type II bursts generated in various initial stellar conditions (densities and magnetic field strengths at base of corona) that could generate a signal in the frequency range for the LBA observation (30-78 MHz) using plausible coronal conditions for YZ CMi.

The bursts shown are formed such that their start frequency is 78 MHz. Each shock will form at the same local density but the radial height that will correspond to this density will depend on the initial base coronal density. Since this shock is required to be super-Alfvénic, the magnetic field (B) will be used to determine the minimum speed the radio burst can travel. The red and the blue lines on this figure are also used as a rough guide for what the true shape of a type II burst is expected to be for YZ CMi. This is intended to show the longest possible duration of a type II burst under various stellar conditions. The duration of the type II burst can span from ∼12 mins to over 100 mins based on the stellar conditions alone. A measured frequency, a frequency drift rate determined from observations, and prior knowledge about the coronal temperature to determine a barometric density scale height can be used with equation 2 to extrapolate the velocity of the shock.

One final factor to be mindful of is the effects of signal dispersion through the Interstellar Medium (ISM). As the signal is propagating through the ISM plasma, a frequency-dependent time delay is introduced. The delay between any two frequencies is given as:

∆τ = e 2 2πcm e 1 ν 2 1 - 1 ν 2 2 L 0 N ( )d ( 3
)
where ∆τ is the time delay between the two frequencies ν 1 and ν 2 , the ISM particle density is N ( ) ≈ 0.03 cm -3 (Shu 1991), the distance L to YZ Cmi is 5.93 pc, c is the speed of light, and the electron charge and mass are e and m e respectively. The HBA observations (115 -163 MHz) will have maximal delays of 0.029 seconds, which is insignificant for this analysis.

The LBA runs 98796 (30 -78 MHz) and 98950 (10-58 MHz) will have maximal delays of 0.7 and 6.7 seconds respectively. Any detections in the lowest frequency observations would require that this effect be accounted for to provide accurate drift rate measurements, however it is a small correction. As it is difficult to see its impact directly, Figure 7 shows the signal delay as a function of frequency to emphasize the impact of the delay on lower frequencies. This curve assumes a starting frequency of 190 MHz for an M dwarf with density at the base of the corona n 0 = 10 10 cm -3 and magnetic field B 0 = 100 G.

Examination of Assumptions

There are three primary assumptions made for this experiment. The first is that the increase in flare rate of an M dwarf relative to the Sun is partnered with an equal increase in rate of CME occurrence. Simply, this is the assumption that solar scaling relations between the Sun and M dwarfs hold. Possibilities such as an overlying magnetic field structure which is preventing the ejection itself, as has been observed for large solar flares lacking a respective CME (Sun et al. 2015), or perhaps the existence of an alternative mode for energy release assuming the role of a CME are not considered.

The second assumption is that a shock will form, given sufficient conditions, as the result of a CME passing through the corona. If overlying fields are preventing the CME from ever reaching super-Alfvénic speeds, or if the densities are being changed in some unknown manner, then no shock will ever form to result in the type II burst required for detection.

Finally, it is assumed that when a CME has occurred that has created a type II burst, that the burst is detectible. As discussed in section 4.2, the signal is dependent on the frequency, brightness temperature T b , and the angular size Ω. Conditions could be such that any CMEs released are only generating MHD shocks further away from the base of the stellar corona resulting in type II bursts at frequencies below the magnetospheric cutoff and would be undetectable using ground based observatories. Another reasonable alternative is that the brightness temperature and/or the angular size are too small to produce a detectable signal. It may be the case that the shocks are all very bright but so angularly small that they still are undetectable, or the reverse case that they cover quite a large area but are too cold to be seen.

Method for Connecting Flares and CMEs

On the Sun, flares and CMEs are closely related phenomena; they both occur as the result of magnetic reconnection and manifest together in eruptive events. Magnetic reconnection occurs when highly stressed and opposing magnetic field lines are pressed together. The field lines break and reconnect into a new configuration allowing the release of energy and/or mass. Yashiro et al. (2006) showed that X-ray solar flare energies above 10 31 erg in the GOES bandpass (1-8 Å) have a 100% correspondence between solar flares and CMEs. Flares occurring on M dwarfs typically have energies between 10 28 to 10 34 erg and have mean flare rates that are 10 4 -10 5 times higher than on the Sun (Lacy et al. 1976;Kowalski et al. 2010). Assuming the relationship holds for M dwarfs, this would imply a large number of M dwarf flares would have the 100% correspondence of flares to CMEs.

What type of characteristics will this associated CME have for a flare of a certain energy? Physically driven arguments show that in the solar case there is a rough equipartition between the total radiated energy of the flare and the mechanical energy of its associated CME (Emslie et al. 2005;Osten & Wolk 2015).

1 2 M CM E v 2 = E rad f ( 4
)
Where f is the fraction of the bolometric flare emission appropriate for the waveband in which the energy of the flare is being measured. In this case, f = f u = 0.11 for a u-band (300-430 nm) optical observation (Osten & Wolk 2015), is a constant of proportionality ≈ 0.3, v is the velocity of the CME, M CM E is the total mass of the CME, and E rad is the radiated optical flare energy in the u-band.

Given a measured flare energy, and this relation, there still exists a mass-velocity degeneracy for a CME occurring with a flare. However, radio measurements of the CME provide a velocity constraint which can be used to break the degeneracy. Assuming this solar relation holds for the stellar case, a method for constraining CME mass, that was previously impossible using radio observations alone, is obtained.

An experimentally found empirical relationship between the solar flare X-ray energy and its associated CME mass (Aarnio et al. 2012;Drake et al. 2013), provides a second method to constrain CME mass.

M CM E = K M E γ GOES [g] (5)
Where M CM E is the CME mass, K M = 2.7 ± 1.2 × 10 -3 in cgs units, γ = 0.63 ± 0.04 (Aarnio et al. 2012), and E GOES is the energy flare energy (erg) in the GOES (1-8 Å) band.

To apply this relation to the stellar case, a conversion from the GOES bandpass to the optical bandpass (u-band) is required,

M CM E = K M f GOES f u γ E γ u [g] (6)
with f u = 0.11 and f GOES = 0.06 (Osten & Wolk 2015) as the fractions of the bolometric flare energy designated to their respective frequency ranges of the optical bandpass and the GOES bandpass of 1-8 Å. As with the previous case, assuming that the solar relation extends to the stellar case, this provides a direct method for constraining the CME mass from observed optical flare energies. There are large uncertainties in estimating mass from flare energies using equations ( 5) or ( 6). The uncertainties in K M and γ together give about 1-2 orders of magnitudes spread. Therefore, its success may be limited when examining events individually, but it should have more so when dealing with larger quantities of detections or while searching for average M CM E .

In future experiments, a solar scaling test is to verify that equations 4 and 6 provide consistent results for M CM E . This test may bolster or refute the assumption that the solar scaling relations extend to the stellar case. Assuming the test holds, equations 4 and 6 may be utilized to solve for flare energy as a function of CME velocity.

E u = f u 2 K M f GOES f u γ v 2 1 1-γ (7)
. This is now a single variable equation of E u (v). Simultaneous observations of radio and optical frequencies are required to evaluate this equation. The radio observations are able to provide drift rates of putative type II bursts produced by the CMEs, and thus velocity constraints. Optical u-band observations are able to measure the associated flare energy.

Equations ( 4), ( 6), and ( 7) provide a theoretical answer to the question of 'What type of characteristics will this associated CME have for a flare of a certain energy?' For this analysis, it is assumed that a detectable shock is caused by a CME traveling with a velocity slightly above the Alfvén speed generating an MHD shock in the traditional manner. This method will provide a theoretical lower bound of flare energy required by a flare to produce a shock as it is just above the threshold of shock formation. Here the Alfvén speed is represented with v a where B is the magnetic field strength and n is the local electron density.

Assuming that the type II burst is slightly super-Alfvénic (1.2v a ) and thus able to create a type II shock, the minimum energy required by the associated flare can be estimated. Zucca et al. (2014) suggests that solar type II shocks preferentially form near the solar surface (∼ 1 -2 R ), and Susino et al. (2015) suggest that the entirety of the shock surface may be super-Alfvénic at early phases (2 -4 R ). Therefore, it will be assumed initially that a shock will form at R = 2R * . Using typical numbers for M dwarfs at the base of the corona of n 0 = 3 × 10 10 cm -3 (Ness et al. 2004) and B 0 = 100 G (Raassen et al. 2007), the minimum flare energy becomes E u ≈ 1.9 × 10 30 erg. This energy represents a moderately sized u-band flare. Table 3's bottom row lists the associated flare energies for type II bursts at various distances with the main trend of requiring less energy to generate shocks at further distances. To compare this energy to that predicted by Yashiro et al. (2006) for where a one-to-one ratio occurs for flares and CMEs on the Sun, this value must be multiplied by f GOES fu or approximately 0.545. This brings the predicted flare energy to E GOES ≈ 1.07 × 10 30 erg which is close to the predictions of Yashiro et al. (2006). The minimum energy to shock is very sensitive to the initial conditions n 0 and B 0 . Lowering the magnetic field strength or increasing the number density all cause the minimum energy to shock to lower and vice versa for an increase.

By rewriting equation ( 2) to v = 2H 0 cosθ ν ν , equation ( 7) can be manipulated into the form:

E u = 2H 2 0 f u K M cos 2 θ f GOES f u γ ν ν 2 1 1-γ (8)
where E u (v) has become E u ν ν . Figure 8 displays equation ( 8) graphically for several specified frequency values assuming that cosθ ≈ 1. Evaluating equation ( 8), with the parameters provided by Aarnio et al. (2012) for K M and γ, leads to the simplification MHz (right). The green/lightgray shaded region represents the parameter space where LOFAR's LBA system would be sensitive and the teal/dark gray region represents the parameter space where the HBA system would be sensitive under various stellar conditions.

The figure uses a magnetic field strength of B 0 = 100 G and n 0 = 3 × 10 10 cm -3 at the base of the corona for the green and teal regions while it uses B 0 = 300 G and n 0 = 10 10 cm -3 for the gray and dark gray regions. For the LBA (green / light gray), the lowest green corners (along the 10 MHz line) correspond to ν = 0.007 MHz/s at energy E u ≈ 5 × 10 27 erg, and The shaded region displays the parameter space that LOFAR is able to cover under various conditions. Light gray and green is used to represent LBA and dark gray and teal is used to represent HBA. Using a lower limit of 1.2v a for the minimum speed required to drive a shock the gray lowest corner ν = 0.046 MHz/s at energy E u ≈ 1.1 × 10 32 erg. For the HBA (teal / dark gray), the lowest teal corners (which would follow a 110 MHz line) correspond to ν = 0.19M Hz/s at energy E u ≈ 5.5 × 10 29 erg, and the dark gray lowest corner ν = 1.2 MH/s at energy E u ≈ 1.2 × 10 34 erg. As an example to read this plot, if a 10 MHz burst was observed to have a frequency drift rate of 3.5 × 10 -2 MHz/s, then the energy of a flare that would have been associated to that drift rate would be about 10 32 ergs.

The bottom edge of each shaded region in Figure 8 is defined as the energy of a flare that would be associated with a shock that is moving at 1.2v a ∝ B √ n . This is just above the minimum speed required a CME to generate an MHD shock. This is a physically based constraint and not a technical limitation for observation. Both increasing the number density or lowering the magnetic field strength will push the left most edge lower down the lines of constant frequency. There is no upper boundary as it is possible to have arbitrarily energetic shocks that would be observable. For clarity this means that the green and teal regions will extend into the gray/dark gray regions and beyond.

Figure 8 is intended to emphasize that the observable flare energies will be constrained to a range of drift rates and can change dramatically based on the conditions of the star.

The green and teal regions represent the theoretical sensitivity for a star similar to YZ CMi, while the gray and dark gray regions represent conditions that differ by 1/3 the density and 3 times the magnetic field strength. Using the flare energy shown by Yashiro et al. (2006) to have a one-to-one correspondence between flares and CMEs (10 31 ergs), the sensitivity to stellar parameters can be shown. A YZ CMi type star with a type II burst emitting in LOFAR's observing frequency range would be fully covering the parameter space (green and teal regions been the lines of constant frequency of 10 and190 MHz for a constant energy of 10 31 ergs).

Alternatively, a star similar to YZ CMi but with 1/3 the base coronal density and 3 times the magnetic field strength with a type II burst emitting in LOFAR's observing frequencies would have a signal that would be completely outside of the parameter space (gray and dark gray regions) that LOFAR would be able to cover. Here the gray and dark gray regions do not extend down to the 10 31 ergs energy region between the lines of constant frequency 10 and 190 MHz and thus would be unable to detect CMEs corresponding to a flare with this energy.

A star with low densities and strong magnetic fields may not produce an observable shock with a flare of certain energy while a star with high density and weak magnetic fields would produce a detectable shock off an identical flare. Preferential stars to observe for shocks are those stars with lower magnetic field strengths and higher coronal electron densities while keeping high flare rates. This predicts detectable drift rates of ∼ 10 -1

MHz/s.

Conclusion

The observations presented in Table 1 have yielded no plausible detections of drifting low-frequency (type II) radio bursts on YZ CMi and thus no CME detections either.

Assuming the rate of shocks is a lower limit to the rate at which CMEs occur, no detections in a total of 15 hours of observation places a limit of ν typeII < 0.0667 shocks/hr ≤ ν CM E for YZ CMi due to the stochastic nature of the events and limits of observational sensitivity.

This rate is for the frequency range between 30-160 MHz, at the sensitivities stated above, assuming the rate of shocks is less than or equal to the rate of CMEs. Due to the lack of detections, and the exploratory nature of this paper, no rigorous constraints can yet be placed on whether flares on M dwarfs have CMEs.

A likely contributor to the lack of detections is that CMEs producing a type II shock in the correct parameter space also need to have a sufficiently high brightness temperature and/or large enough angular size to produce a flux density detectable by LOFAR. More sensitive observations (with more stations) and longer observation periods would increase the likelihood of detections.

Based on Section 5, simultaneous observations between radio and u-band would allow for a deeper investigation of M dwarf CMEs. Stars with weaker magnetic field strengths and higher electron densities of the corona while maintaining high flare rates are preferred for type II burst detections. Weaker flares are related to slower CMEs via equation ( 7).

This means that the observable drift rates are slower due to the shock forming further away from the star. Generating high drift rate shocks require shock formation closer to the star in its higher density regions and require more energy reducing the frequency of occurrence.

Simultaneous observation provides direct flare measurements which would allow for CME property predictions that could be tested against CME detections if found.

Given the assumptions and limitations of this current dataset, rigorous limits can not

  Figure4shows a comparison of on and off-beams over the entire frequency range for a 35 min portion of the data set. The points displayed are only those remaining after a 5σ threshold has been enforced in order to find outliers and a linear drift rate slope has been superimposed on top. For reasons explained below, the type II burst is expected to have a duration of tens of minutes. The time range imposed on Figure4is chosen to be long enough so that the entire duration of an expected type II on YZ CMi would be observed. This time is also long enough such that if the conditions differ from what is expected, a partial signal should still be present which would suggest the need to use a longer time segment.

  Susino et al. (2015) studied a GOES M2.6 class solar flare (medium-large flare) and found that in the early phases (2 -4 R ) the whole shock surface is completely super-Alfvénic, while later on it becomes super-Alfvénic only at the nose.Byrne (2012) describes the angular spread of a CME as a function of the radial height (r) above the Sun with a power-law expansion (∆θ(r) = ∆θ 0 r 0.22 , where ∆θ 0 = 26 • ) . The estimated angular spread of a fully super-Alfvénic CME at 4R is ∼ 35 • . The angular scale of the transient shock at this height (4R ) would be ∼1.6 times the size of the solar disc. Assuming that this relation holds true for M Dwarfs as well, it can be expected that a CME would have an area 1.6 times the stellar surface. Assuming a head on CME (where the entirety of

E

  u ≈ 4.75 × 10 44 ν ν 5.4 . The two solid black lines represent lines of constant frequency 10 MHz (left) and 190

  Fig. 1.-98796 LBA beam comparison to display RFI and potential signal detection. Panel a) shows the "on beam" while panel b) is the "off beam". Each frequency subband is scaled by the median value of that subband providing a scaled dimensionless intensity.

Fig. 2 .

 2 Fig. 2.-LBA data set 98796 is used to display an example of creating a line of best fit for determining inherent frequency-dependant noise levels. The displayed frequency is for ν = 74.90 MHz. Gaps occur from prior RFI flagging.

Fig. 3 .

 3 Fig. 3.-A histogram of Figure 2's deviations from the line of best fit is displayed with its modeled Gaussian fit. The Gaussian is used in order to obtain a noise level for this frequency channel (ν = 74.90 MHz). This particular example has a standard deviation of 0.068 relative intensity.

Fig. 5 .

 5 Fig. 5.-The detectable levels of radio emission at several frequencies covered by the observations for combinations of brightness temperature and radio source extent. The noise parameters in Table 2 are used to determine sensitivities. Blue points represent a frequency of 163 MHz and cyan points represent 115 MHz, and purple points represent 77.6 MHz. The points presented are ones that have a signal-to-noise (S/N) ratio greater than 5 and the size of the point is scaled logarithmically to the S/N ratio. The x-axis ranges from values of 0.04 to 5 times the area of the stellar disk.

Fig. 6 .

 6 Fig. 6.-The frequency-time signal expected for a Type II burst under various parameters.Only the intrinsic motion through the outer stellar atmosphere is plotted as ISM signal dispersion effects are not taken into account as they are negligible. The shocks are calculated for the LBA window between 78 to 30 MHz akin to observation 98796. The shock is formed at a height required to have an initial frequency of 78 MHz (initial shock height will be different for each set of initial conditions) and thus the path described will be the slowest potential path through the frequency range possible.

Fig. 7 .

 7 Fig. 7.-The signal delay as a function of frequency. The delay time is measured relative to the arrival time of 190 MHz. This curve assumes a starting frequency of 190 MHz for an M dwarf with density at the base of the corona of n 0 = 10 10 cm -3 and magnetic field B 0 = 100 G.

Fig. 8 .

 8 Fig. 8.-The associated flare energies as a function of drift rate for a potential type II shock.

(

  cosθ ≈1) in LOFAR's frequency range defines the bottom edge. The solid lines represent lines of constant frequency at 10 (left) and 190 (right) MHz. The dashed lines correspond to the uncertainty of this region deriving from determining CME mass from flare energies. The gray and dark gray shaded regions use values of B 0 = 300 G and n 0 = 10 10 cm -3 . The green and teal regions use the values B 0 = 100 G and n 0 = 3 × 10 10 cm -3 . For clarity the green and teal regions continue into and cover the same areas as the gray and dark gray regions.

  The two dashed lines inFigure (8) are representative of the uncertainty in the parameter space induced by the uncertainties of K M and γ fromAarnio et al. (2012).The large spread inAarnio et al. (2012) leading to equation (5) has the majority of its impact on equation (8) via the uncertainties in the power γ despite the fact that it has a lower fractional-uncertainty δγ γ than K M . This creates the potential for several orders of magnitude uncertainty in either directions about any specific measurement. Since the empirical relations inAarnio et al. (2012) andDrake et al. (2013) provide consistent values for CME mass, even given wide spread, a large number of dual measurements are needed to test accurately the validity of solar-stellar scaling.

Table 1 .

 1 LOFAR Observations and InformationComprehensive table of various aspects of the observations used for the analysis in this paper. The HBA filter is between 110-190 MHz, and the LBA filter is between 10-90 MHz. The pointing is for the target YZ CMi. For the off beam, HBA observations have an offset of 10 degrees in Declination while LBA observations have an offset of 15 degrees in Declination. LBA observations were taking with the LBAOUTER layout. HBA observations used the HBAdual configuration.

	# of (Inner) Total Time Sampling Number of	Core Stations Observed [s] Time [s] Channels	24 10799.81 0.131072 62464	23 10799.81 0.131072 62464	2 × 21 10779.91 0.131072 15616	2 × 20 10779.91 0.131072 15616	2 × 20 10779.91 0.131072 15616
	Start Time	[UT]	2013-02-25 03:45:00	2013-02-26 03:45:00	2013-02-28 03:45:00	2013-11-01 03:45:00	2013-10-01 03:45:00
	Observed Width per	Frequencies [MHz] channel [KHz]	29.980-77.636 0.763	10.448-58.105 0.763	114.939-162.596 3.052	114.939-162.596 3.052	114.939-162.596 3.052
	Observation Antenna	Set ID	LBA 98796	LBA 98950	HBA 99198	HBA 184329	HBA 184330	Note. -

to form one 24 element field

(van Haarlem et al. 2013

). The HBAdual configuration was used for these observations.

Table 1 .

 1 Table 2 lists the approximate sensitivities, as determined by the LOFAR sensitivity web page and number of antennas presented in The values in Table 2 are approximate, and represent an underestimate of the true limits.

Table 2 .

 2 LOFAR Sensitivity

	Observation	Antenna	Sensitivity 1
	ID	Set	[mJy]
	98796	LBAOUTER 327 -474
	98950	LBAOUTER 449 -251
	99198	HBAdual	12 -10
	184329	HBAdual	13 -10
	184330	HBAdual	13 -10

Table 3 .

 3 LOFAR Sensitivity type II burst. E u is the energy required of the associated flare to the CME in order to create a shock at that distance. See text for more details.

	Distance in R *	2	3	4	5	6	7	8
	Drift Rate [MHz/s]	0.76 0.27 0.059 0.035 0.013 0.0046 0.0016
	Initial Frequency [MHz] 347 164	53	37	17	8	4
	E u [10 28 ergs]	198 45.3	5	2.4	0.55	0.13	0.0029
	Note. -This calculation assumes that the shocks initially form at the distance
	specified, therefore assuming the slowest possible shock (and lowest drift rate) which
	would generate a

https://www.astron.nl/radio-observatory/astronomers/lofar-imaging-capabilitiessensitivity/sensitivity-lofar-array/sensiti


Vítor Silva
email: silva@cos.ufrj.br

José Camata
email: camata@nacad.ufrj.br

Daniel De Oliveira
email: danielcmo@ic.uff.br

Alvaro L G A Coutinho
email: alvaro@nacad.ufrj.br

Patrick Valduriez
email: patrick.valduriez@inria.fr

Marta Mattoso

In Situ Data Steering on Sedimentation Simulation with Provenance Data

Keywords: libMesh, in situ data analysis, raw data analysis, dataflow management, provenance. SELECT time_step, x, y, z, d FROM horizontal_line_extraction_1 WHERE f.final_non_linear_residual, s.final_linear_residual, s.final_non_linear_residual

Parallel adaptive mesh refinement and coarsening (AMR) are optimal strategies for tackling large-scale simulations. libMesh is an open-source finite-element library that supports parallel AMR and is used in multiphysics applications. In complex simulation runs, users have to track quantities of interest (residuals, errors estimates, etc.) to control as much as possible the execution. However, this tracking is typically done only after the simulation ends. This paper presents DfAnalyzer, a solution based on provenance data to extract and relate strategic simulation data for online queries. We integrate DfAnalyzer to libMesh and ParaView Catalyst, so that queries on quantities of interest are enhanced by in situ visualization.

I. OVERVIEW

Parallel adaptive mesh refinement and coarsening (AMR) are optimal strategies for tackling large-scale sedimentation simulations. libMesh [START_REF] Kirk | libMesh : a C++ library for parallel adaptive mesh refinement/coarsening simulations[END_REF] is an open-source finite-element library that supports parallel AMR and has been used in several multiphysics applications. As an application built upon the libMesh library, libMesh-sedimentation 1 simulates turbidity currents typically found in geological processes. The sediment transported due to fluid motion is described by a Eulerian framework in which the mathematical models result from the incompressible Navier-Stokes equation (fluid) combined with an advection-dominated transport equation (sediment concentration). libMesh-sedimentation employs a variational multiscale finite element method [START_REF] Guerra | Uncertainty quantification in numerical simulation of particle-laden flows[END_REF] where a staggered approach is used to evolve in time the coupled fluid-sediment equations.

In complex simulation runs, users need in situ data steering to track quantities of interest such as residuals, errors estimates, etc. to control as much as possible the execution. Often users need to analyze raw data from multiple related files generated in different solver steps at runtime. To help on this data steering, in a previous work [START_REF] Silva | Analyzing related raw data files through dataflows[END_REF], we defined a dataflow provenance aware database schema, PROV-Df, as an extension of W3C PROV [START_REF] Groth | W3C PROV -An Overview of the PROV Family of Documents[END_REF], to represent selected domain raw data captured by a Scientific Workflow Management System (SWMS) with provenance data [START_REF] Davidson | Provenance and scientific workflows: challenges and opportunities[END_REF]. In this SWMS approach, the simulation has to be modeled as a workflow, and we use 1 https://bitbucket.org/nacadhpc4e/libmesh-sedimentation the Chiron SWMS [START_REF] Ogasawara | Chiron: A Parallel Engine for Algebraic Scientific Workflows[END_REF] to control the execution while extracting domain data from the files generated by the simulation. Chiron is a non-intrusive solution with respect to the simulation code. However, when using a SWMS to manage the parallel execution of the workflow, it often conflicts with the simulation solver scheduling and data parallelism. Most SWMS have provenance data support, but do not manage workflow activities that invoke highly parallel codes and, consequently, they may jeopardize performance. Therefore, we propose online provenance and dataflow queries for solver libraries or other workflows without using a SWMS.

Different from the SWMS approach, we present ARMFUL2 , a component-based architecture that enables the Analysis of Raw Data from Multiple files. ARMFUL allows for raw data analysis by extracting raw data from files and relating them through dataflow provenance data. Furthermore, ARMFUL data analysis support can be coupled to a simulation code, library, script or even a SWMS. In this paper we instantiate ARMFUL as DfAnalyzer, which is coupled to the libMesh-sedimentation solver to extract and relate strategic simulation data for online queries. The DfAnalyzer approach keeps all the parallel execution control at the simulation solver.

In the same way, users plug visualization components to their simulation code, in DfAnalyzer, we provide for raw data extraction components to be plugged in strategic places of the simulation code, e.g. registering time step control data, solver convergence, AMR behavior, and data dependencies between solver steps. DfAnalyzer can also take advantage of in situ data managed by ParaView Catalyst [START_REF] Ayachit | ParaView Catalyst: Enabling In Situ Data Analysis and Visualization[END_REF]. In this case, DfAnalyzer communicates with ParaView to extract data from Catalyst and relate it to its provenance database, benefiting from data that is already in memory. The resulting provenance database of DfAnalyzer complements the data coverage supported by Catalyst queries. While Catalyst supports query processing on a specific dataset (e.g., XDMF file addressing HDF5 files), DfAnalyzer allows for users to monitor the file flow generation during the simulation execution and to relate data elements from different solver steps (dataflow management). Fig. 1 shows an overview of DfAnalyzer, where the query results, obtained at any time during the simulation, allow for in situ data steering. Queries can show, for example, the appearance of sediments for a specific time step or solver convergence difficulties, which helps to change set-up parameters such as tolerances, at runtime. DfAnalyzer has two main components, named as Provenance Gatherer (PG) and Query Processor (QP), to manage dataflow at physical (i.e., file flow) and logical (i.e., data flow) levels [START_REF] Silva | Analyzing related raw data files through dataflows[END_REF]. PG component captures provenance and domain-specific data from application code, generating a JSON file with the extracted data and their dependencies. Then, all captured provenance and domain-specific data are loaded into the database that follows PROV-Df schema and is managed by the column-store Database Management System MonetDB [START_REF] Boncz | Breaking the memory wall in MonetDB[END_REF]. QP helps users to submit SQL queries to the provenance database and to show results. To complement in situ provenance database analysis with visual information at runtime, ParaView Catalyst is embedded into libMesh-sedimentation. Catalyst can access libMesh data structures, extracting additional data information from solution fields and performing live visualization. Catalyst can provide profiles of the deposited sediments that are the relevant quantities of interest for geologists. Those profiles are generated through filters applied on the original data and are defined in Python scripts at the beginning of the simulation. Comparing this solution to the ad-hoc programs for raw data analysis, our approach only needs to export Python scripts using analyses already developed in ParaView UI, without modifying the libMesh-sedimentation code. Then, Catalyst will perform the parallel raw data analysis and visualization based on those scripts. II. EXPERIMENTAL EVALUATION We ran DfAnalyzer integrated to libMesh-sedimentation and ParaView Catalyst using 1,040 cores in Stampede cluster 3 , at Texas Advanced Computing Center (TACC). This proof-ofconcept experiment executes the sedimentation solver with adaptive mesh refinement and consumes an initial input mesh with 480×80×80 hexahedra. The adaptive mesh refinement is applied every 10 time steps. The simulation took 137.75 minutes running 200 time steps. Catalyst is also invoked every 10 time steps. Considering the total elapsed time, 136.80 minutes correspond to the solver elapsed time (99.31% of the total elapsed time), 0.60 minute to the in situ data analysis (0.44%), and 0.35 minute to the provenance data capture (0.25%) -and manipulates approximately 3.17 GB of data -5.7 MB corresponds to the provenance database, 3.15 GB to raw data files, and 15.0 MB to log files. Even in real-life larger runs, where number of files and their sizes are much larger, 3 https://portal.tacc.utexas.edu/user-guides/stampede execution time overhead of DfAnalyzer will remain negligible compared to the solver time.

Figs. 2 and3 present two analytical queries that users are able to run using DfAnalyzer's provenance database. The query in Fig. 2 monitors the appearance of sediments in the domain bottom layer for a specific time step. Fig. 3 shows how users can monitor the progress of their solver at runtime. For example, in Fig. 3 the query returns solver data in a specific time step that allows for users to steer their simulation based, for instance, on the solver convergence. When solver convergence or AMR difficulties are detected, users can benefit from provenance support to trace simulation data from previous non-linear iterations and time steps or meshes, since our provenance database represents the dataflow while being generated by the sedimentation simulation execution and acting in runtime.

Fig. 1 .

 1 Fig. 1. Overview of our solution with DfAnalyzer.

Fig. 2 .

 2 Fig. 2. Monitoring query: analysis of a horizontal line in the mesh associated to the simulation identifier 1, when the deposited sediment (d) in the bottom layer is different from zero.

Fig. 3 .

 3 Fig. 3. Query for debugging and user steering: analysis of the algorithm output parameters after the convergence of the simulation solver in the fluid and sediments loops in a specific simulation.

https://hpcdb.github.io/armful

ACKNOWLEDGMENTS

This work was partially funded by CNPq, CAPES, FAPERJ and Inria (MUSIC project) and the EU H2020 Programme jointly with MCTI/RNP-Brazil (HPC4E project). Computer time on Stampede is provided by TACC at UT-Austin.