Observational Signatures of Coronal Loop Heating and Cooling
Driven by Footpoint Shuffling
The evolution of a coronal loop is studied by means of numerical simulations of the fully compressible three-dimensional magnetohydrodynamic equations using the HYPERION code. The footpoints of the loop magnetic field are advected by random motions. As a consequence the magnetic field in the loop is energized and develops turbulent nonlinear dynamics characterized by the continuous formation and dissipation of field-aligned current sheets: energy is deposited at small scales where heating occurs. Dissipation is non-uniformly distributed so that only a fraction of the coronal mass and volume gets heated at any time. Temperature and density are highly structured at scales which, in the solar corona, remain observationally unresolved: the plasma of our simulated loop is multi-thermal, where highly dynamical hotter and cooler plasma strands are scattered throughout the loop at sub-observational scales. Numerical simulations of coronal loops of 50000 km length and axial magnetic field intensities ranging from 0.01 to 0.04 Tesla are presented. To connect these simulations to observations we use the computed number densities and temperatures to synthesize the intensities expected in emission lines typically observed with the Extreme ultraviolet Imaging Spectrometer (EIS) on Hinode. These intensities are used to compute differential emission measure distributions using the Monte Carlo Markov Chain code, which are very similar to those derived from observations of solar active regions. We conclude that coronal heating is found to be strongly intermittent in space and time, with only small portions of the coronal loop being heated: in fact, at any given time, most of the corona is cooling down.
Subject headings:magnetohydrodynamics (MHD) — Sun: activity — Sun: corona — Sun: magnetic fields — turbulence
The solar magnetic field has long been recognized as playing a key role in the transport, storage and release of energy from the photosphere to the corona (Gold & Hoyle, 1960; Sturrock & Uchida, 1981). Parker (1972, 1994) proposed that photospheric motions set the coronal magnetic field in “dynamic non-equilibrium”, that leads to the formation of current sheets on fast ideal timescales (Rappazzo & Parker, 2013) where magnetic reconnection releases energy in small impulsive heating events termed “nanoflares” (Parker, 1988). This process has been shown to have the characteristics of magnetically dominated MHD turbulence (Einaudi et al., 1996; Dmitruk & Gómez, 1997; Dmitruk et al., 2003; Rappazzo et al., 2007, 2008), where the out-of-equilibrium magnetic field generates a broad-band small velocity that creates small scales distorting the magnetic islands and pushing field lines together (Rappazzo & Velli, 2011). Similar dynamics are also displayed in cold plasma (Hendrix & van Hoven, 1996) and full MHD simulations (Galsgaard & Nordlund, 1996; Dahlburg et al., 2012).
A first connection to observations has been provided by the statistics of these bursty dissipative events, that have been shown to follow a power-law behavior in total energy, peak dissipation and duration with indices not far from those determined observationally in X-rays (Georgoulis et al., 1998; Dmitruk et al., 1998).
But to constrain any model, advance our understanding of coronal heating and correctly interpret observations it is crucial to study the thermodynamics of such a system. Simulations of entire active regions allow the investigation of the geometric properties of radiative emission and thermodynamical quantities (e.g., temperature, mass flows and average volumetric heating rates, Gudiksen & Nordlund, 2002; Zacharias et al., 2011; Bourdin et al., 2013), but their coarse resolution at scales below energy injection (about the granular scale km), necessary to include an entire active region, do not allow the full development of nonlinear dynamics leading to the formation of strong current sheets where energy is deposited.
Magnetic reconnection is not directly observable in the corona because it has become increasingly clear that the effective heating and particle acceleration occurs at scales of the order of the ion (proton) inertial length , which for an ion density cm becomes m (the proton plasma frequency is , the speed of light, the electron charge, and the proton mass), well below the resolution limits of present instrumentation — to date the highest spatial resolution achieved for direct observations of the corona is approximately km by the Hi-C imager (Cirtain et al., 2013). Additionally for typical active region temperatures K and magnetic field intensities G the ion gyroradius is of same order of magnitude as .
What can be observed directly is radiation. By analyzing the spectral properties of the observed radiation it is possible to infer some of the physical properties of the plasma in the solar upper atmosphere, such as the number density and temperature distribution along the line of sight. Thus comparisons between observations and models must focus on the analysis of the spectral properties of the plasma.
Here we analyze results from the HYPERION compressible MHD code. HYPERION is a parallelized Fourier collocation finite difference code with Runge-Kutta time discretization that solves the compressible MHD equations with parallel thermal conduction and radiation included (Dahlburg et al., 2010, 2012). HYPERION is able to produce temperatures and number densities obtained in a framework where the “heating function” is due only to the resistive and viscous dissipation induced in the corona by the footpoint shuffling. Recent simulations (Dahlburg et al., 2012) have shown that temperature is highly structured at scales below observational resolution in loops whose magnetic field lines are shuffled at their footpoints by random photospheric-mimicking motions: temperature peaks around current sheets forming similarly shaped structures, approximately elongated in the strong guide field direction, surrounded by cooler plasma.
In this paper we use our simulations of resolved loops to return predictions for simulated “observables”, such as the number density and differential emission measure distribution, that can be compared with observations. There has been considerable interest in the temperature distribution observed in coronal loops (e.g., Del Zanna & Mason, 2003; Aschwanden et al., 2007; Schmelz, 2002; Warren et al., 2008, 2012). Many of these studies have found relatively narrow emission measure distributions, and it has been unclear how these observations could be reconciled with theory.
We simulate loops of km length and axial magnetic fields of 100, 200, and 400 G. The resulting temperatures and densities are used to synthesize the emission line intensities that the Extreme Ultraviolet (EUV) Imaging Spectrometer (EIS, Culhane et al. 2007) would observe. These intensities are input into the same analysis software used in many observational studies. For these first calculations we find very good agreement between the emission measure distributions derived from the simulations and the general trends in the distributions derived from data. The distributions are relatively narrow, peak at temperatures between and 6.4, and show very little emission at flare-like temperatures (). The mean temperature in the distribution, along with its width, also rises with increasing field strength, consistent with observations.
2. Formulation of the problem
In this section we describe our extension of the Parker coronal heating model from RMHD to a formulation that includes many more significant physical processes. We first describe our magnetohydrodynamic model in which physical augmentations, such as thermal conduction and optically thin radiation, are contained. Line-tied boundary conditions appropriate to the upper chromosphere are then given. The velocity forcing function at the boundaries is also described. The formulations for the elliptical gravity model, initial temperature and initial number density are also given.
2.1. Governing equations
We model the solar corona as a compressible, dissipative magnetofluid with nonlinear thermal conduction and optically thin radiation losses. The governing equations, written here in dimensionless form, are:
with the solenoidality condition . The system is closed by the equation of state
The non-dimensional variables are defined in the following way: is the number density, is the flow velocity, is the thermal pressure, is the magnetic induction field, is the electric current density, is the plasma temperature, is the viscous stress tensor, is the strain tensor, and is the adiabatic ratio.
To render the equations dimensionless we set characteristic values at the walls of the computational box: a number density , vertical Alfvén speed at the boundaries , the orthogonal box width , and the temperature . Therefore time () is measured in units of the Alfvén time ( – note that this is not the axial loop length transit time.). The parallel thermal conductivity is given by , while the perpendicular thermal conduction is considered negligible and hence is set to zero.
The magnetic resistivity , and shear viscosity are assumed to be constant and uniform, and Stokes relationship is assumed so the bulk viscosity . In our previous paper (Dahlburg et al., 2012) the function that describes the temperature dependence of the radiation was evaluated in the same way as Hildner (1974). Here we use instead the radiation function based on the CHIANTI atomic database (Landi et al., 2012), normalized by its value at the base temperature . The Newton cooling term is described in section 2.2.3.
The important dimensionless numbers are: viscous Lundquist number ( kg is the proton mass), Lundquist number ( Henrys / meter is the magnetic permeability), pressure ratio at the wall, Prandtl number, and , the radiative Prandtl number . is the specific heat at constant volume. The magnetohydrodynamic Froude number () is equal to , where m s is the solar surface gravity.
In what follows we assume normalizing quantities representative of the upper solar chromosphere: m, K, and m. is the only quantity that is varied in the three numerical simulations 0.01 Teslas; 0.02 Teslas and 0.04 Teslas; see Table 1). We set . A loop length of = 50000 km is used in all of the simulations. The normalized time scale of the forcing, , is set to represent a five minute convection time scale. The normalized velocity is m s. This velocity is expressed in dimensionless form as .
2.2. Boundary and initial conditions
We solve the governing equations in a Cartesian domain of size , where is the loop aspect ratio determined by the loop length and the characteristic length (). The system has periodic boundary conditions in and , line-tied boundary conditions at the top and bottom -plates, and it is threaded by a strong guide magnetic field in the -direction.
As explained later in subsection 3.1 we utilize the magnetic vector potential rather than the magnetic induction field. In addition, our implementation of a staggered mesh in is explicated in subsection 3.1 Using the normalizing quantities given above, the dimensionless line-tied boundary conditions which are enforced at the top and bottom walls of the simulation take the following form:
The velocity stream function () is described in section 2.2.1. The magnetic field is expressed as with , where is the vector potential associated with the fluctuating magnetic field. At the top and bottom -plates , and are kept constant at their initial values , and , while the magnetic vector potential is convected by the resulting flows.
2.2.1 Velocity forcing function
At the boundaries we employ a time-dependent forcing function analogous to those used in previous studies (Hendrix & van Hoven, 1996; Einaudi et al., 1996; Einaudi & Velli, 1999), i.e., at the top boundary we evolve a function
and at the bottom boundary we evolve a similar function
in which all wave-numbers with are excited, so that the typical length-scale of the eddies is . and are random numbers chosen such that . Every , the coefficients and are randomly changed alternatively for eddies 1 through 4.
At each timestep a provisional wall velocity is computed from:
To ensure that the kinetic energy at the wall remains constant, we compute
separately at the top and bottom boundaries (these are denoted by and ). To achieve the desired velocity we then have the following stream functions at the top and bottom boundaries:
where . Based on these stream functions, the top boundary velocity is given by:
and the bottom boundary velocity is given by:
2.2.2 Initial condition for dynamical variables
As explained later in subsection 3.1 we utilize the magnetic vector potential rather than the magnetic induction field. The initial values for the momentum and magnetic vector potentlal are given by:
2.2.3 Initial temperature, number density and gravity specification
Here we describe how we initialize the temperature and number density, as well as specify the gravity function. The loop gravity is determined by an elliptical model, with
where is the semi-major and is the semi-minor axis of the ellipse. The elliptical model decouples the loop height from the loop length, since and can be specified independently. The footpoints of the loop are located where , i.e., the loop length is given by .
We impose as initial condition a temperature profile () with the dimensionless temperature 1 at the boundaries and 100 in the center (this corresponds to dimensional values of K and K. Let , then:
The parameter determines the steepness of the temperature profile at the boundaries, as well as the flatness of the temperature profile in the center of the system. We set to ensure a rapid increase of temperature away from the boundaries. The number density can then be solved for in the usual manner, i.e.,
Rearranging this equation, we have:
We solve this numerically with a shooting method. Calculating the number density in this way allows us to consider longer loops. In this paper we choose , consistent with (and a dimensional loop length of 50000 km). Combined with our choices for , this places our loop within the range of what is typically observed in the solar corona. To illustrate this, in Figure 1 we show active region NOAA 11082 with field lines in the range 45000 – 55000 km computed from a potential extrapolation of a Helioseismic and Magnetic Imager (Schou et al., 2012) magnetogram.
The term in equation 3 denotes a Newton cooling function which is enforced close to the boundaries (Dorch & Nordlund, 2001; Bingert and Peter, 2011). In dimensionless form we use at the lower boundary and at the upper boundary. Here is the Newton cooling time and is the Newton cooling height. We use and = 1/4. In dimensional terms this corresponds to times between 0.145 s and 0.58 s for the various magnetic field cases (see Table 1), and a height of 1000 km. The Newton cooling term is only effective over the first few points in at each boundary. Note as well that the radiation function is exponentially decreased in the inverse manner near the boundary to account for the increasing optical thickness of the upper chromosphere.
3. Numerical considerations
3.1. Numerical method
With previous definitions equation (4) and the magnetic field solenoidality condition () can be replaced by the magnetic vector potential equation:
A staggered mesh is employed in the -direction (Schnack et al., 1987). The fields that are defined at the boundaries are advanced in time on the standard mesh. Other quantities of interest are defined and advanced in time on the staggered mesh. That is, on the standard mesh we evaluate and . Some derived fields such as and are also defined on the standard mesh. On the staggered mesh we evaluate and . Note that for plotting purposes we interpolate these latter fields onto the standard mesh (at the boundaries an extrapolation is performed).
We solve numerically equations (1)-(3) and (34) together with equation (5). When solving for the magnetic field we add the DC magnetic field contribution to . Space is discretized in and with a Fourier collocation scheme (Dahlburg & Picone, 1989) with isotropic truncation dealiasing. Spatial derivatives are calculated in Fourier space, and nonlinear product terms are advanced in configuration space. A second-order central difference technique on a uniform mesh is used for the discretization in (Dahlburg et al., 1986). Variables are advanced in time by a low-storage Runge-Kutta scheme. Several options are available: two-step second-order, three-step third-order, four-step third-order and five-step fourth-order (Carpenter & Kennedy, 1994). Results presented in this paper use the last option, as it permits the largest time step. Thermal conduction is advanced with second-order Super TimeStepping (Meyer et al., 2012).
HYPERION, which previously used only MPI for parallel execution, was modified for hybrid parallelization using a combination of OpenMP and MPI. The code retains its original MPI-only strategy of assigning groups of – planes to each MPI rank by decomposing the three-dimensional simulation domain along the direction. This keeps all of the data needed for FFTs in the periodic and directions local to each MPI rank. Scalability of the original MPI-only code was limited, however, because the maximum number of MPI ranks that could be used in a given simulation could not exceed the number of – planes in the domain. In the hybrid code, OpenMP multithreading is used to exploit parallel work within the groups of – planes assigned to each MPI rank, for example by computing one-dimensional FFTs in the – planes in parallel. This allows more CPU cores to be utilized than was possible with the MPI-only version and, for a fixed number of cores, reduces the overhead of MPI communication relative to the MPI-only code.
3.2. Simulation rescaling
The dimensionless numbers based on the physical parameters are given in Table 1. Note that physical Lundquist and Reynolds numbers are far too large for present day computations. Consider, for example, case A which has a characteristic flow velocity given by m s and a characteristic Alfvén velocity given by m s. We have a physical Reynolds number equal to:
and a physical magnetic Reynolds number equal to:
(here can be thought of as an Alfvén Mach number). Rather than use these numerically unresolvable Reynolds numbers we present the results obtained running the code with smaller Reynolds numbers that can be used with the currently achievable numerical resolution. For example in case A they are and , with a horizontal resolution of , i.e., for case A we use
These somewhat conservative values of the Reynolds numbers are taken based on previous numerical simulations of turbulent magnetofluids (Dahlburg & Picone, 1989). In order to keep the same relative efficiency of the radiative and conductive terms in the energy equation as in the real corona, we have rescaled and accordingly with the choice of , i.e., we set
so that for case A: and . This rescaling is motivated by the result found in the RMHD model (Rappazzo et al., 2008) that turbulent dissipative processes are independent of viscosity and resistivity when an inertial range is well resolved. The rescaled values are given in Table 2.
Numerical resolutions for all of the simulations are given in Table 2. These resolutions are smaller than our previous RMHD simulations. However, the present simulations integrate more complex governing equations, evolving eight different field components (number density, temperature, the magnetic vector potential field and the velocity field) compared to only two scalar fields in RMHD. In addition, the density stratification from the upper chromosphere to the corona constrains us, at present, to compute with a very small time step due to the large variation in the Alfvén speed along the loop.
We here discuss the transmission, storage and release of energy in the simulated coronal loop. At the system starts out in a ground state, defined by the constant initial axial magnetic field , zero magnetic field fluctuations , while the initial number density and temperature profiles are as described in section 2.2.3. The velocity field vanishes everywhere initially except at the top and bottom boundaries () as described in section 2.2.1. Radiation and thermal conduction are ramped up linearly until they attain their full values at .
4.1. Loop energization
The random velocity fields at the top and bottom boundaries twist the field lines in a disordered way (since the forcing velocity is not symmetric), creating a magnetic field component predominantly orthogonal to the DC magnetic field. Initially evolves quasi-statically thus growing linearly with time (Rappazzo et al., 2008; Rappazzo, 2015). But as soon as the intensity of grows beyond a certain threshold, that depends on the loop parameters, current sheets form on a fast ideal timescale, with their width thinning down to the dissipative scale in about an axial Alfvén transit time (Rappazzo & Parker, 2013). Furthermore thinning current sheets have been recently shown to be unstable to tearing modes with “ideal” (i.e., of the order of ) growth rates even for thicknesses larger than Sweet-Parker (Pucci & Velli, 2014). Overall this implies that once the field lines are twisted beyond a certain threshold, or equivalently once the magnetic field intensity grows beyond a corresponding threshold, the magnetic field is no longer in equilibrium and transitions on the ideal timescale to a magnetically dominated MHD turbulence regime, where magnetic fluctuations are stronger than velocity fluctuations (Einaudi et al., 1996; Dmitruk & Gómez, 1997; Rappazzo & Velli, 2011).
The work done by boundary motions on the magnetic field line footpoints corresponds to a Poynting flux whose axial component gives the energy flux entering the system from the -boundaries (e.g., see Rappazzo et al., 2008), where is the velocity at the -boundary and the magnetic field at the -boundary. Because the characteristic -boundary velocity timescales are much longer than the Alfvén transit time , initially grows linearly in time akin to . But once the dynamics transition to a fully turbulent regime the system attains a statistically steady state where the Poynting flux is on average balanced by energy dissipation, so that also velocity and magnetic field saturate fluctuating around their mean values.
Figure 2 shows the Joule heating and Poynting flux in dimensional form as functions of time for case A, integrated respectively over the entire volume and over both -boundary surfaces. Akin to our previous reduced MHD simulations the Poynting flux exhibits large fluctuations about its average value. This occurs because the Poynting flux contains the scalar product of the velocity at the -boundary , a given quantity, and the perpendicular component of the magnetic field that is determined by the nonlinear turbulent dynamics of the system. This input energy flux is therefore also a turbulent quantity with large fluctuations in time. Note that because the -boundary velocity field changes only slowly in time, the correlation between the velocity and the magnetic field at the -boundaries is always strong so that the Poynting flux is always positive, i.e., energy is never removed from the loop by the boundary motions (for a study of the correlation between boundary velocity and magnetic field see Rappazzo et al., 2010). For the latter to occur, the -boundary velocity field should change over time-scales comparable to or faster than the Alfvén transit time along the loop.
In addition, the random forcing of the kind we employ is not conducive to the formation of loop structures capable of storing a large amount of energy. Our forcing does not inject a net magnetic helicity - associated with inverse cascades and therefore potentially large energy storage - into our loop. With our forcing, the injected energy is significant enough to power a hot corona, but clearly not a major solar flare: most of the injected Poynting flux is efficiently converted into thermal energy as well as kinetic energy (the remaining injected energy persists as magnetic energy of the perturbed field).
Previous reduced MHD investigations have shown that the time-averaged Poynting flux varies approximately quadratically with the strength of the guide magnetic field () (e.g., Rappazzo et al., 2007, 2008). Figure 3 shows the Joule heating and Poynting flux in dimensional form as functions of time for case C, integrated respectively over the entire volume and over both -boundary surfaces. Recall that = 0.01 Tesla for Case A and =0.04 Tesla for Case C. From Figure 2 it is seen that the Poynting flux is of order J m s for Case A. From Figure 3 it is seen that the Poynting flux is of order J m s for Case C. Hence the Poynting flux increases by a factor of about twenty as the guide magnetic field increases by a factor of four, which within the error due to the short duration of the simulations is consistent with a quadratic relation.
Figure 2 also shows the Joule heating as a function of time for case A. It can be seen that the Joule heating is somewhat correlated with the Poynting flux – it exhibits the same pattern of relative maxima and minima but with a time lag (in Case A this time lag is about 200 seconds). This lag represents the time it takes the energy to propagate in to the loop and for the appropriate magnetic structure, i.e., electric current sheets, to form to permit dissipation. Similar remarks apply to Figure 3 that considers Case C.
4.2. Three-dimensionality and intermittency
The fluctuations seen in the Joule heating in Figure 2 are also evidence of temporal intermittency. Although the numerical simulations presented here have a relatively low spatial resolution they do present some level of intermittency. As expected, and as shown for this problem in our previous reduced MHD simulations, both temporal and spatial intermittency increase at higher resolutions, i.e., with Reynolds number (Rappazzo et al., 2008, 2010, 2013). Evidence of spatial intermittency for current density and temperature was already shown in our previous fully compressible simulations(Dahlburg et al., 2012) It was found that temperature is not uniform in space, rather it strongly increases in and around electric current sheets, forming similarly shaped spatial structures elongated in the direction of the strong guide field (Dahlburg et al., 2012). Note that both temporal and spatial intermittency should increase as the Lundquist numbers increase.
Our 3D compressible MHD simulations allow exploration of some of the thermodynamic implications of this turbulent and intermittent type of heating. The coronal loop in our simulation is in a self-consistent state, energetically determined by the balance between boundary-forcing, nonlinear dynamics, heating and cooling. To this must be added the non-trivial caveat that the energy flux entering the system is not determined simply by the -boundary velocity , but also by the nonlinear, turbulent dynamics developing in the loop, This is a consequence of the Poynting flux being given by the scalar product between the -boundary velocity and the orthogonal magnetic field component generated by the nonlinear dynamics . The heating is only due to resistive and viscous dissipation which happens at different locations at different times where small scales are produced, i.e., within current sheets continuously forming and disrupting. The behavior of the volume-averaged quantities, such as kinetic and fluctuating magnetic energies and resistive and viscous dissipation show a temporal behavior similar to previous RMHD results (e.g., Rappazzo et al., 2007, 2008, 2010).
Fully compressible simulations with HYPERION show the time evolution of the maximum electric current, as seen in Figure 4 for case D, which shows already some fluctuations. Figure 4 also shows the maximum temperature as a function of time, which correlates strongly with . Though not shown here, this correlation is seen to strengthen in simulations where the axial magnetic field strength is increased. Indeed, increasing the axial field brings our 3D MHD simulations closer in nature to the RMHD case.
Dahlburg et al. (2012) showed that the temperature is spatially structured, i.e., it is spatially intermittent. Figure 5 shows the and positions of in space for case D at selected times. It can be seen that wanders about, observationally resulting in a changing radiation emission pattern that can easily give the mistaken impression of an oscillating loop.
As seen in Figure 1, there is considerable variation in loop lengths and magnetic field strengths in the solar corona. Here we briefly consider the influence of the axial magnetic field strength on our results. Figure 6 shows how the maximum temperature depends on the axial magnetic field strength (cases A, B, and C). It can be seen that the maximum temperature increases with the magnetic field strength, with a slightly weaker than linear dependence on field strength.
How do the results change with horizontal ( and ) resolution and Lundquist numbers? Case D has twice the horizontal resolution as Case A. In addition, the Lundquist numbers are doubled and the Prandtl numbers are halved for Case D. Figure 7 shows how the maximum temperature depends on numerical resolution and the Lundquist numbers (cases A and D). Note that the RMS temperatures are not too different for the two cases, but the temperature oscillations in the higher Lundquist number case are somewhat stronger. Of course as in all turbulent systems, the full understanding of the high Reynolds number regime is non trivial, and it will be investigated in future work. Nevertheless our previous reduced MHD simulations indicate that dissipation rates and Poynting flux saturate at resolutions of about , and as shown here maximum temperature variation is weak at .
4.3. Effects of vertical () numerical resolution
The energization and response of the system depend on gradients at the boundaries (Bradshaw & Klimchuk, 2011; Mikíc et al., 2013). The most significant are the gradients of the magnetic vector potential, the temperature, and the number density. The Poynting flux depends on the magnitude of the and magnetic fields. In HYPERION these fields are computed as the curl of the magnetic vector potential. For the and magnetic fields there is a component due to the gradients of the and magnetic vector potential, hence the energization of the system depends on the accurate computation of these gradients. In the same way the response of the system to heating depends on the evolution of thermodynamic gradients near the boundaries. At first glance it might appear that we have under-resolved these gradients. The scale height of the initial temperature  can be estimated in nondimensional terms from equation 31. At the boundaries the temperature scale length is found to be 0.00789 (31.56 km in dimensional terms). For our system with (or 50000 km), this is resolved using 1585 uniformly spaced mesh points. Note that the initial number density scale height will be approximately the same. We will determine a posteriori how the resolution affects the energization and plasma response. Case A will be used as the baseline. In these simulations all of the physical parameters are the same as in case A; only the resolution is changed. Case A has 144 points in , case E has 288 points, case F has 576 points and case G has 1620 points (sufficient to resolve the temperature and number density scales at the boundaries). All of these cases have a dimensional magnetic field of 0.01 Tesla, so the stiffening effect of the DC magnetic field is weak relative to the other runs. Case G was only simulated for approximately 1800 seconds to allow for the computation of synthetic emissions and an emission measure, to be shown in a subsequent part of this paper.
A comparison of the Poynting flux for the simulations with different resolutions is shown in figure 8. It can be seen that all of the cases oscillate about approximately the same average value in time. In HYPERION the magnetic vector potential is advanced in time. Recall that the Poynting flux depends on the perpendicular component of the magnetic field. Hence the perpendicular magnetic field is a derived quantity – in particular it will depend on derivatives. The value of these derivatives will vary somewhat with the resolution. The nonlinearity of the system is reflected in the temporal variability shown in figure 8.
A comparison of the maximum temperatures for the simulations with different resolutions is shown in figure 9. The behavior here is similar to that exhibited by the Poynting flux. The maximum temperature oscillates about approximately the same value for all of the cases, with some variability seen in the details of the fluctuations. We conclude that, for the range of resolutions we have considered here, there is not a significant change in the numerical results.
4.4. Emission measure distribution
Since an emission line is formed over a relatively narrow temperature range, spectrally resolved observations can be used to infer the temperature structure of the solar atmosphere. This is often achieved by computing the differential emission measure distribution (DEM), which is a solution to the equation
Here and are the intensity and emissivity of the emission line. The emissivity includes all of the information specific to the atomic transition. The quantity is the DEM, which describes the conditions in the solar atmosphere, and is written
where is the electron density and is a coordinate along the line of sight. Further details and an application to solar observations can be found, for instance, in Warren et al. (2008).
The density and temperature for each voxel, that is each element of the simulation volume, were used to calculate the intensity of EUV spectral lines. We chose a set of 25 EUV lines ranging from K to K in temperature of formation. With the exception of Fe xviii 974.86 Å, the lines selected are all in the observed wavelength range of the EIS instrument on board Hinode (Culhane et al., 2007) and cover a variety of ionization stages of Mg, Si, Fe, S, Ar and Ca. Data from the EIS instrument have been routinely used to calculate emission measure distributions in different coronal conditions. The Fe xviii 974.86 Å was added to improve the constraints on the high temperature end and mimics the use of AIA 94 Å, which images Fe xviii (Warren et al., 2012). The emissivities for each line was calculated using the CHIANTI atomic database (Dere et al., 1997; Landi et al., 2013) assuming coronal abundances (Feldman et al., 1992) and the CHIANTI ionization equilibrium tables.
Figure 10 shows the density and temperature distributions along and across sections of the simulation domain, for time from 1770 s to 1830 s for the simulated cases A through D. Figure 11 shows at the top and bottom panels the synthesized intensities of a set of seven spectral lines integrated along the perpendicular direction to the loops’ axis, similar to observing the loops side-on. The panels in the center are the intensities of the loops’ mid-section integrated along five voxels in the direction. The integration times are 60 s in all four cases.
The emitting volume selected for the EM exercise corresponds to the apex of the loops. To compute line intensities we integrate the emissivities over a region 1750 km wide centered at the mid-plane of the computational domain. The volume of integration corresponds therefore to a area on a hypothetical plane of the image, 4000 km deep, namely a cross-section of about 5″, typical of loop observations in the corona. The integrated intensities in each spectral line, with an assumed uncertainty of 20%, serve as input to a Differential Emission Measure calculation algorithm. We used the Monte Carlo Markov Chain (MCMC) code (Kashyap & Drake, 1998), applied in the manner described in Warren et al. (2012). The MCMC algorithm calculates multiple (250) solutions with perturbed values for the intensities, providing an estimate of the error in the EM distribution calculation.
Figure 12 shows the EM solutions for the A, B and C runs, where the red line corresponds to the best-fit solution. The plot also shows, for context, color-coded lines representing the emission measure loci for the different atomic species of the EIS spectral lines. They illustrate the range of temperature dependent emission measures compatible with the intensity of a particular spectral line. A set of intensities at different temperatures constrain the EM distribution compatible with the complete dataset.
The computed emission measure distributions are qualitatively similar to the emission measure distributions computed from observed intensities, with a characteristic Gaussian-like distribution in the 1–4 MK temperature range. The weighted mean temperatures for the the three cases are respectively: 6.00, 6.13 and 6.22 MK. These are characteristic temperatures for the spectral windows (e.g., Fe xii) and filter bandpasses (171 Å, 195 Å) where we observe a significant fraction of the loop emission in the corona and are consistent with the peak emission measure temperatures of some of these loops (Warren et al., 2008). The bottom panel of Figure 12 shows as a comparison of the EM distribution for a sample area in active region NOAA 11082 (Figure 13), illustrating the similarities of the simulation emission measures with regions of the corona. The EM distributions are perhaps the easiest way to compare the simulations with the observations. The simulated intensities for individual spectral lines can differ from typical observed values by factors of 2–10. Such line-by-line comparisons are beyond the scope of this work, but will be considered in the future. Temperatures can be higher at the core of active regions, with emission measure distributions peaking at 4 MK and exhibiting asymmetric profiles with a steeper drop in the high temperature end (Warren et al., 2012).
Fig. 14 shows that the emission measure analysis is able to restore the true line-of-sight emission measure, that is the true density distributions as a function of temperature in the volume of integration (). This is applicable for the four cases A through D which are shown. The figure also demonstrates that we do not find significant differences in the EM distribution between the low and high resolution runs.
5. Discussion and conclusions
In this paper we have examined the dynamics of a coronal loop threaded by a strong axial magnetic field, where the field line footpoints are advected by random motions. Consistent with previous two-dimensional and 3D reduced MHD simulations (Einaudi et al., 1996; Dmitruk & Gómez, 1997; Dmitruk et al., 2003; Rappazzo et al., 2007, 2008; Rappazzo & Velli, 2011), and our previous fully compressible work (Dahlburg et al., 2012), the loop dynamics are found to be nonlinear, with a turbulent MHD cascade transporting the energy injected by boundary motions at large perpendicular scales to smaller scales. This leads to the formation of approximately field-aligned current sheets where energy is dissipated and around which temperature strongly increases, with small scales mainly in planes perpendicular to the axial magnetic field, along which both current and temperature structures are elongated. These small scales are not uniformly distributed in the loop, but rather the dynamics become increasingly more intermittent at higher Lundquist numbers both in space and in time. Localized electric current sheets continuously form and disrupt, leading to localized heating of the plasma on short time scales.
Our results show that the loop is the site of a continuous occurrence of reconnection events that present observations are unable to resolve both spatially and temporally. In the presence of a strong guide field magnetic reconnection occurs at the X-points of the orthogonal magnetic field component, leading to a continuous change of connectivity of the field lines that cross the reconnection sites (“interchange”) where the heating occurs (Schrijver, 2007). These many sub-resolution “heating” events add up to produce the observed emission, giving the impression at larger (observational) scales of a continuous diffuse heating. What is called “coronal heating” is actually the superposition of all events due to localized energy deposition along the subsequent different field lines that cross the reconnection sites, at the many current sheets elongated along the axial direction present in the loop volume at any give time (for a visualization of such current sheets, see, e.g., Rappazzo et al., 2008). Clearly the heating deposited along “strands” (small elemental flux tubes) is much smaller than the total dissipated power in the heating peaks shown in Figures 2–5, which is of the order of to Watts with a duration of about 1000 s. This suggests that the energy released along each strand is reasonably expected to be much smaller than J, which is about times the typical energy released in a flare. We expect the energy deposited along strands in typical heating events to exhibit a distribution with a peak at energy smaller than J, and plan to investigate more in depth the energy release mechanism and statistics in future work.
Recall that in our calculations we have used values of resistivity and viscosity that are much larger than the real ones. In the real Sun even smaller spatial and temporal scales are attained leading to even smaller energies being involved in each event. Evidence for this has been seen in RMHD calculations (Rappazzo et al., 2008; Rappazzo & Parker, 2013). Considering the values of resistivity and viscosity we adopt are much higher than the solar values, we expect in each event an energy release much smaller than that for a nanoflare. We will study this point in detail with higher resolution simulations in future.
We have employed an emission measure analysis to investigate whether the simulated intensities of the computational loops are representative of plasma in the corona and find great similarities both in peak temperature and distribution. We find that the simulated intensities and corresponding emission measures are in excellent agreement and they are accurate representations of the true emission measures. Testa et al. (2012), looking into 3D simulations of active regions, found that this method can be inaccurate when structures with significantly different density overlap along the same line-of-sight. The temperatures, which increase as the value of the axial magnetic field increases, are characteristic of warm loop structures visible in EUV channels. The loop is found to be a multi-temperature structure with isolated regions at temperatures of several million degrees and most of the loop at much lower temperatures. For each case presented in the previous sections the emission measure retains the same form for the entire hour of the computation in spite of the strong spatial and temporal intermittency.
In this paper we have adopted a Cartesian model of a coronal loop. The random motions at the boundaries shuffle the magnetic footpoints such that there is no ordered twisting of the field-lines. This random twisting does not facilitate the formation of magnetic structures that can store large amount of energy. Rather the system reaches a statistically steady state where integrated physical quantities (magnetic energy, Poynting flux, dissipation rates and radiative losses) fluctuate around their time-average values, so that the injected energy per unit time is entirely dissipated on the average (i.e., considering a long enough time interval). In a Cartesian model the only way to store a large amount of energy, that can subsequently give rise to larger magnetic energy release events, is to apply a spatially isolated and symmetric -boundary velocity. For instance, a vortex with intensity stronger than surrounding -boundary motions can twist the coronal magnetic field lines quasi-statically, thus storing magnetic energy, until a kink instability develops (Rappazzo et al., 2013). Similarly, an isolated -boundary vortex, even in the presence of strong nonlinear dynamics in the corona, can store energy at large spatial scales via an inverse cascade of energy, with subsequent energy release events in the micro-flare range (Rappazzo et al., 2013). Similarly also -boundary shear motions, isolated or stronger than surrounding motions, can store a large amount of energy as sheared magnetic field lines that can subsequently be released impulsively (Dahlburg et al., 2005, 2009; Rappazzo et al., 2010).
We want to stress the fact that the amount of energy entering from the footpoints is an outcome of the simulation since such energy depends on that cannot be specified as a boundary condition. That means that the loop nonlinear dynamics determines how much energy can be injected into the system and that the “heating function” cannot be assigned a priori. Furthermore the almost perfect correlation between and , confirming that the peaks in temperature are due to the local enhancement of the current, shows that the heating is due to local phenomena. These phenomena are the results of the complex perpendicular dynamics driven by -boundary motions which induce a local increase of the heating which in this framework is due exclusively to magnetic reconnection. Most of the dissipation occurs within localized current sheets which disrupt rapidly on Alfvén time-scales when their perpendicular size decreases to the smallest spatial scale present in our simulation. It is interesting to notice the good correlation between the behavior of the Poynting Flux and the energy dissipation. The two curves are very similar and shifted in time which means that when the system admits a bigger average energy flux from the two bases, current starts piling up locally leading to an increase of total energy dissipated and to a formation and disruption of localized current sheets. The average Poynting flux depends on the length of the loop. The time averaged flux is of the order of J m s for the hotter loop (case C) and J m s for the cooler one (case A ). The Poynting flux thus increases almost quadratically with magnetic guide field intensity as in previous reduced MHD studies (Rappazzo et al., 2008).
As already mentioned the resolution of our simulations is coarse compared with the real scales present in the corona and consequently we are using values of resistivity and viscosity which are much higher than the real ones, which are unachievable using present computers. We have verified that doubling the numerical resolution, and therefore halving the resistivity and viscosity, changes the significant results only weakly (as expected ). Previous reduced MHD simulations have shown that total dissipation rates and Poynting flux increase with increasing resolutions, saturating at resolutions of about . In a fully developed turbulent regime dissipation rates are not expected to depend on the Reynolds numbers beyond a certain threshold. Previous simulations suggest that the resolutions adopted in this paper are below but not too far from such a threshold (Rappazzo et al., 2008, 2010), confirming that the results presented here for the radiative losses are realistic. The challenging investigations of the dynamics in the high Reynolds regime and their impact on observations - expected to increase intermittency effects, and ultimately require kinetic calculations on the dissipation scale, is left to dedicated future works. Most phenomenological studies of coronal heating have so far concentrated on the thermodynamic response of the coronal plasma using one-dimensional hydrodynamic loop models with a prescribed heating function. These models have been in common use in coronal physics for over thirty years, and have been fundamental in providing a basic framework for a wide variety of coronal phenomena, including loop temperature distributions (Rosner et al., 1978), prominence formation (Oran et al., 1982) and, more germane to this paper, coronal heating (Reale, 2014; Klimchuk, 2006; Cargill & Klimchuk, 2004).
In the one-dimensional model heating is often represented as a constant (see, e.g., Chiuderi et al., 1981; Torricelli-Ciamponi et al., 1982). It also can be generalized to be a function of some combination of mass density, temperature and thermal pressure. In the latter cases the heating can be both spatially and temporally dependent. A particular case of interest is that of impulsive heating, in which the heating function is turned on for some span of time, and then shut off to allow the loop to evolve to a new equilibrium. The heating can be localized in the photosphere or appear as bursts in the corona.
The main limitation of one-dimensional models is that whatever functional dependence is chosen, the heating remains an ad hoc function and the main task of the researcher is to see which of these dependencies provides the best fit to observations. The chosen functional dependence is thought to lead to some understanding of which mechanism heats the corona, but the link with coronal heating theories and models remains essentially undetermined. For instance in the case of the Parker model investigated in this paper it is not obvious at all what heating function in 1D models would represent it better. One might be tempted, for example, to use a heating function varying very slowly with time, since photospheric motions have a very low frequency compared to the fast Alfvén crossing time. But previous reduced MHD and our simulations show that the system develops turbulent dynamics, with the timescale of the system strongly decreasing at smaller scales (Rappazzo et al., 2007, 2008; Rappazzo & Velli, 2011). This is independently confirmed by the recent finding that in such systems current sheets form on very fast ideal timescales, with their thickness thinning at least exponentially and reaching the dissipative scale in about an Alfvén crossing time (Rappazzo & Parker, 2013). Such thinning current sheets have also been shown to be unstable to tearing modes with ideal growth rates (Pucci & Velli, 2014), with the formation of many magnetic islands and X-points and the complex dynamics of so-called super-tearing or plasmoid instability (Bulanov et al., 1978; Biskamp, 1986; Loureiro et al., 2007; Lapenta, 2008) ensuing. Determining the equivalent heating function for 1D simulations from this framework of coronal heating is therefore a complex task. In particular such heating function for the Parker model has never been investigated, and therefore the 1D hydrodynamic models have not been de facto able to test it (Klimchuk, 2015).
As observational evidence has accumulated that many loops are not isothermal, it has become apparent that coronal loops cannot be modeled using a single flux tube (Schmelz et al., 2010). The narrow temperature distributions (Warren et al., 2008) and their transient nature (Ugarte-Urra et al., 2009) point to multiple structures and coherence. In an effort to account for these observations, refined multi-strand models (e.g Klimchuk, 2009) have been developed, in which an ensemble of one dimensional loops is assembled in an attempt to construct a three-dimensional loop. Our numerical simulations give strong support to a multi-temperature coronal loop structure whose specific temperature distribution is likely to depend on the loop parameters, similar to the Emission Measure Distribution shown in Figure 10, that we plan to further investigate in future work.
It is important to emphasize that, as far as the thermodynamics is concerned, we are solving the same equations that are used in a reduced form in one-dimensional models. Looking at the big differences in temperature appearing in the 2D plots in the mid-plane, it is easy to understand that the temperature profiles along different field lines originating in different points of the mid-plane can differ in a very substantial way since for all field lines the temperature at the footpoints is K. No field line can be considered representative of what happens in the loop. The limitation to one spatial dimension would leave out the self-consistent nonlinear dynamics with the most significant energy transfers, responsible for the formation of current sheets and thus the energy deposition at small scales, occurring in the perpendicular directions. Additionally for energy to be transferred from the magnetic field to the plasma magnetic reconnection must occur, hence magnetic field lines are constantly being broken and reconnected (Schrijver, 2007) strongly impacting the energy distribution along different strands.
- Aschwanden et al. (2007) Aschwanden, M. J., Nightingale, R. W., & Boerner, P. 2007, ApJ, 656, 577
- Bingert and Peter (2011) Bingert, S, & Peter H. 2011, A&A, 530, A112
- Biskamp (1986) Biskamp, D. 1986, PhFl, 29, 1520
- Bulanov et al. (1978) Bulanov, S. V., Syrovatskii, S. I., & Sakai, J. 1978, JETPL, 28, 177
- Bourdin et al. (2013) Bourdin, P.-A., Bingert, S., & Peter, H. 2013, A&A, 555, A123
- Bradshaw & Klimchuk (2011) Bradshaw, S. J., & Klimchuk, J. A. 2011, ApJS, 194, 26
- Carpenter & Kennedy (1994) Carpenter, M. H., & Kennedy, C. A. 1994, Technical Report 109112
- Cargill & Klimchuk (2004) Cargill, P. J., & Klimchuk, J. A. 2004, ApJ, 605, 911
- Cirtain et al. (2013) Cirtain, J. W., Golub, L., Winebarger, A. R., et al. 2013, Nature, 493, 501
- Chiuderi et al. (1981) Chiuderi, C., Einaudi, G., & Torricelli-Ciamponi, G. 1981, A&A, 97, 27
- Culhane et al. (2007) Culhane, J. L., Harra, L. K., James, A. M., et al. 2007, Sol. Phys., 243, 19
- Dahlburg et al. (2012) Dahlburg, R. B., Einaudi, G., Rappazzo, A. F., & Velli, M. 2012, A&A, 544, L20
- Dahlburg et al. (2005) Dahlburg, R. B., Klimchuk, J. A., & Antiochos, S. K. 2005, ApJ, 622, 1191
- Dahlburg et al. (2009) Dahlburg, R. B., Liu, J.-H., Klimchuk, J. A., & Nigro, G. 2009, ApJ, 704, 1059
- Dahlburg & Picone (1989) Dahlburg, R. B., & Picone, J. M. 1989, Physics of Fluids B, 1, 2153
- Dahlburg et al. (2010) Dahlburg, R. B., Rappazzo, A. F., & Velli, M. 2010, Twelfth International Solar Wind Conference, 1216, 40
- Dahlburg et al. (1986) Dahlburg, R. B., Zang, T. A., & Montgomery, D. 1986, Journal of Fluid Mechanics, 169, 71
- Del Zanna & Mason (2003) Del Zanna, G., & Mason, H. E. 2003, A&A, 406, 1089
- Dere et al. (1997) Dere, K. P., Landi, E., Mason, H. E., Monsignori Fossi, B. C., & Young, P. R. 1997, A&AS, 125, 149
- Dmitruk & Gómez (1997) Dmitruk, P., & Gómez, D. O. 1997, ApJ, 484, L83
- Dmitruk et al. (1998) Dmitruk, P., Gómez, D. O., & DeLuca, E. E. 1998, ApJ, 505, 974
- Dmitruk et al. (2003) Dmitruk, P., Gómez, D. O., & Matthaeus, W. H. 2003, Physics of Plasmas, 10, 3584
- Dorch & Nordlund (2001) Dorch, S. B. F., & Nordlund, Å. 2001, A&A, 365, 562
- Einaudi & Velli (1999) Einaudi, G., & Velli, M. 1999, Physics of Plasmas, 6, 4146
- Einaudi et al. (1996) Einaudi, G., Velli, M., Politano, H., & Pouquet, A. 1996, ApJ, 457, L113
- Feldman et al. (1992) Feldman, U., Mandelbaum, P., Seely, J. F., Doschek, G. A., & Gursky, H. 1992, ApJS, 81, 387
- Galsgaard & Nordlund (1996) Galsgaard, K., & Nordlund, Å. 1996, J. Geophys. Res., 101, 13445
- Georgoulis et al. (1998) Georgoulis, M. K., Velli, M., & Einaudi, G. 1998, ApJ, 497, 957
- Gold & Hoyle (1960) Gold, T., & Hoyle, F. 1960, MNRAS, 120, 89
- Golub et al. (2007) Golub, L., Deluca, E., Austin, G., et al. 2007, Sol. Phys., 243, 63
- Gudiksen & Nordlund (2002) Gudiksen, B. V., & Nordlund, Å. 2002, ApJ, 572, L113
- Hendrix & van Hoven (1996) Hendrix, D. L., & van Hoven, G. 1996, ApJ, 467, 887
- Hildner (1974) Hildner, E. 1974, Sol. Phys., 35, 123
- Kashyap & Drake (1998) Kashyap, V., & Drake, J. J. 1998, ApJ, 503, 450
- Klimchuk & Cargill (2001) Klimchuk, J. A., & Cargill, P. J. 2001, ApJ, 553, 440
- Klimchuk (2006) Klimchuk, J. A. 2006, Sol. Phys., 234, 41
- Klimchuk (2009) Klimchuk, J. A. 2009, The Second Hinode Science Meeting: Beyond Discovery-Toward Understanding, 415, 221
- Klimchuk (2015) Klimchuk, J. A. 2015, Phil. Trans. R. Soc. A, 373, 20140256
- Landi et al. (2012) Landi, E., Del Zanna, G., Young, P. R., Dere, K. P., & Mason, H. E. 2012, ApJ, 744, 99
- Landi et al. (2013) Landi, E., Young, P. R., Dere, K. P., Del Zanna, G., & Mason, H. E. 2013, ApJ, 763, 86
- Lapenta (2008) Lapenta, G. 2008, PRL, 100, 235001
- Lemen et al. (2012) Lemen, J. R., Title, A. M., Akin, D. J., et al. 2007, Sol. Phys., 275, 17
- Loureiro et al. (2007) Loureiro, N. F., Schekochihin, A. A., & Cowley, S. C. 2007, PhPl, 14, 100703
- Meyer et al. (2012) Meyer, C. D., Balsara, D. S., & Aslam, T. D. 2012, MNRAS, 422, 2102
- Mikíc et al. (2013) Mikíc, Z., Lionello, R., Mok, Y., Linker, J. A., & Winebarger, A. R. 2013, ApJ, 773, 94
- Oran et al. (1982) Oran, E. S., Mariska, J. T., & Boris, J. P. 1982, ApJ, 254, 349
- Parker (1972) Parker, E. N. 1972, ApJ, 174, 499
- Parker (1988) —. 1988, ApJ, 330, 474
- Parker (1994) —. 1994, Spontaneous current sheets in magnetic fields (New York: Oxford University Press)
- Patsourakos & Klimchuk (2006) Patsourakos, S., & Klimchuk, J. A. 2006, ApJ, 647, 1452
- Pucci & Velli (2014) Pucci, F., & Velli, M. 2014, ApJ, 780, L19
- Rappazzo (2015) Rappazzo, A. F. 2015, ApJ, 815, 8
- Rappazzo & Parker (2013) Rappazzo, A. F., & Parker, E. N. 2013, ApJ, 773, L2
- Rappazzo & Velli (2011) Rappazzo, A. F., & Velli, M. 2011, Phys. Rev. E, 83, 065401
- Rappazzo et al. (2010) Rappazzo, A. F., Velli, M., & Einaudi, G. 2010, ApJ, 722, 65
- Rappazzo et al. (2013) —. 2013, ApJ, 771, 76
- Rappazzo et al. (2007) Rappazzo, A. F., Velli, M., Einaudi, G., & Dahlburg, R. B. 2007, ApJ, 657, L47
- Rappazzo et al. (2008) —. 2008, ApJ, 677, 1348
- Reale (2014) Reale, F. 2014, Living Reviews in Solar Physics, 11, 4
- Rosner et al. (1978) Rosner, R., Tucker, W. H., & Vaiana, G. S. 1978, ApJ, 220, 643
- Schmelz (2002) Schmelz, J. T. 2002, ApJ, 578, L161
- Schmelz et al. (2010) Schmelz, J. T., Saar, S. H., Nasraoui, K., et al. 2010, ApJ, 723, 1180
- Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229
- Schnack et al. (1987) Schnack, D. D, Barnes, D. C., Mikíc, Z., Harned, D. S., & Caramana, E. J.1987, J. Comp. Phys., 70, 330.
- Schrijver (2007) Schrijver, C. J. 2007, ApJ, 662, L119
- Sturrock & Uchida (1981) Sturrock, P. A., & Uchida, Y. 1981, ApJ, 246, 331
- Testa et al. (2012) Testa, P., De Pontieu, B., Martínez-Sykora, J., Hansteen, V., & Carlsson, M. 2012, ApJ, 758, 54
- Torricelli-Ciamponi et al. (1982) Torricelli-Ciamponi, G., Einaudi, G., & Chiuderi, C. 1982, A&A, 105, L1
- Ugarte-Urra et al. (2009) Ugarte-Urra, I., Warren, H. P., & Brooks, D. H. 2009, ApJ, 695, 642
- Warren et al. (2008) Warren, H. P., Ugarte-Urra, I., Doschek, G. A., Brooks, D. H., & Williams, D. R. 2008, ApJ, 686, L131
- Warren et al. (2012) Warren, H. P., Winebarger, A. R., & Brooks, D. H. 2012, ApJ, 759, 141
- Winebarger et al. (2003) Winebarger, A. R., Warren, H. P., & Seaton, D. B. 2003, ApJ, 593, 1164
- Zacharias et al. (2011) Zacharias, P., Peter, H., & Bingert, S. 2011, A&A, 531, A97