On the occurrence of thermal non-equilibrium in coronal loops

Long-period EUV pulsations, recently discovered to be common in active regions, are understood to be the coronal manifestation of thermal non-equilibrium (TNE). The active regions previously studied with EIT/SOHO and AIA/SDO indicated that long-period intensity pulsations are localized in only one or two loop bundles. The basic idea of this study is to understand why. For this purpose, we tested the response of different loop systems, using different magnetic configurations, to different stratifications and strengths of the heating. We present an extensive parameter-space study using 1D hydrodynamic simulations (1,020 in total) and conclude that the occurrence of TNE requires specific combinations of parameters. Our study shows that the TNE cycles are confined to specific ranges in parameter space. This naturally explains why only some loops undergo constant periodic pulsations over several days: since the loop geometry and the heating properties generally vary from one loop to another in an active region, only the ones in which these parameters are compatible exhibits TNE cycles. Furthermore, these parameters (heating and geometry) are likely to vary significantly over the duration of a cycle, which potentially limits the possibilities of periodic behavior. This study also confirms that long-period intensity pulsations and coronal rain are two aspects of the same phenomenon: both phenomena can occur for similar heating conditions and can appear simultaneously in the simulations.


INTRODUCTION
Solving the coronal heating problem remains one of the biggest challenge in astrophysics. How can the tenuous plasma that constitutes the highest layer of the solar atmosphere be maintained at temperatures two orders of magnitude higher than that of the solar surface? One of the fundamental facets of this problem is to determine the spatial and temporal distribution of the heating.
Thermal non-equilibrium (TNE) is a phenomenon that can occur in the solar atmosphere when the heating is highlystratified (e.g., Mendoza-Briceño et al. 2005;Karpen & Antiochos 2008;Mok et al. 2008;Susino et al. 2010;Antolin et al. 2010;Lionello et al. 2013;Mok et al. 2016). This parclara.froment@astro.uio.no ticular localization of the heating produces chromospheric evaporative upflows that supply the coronal structure with dense and hot material. A thermal runaway is eventually triggered when the radiative losses overcome the limited heating at coronal heights. Condensations are formed locally in the corona and fall down to the loop footpoints along the magnetic field lines. Furthermore, if the heating is quasi-steady, i.e. with a high heating frequency compared to the typical cooling time, this phenomenon can be cyclic. Such a system has no existing thermal equilibrium and will undergo evaporation and condensation cycles with periods of a few hours. This highly nonlinear behavior is what we call TNE. The limit cycle solutions in coronal loops were first explored by Kuin & Martens (1982).
TNE has received an increasing interest in the recent years. The thermal runaway triggered by a local excess of density and leading to cool condensations is one of the standard explanations for the existence of cool materials in the corona. Such catastrophic cooling events can, for example, end up in the formation of prominences (Antiochos & Klimchuk 1991;Antiochos et al. 1999Antiochos et al. , 2000Karpen et al. 2006;Xia et al. 2011Xia et al. , 2014 or coronal rain (Schrijver 2001;De Groof et al. 2004Müller et al. 2003Müller et al. , 2004Müller et al. , 2005Antolin et al. 2010Vashalomidze et al. 2015). Coronal rain is widely observed in off-limb active regions (e.g., Antolin & Rouppe van der Voort 2012; Antolin et al. 2015), even if a proper quantification of the proportion of loops experiencing episodes of coronal rain is still lacking.
Nevertheless, the widespread existence of TNE in loops and, consequently, the widespread contribution of quasisteady footpoint heating, has been questioned (Klimchuk et al. 2010). However, recent modeling studies have shown that the role of TNE in the dynamics of loops may need to be revisited Mikić et al. 2013;Winebarger et al. 2014;Lionello et al. 2016;Mok et al. 2016). These modeling studies support the idea that TNE is probably common in coronal loops, with two main types of condensations involved. Mikić et al. (2013) in particular suggest that different regimes of TNE cycles could exist in loops. They differentiate cycles with complete condensations (CCs) where the temperature, locally in the corona, drops to chromospheric temperatures to form dense (up to ∼ 10 17 m −3 ) and cool blobs, related to the observed coronal rain, and cycles with incomplete condensations (ICs). For this other regime of TNE, the temperature stays at coronal temperatures, and the density remains relatively low (∼ 5 × 10 15 m −3 ). These two different regimes of evaporation and condensation cycles are obtained with different combinations of parameters of the loop geometry and of the heating strength and spatial distribution.
The early statistical study of long-period intensity pulsations by Auchère et al. (2014) brings new impetus to this debate. Using 13 years of data in the 195 Å channel of the Extreme-ultraviolet Imaging Telescope (EIT; Delaboudinière et al. 1995) on board the Solar and Heliospheric Observatory (SOHO; Domingo et al. 1995), the authors found that at least half of the active regions likely undergo these intensity pulsations with periods ranging from 2 to 16 hours. In particular, these pulsations are very common in coronal loops. They have also been observed with the coronal channels of the Atmospheric Imaging Assembly (AIA; Boerner et al. 2012;Lemen et al. 2012) on board the Solar Dynamics Observatory (SDO; Pesnell et al. 2012), during the six first years of the AIA archive (Froment 2016). Froment et al. (2015) studied three examples of such events in detail, with respectively periods of 3.8, 5.6 and 9.0 hours. The authors concluded that these events are related to TNE cycles. This study was focused on the thermal structure evolution of these loops, using simultaneously the six coronal passbands of AIA. The authors used both Differential Emission Measure (DEM), with diagnostics developed by Guennou et al. (2012aGuennou et al. ( ,b, 2013, and time lag analysis, using the same method as presented in Viall & Klimchuk (2012). Auchère et al. (2016a) have recently critically re-examined and confirmed the statistical significance of the detections used in Froment et al. (2015). Furthermore, the pulse-train nature of the observed signals highlighted by Auchère et al. (2016b) reinforced the conclusion that TNE is the cause of the long-period intensity pulsations observed. Finally, we recently conducted a modeling study in order to compare directly the simulations with these observations. In this 1D hydrodynamic simulation study (Froment et al. 2017, hereafter Paper I), we showed that with a highly-stratified and steady heating we can reproduce the main characteristics (average behavior integrated along the line-of-sight and longterm temporal variations) of the thermal behavior of the loop bundle of event 1 studied in Froment et al. (2015).
As already mentioned in Froment et al. (2015), in the active regions where long-period intensity pulsations are detected, only one loop bundle (or in some rare cases two) shows this kind of behavior. The automatic detection algorithm used may have missed events due to the rather strict detection thresholds that we used. It is also likely that some events are missed by the Fourier detection if they are not strictly periodic. Furthermore the background and foreground emission could mask the time-varying signal in certain cases. However, some properties, such as the geometry of the magnetic field lines and the heating characteristics, could favor the TNE cycles for only some loop bundles. In the present paper we explore the sensitivity of TNE occurrence as a function of the heating strength and stratification, for several loop geometries.
The heating parameters used in the simulation presented in Paper I have been chosen among hundreds of heating combinations tested for a single loop geometry, which corresponds to a pulsating loop bundle observed with AIA. We present here an extensive parameter-space study which justifies this particular choice. To extend our analysis, we also pick another loop geometry, corresponding to a non-pulsating loop from a neighboring region, and use a semi-circular one as a test sample to do the same analysis and thus test the influence of the loop geometry.
This parameter-space exploration is presented in Section 2. In Section 3 we discuss our results regarding the occurrence of TNE cycles within the parameter space explored. Then, in Section 4 we examine the properties of the different types of loop systems produced in our simulations (with TNE or not) in regard of the observed loops properties. Finally, we summarize our results in Section 5.

PARAMETER-SPACE SCAN
For this parameter-space study, we use the same 1D hydrodynamic code as in Mikić et al. (2013) and in Paper I. The 1D description is particularly suited for this kind of study as multiple configurations of loops can be easily tested. The loop geometries used in these simulations, except for one loop (loop A, see Section 2.2), are from the Linear Force-Free (LFFF) extrapolations presented in Section 2.1 of Paper I. These extrapolations are made using magnetograms from the Helioseismic and Magnetic Imager (HMI; Scherrer et al. 2012) corresponding to the active regions of event 1 of Froment et al. (2015), i.e. NOAA AR 11499, and NOAA AR 11501, a smaller adjacent active region. Some selected extrapolated field lines are presented in Figure 1. We detected intensity pulsations (with a period of 9.0 hr) in a large loop bundle of NOAA AR 11499, with a very clear signal; the probability that this detection was caused by noise is 1.7 × 10 −8 (Auchère et al. 2016a). This area is delineated by the orange contour in Figure 1.
In the field-of-view presented in Figure 1, outside of the orange-contoured bundle, no other loop bundle shows a longperiod pulsating behavior. The evolution of the pulsating loops has been followed from June 03, 2012 18:00 UT to June 10, 2012 04:29 UT using AIA data.

Method and parameters explored
We choose to focus on three different loop geometries and to scan various heating configurations for these loops. In addition of the loop geometry that matches the pulsating loop bundle already used in Paper I (noted here as loop B), we use a semi-circular geometry as a control sample (loop A) and we picked another loop geometry from the LFFF extrapolation that matches with a non-pulsating loop bundle as observed with AIA (loop C). The field lines corresponding to loop B and loop C, and matching visible loops bundles in the 171 Å AIA image, are indicated in Figure 1. Loop A is an ad-hoc loop and is therefore not present in the field-of-view.
The 1D hydrodynamic simulations are made using the same initial conditions and assumptions as in Paper I (see description in Section 2.2). The only difference here is that instead of using the Spitzer thermal conductivity as in Paper I, we use an option present in the code  to artificially broaden the transition region at low temperatures. This modification of κ allows to reduce the steep gradients below a cut-off temperature (chosen as T c = 250, 000 K here) with a minimal effect on the coronal solutions, as described in Lionello et al. (2009) andMikić et al. (2013). In that way we can afford to use bigger mesh cells and thus fewer mesh points, i.e. 10,000 points, than in Paper I. This technique is particularity suited for this study, since we wish to scan a large area of the space of parameters. Some of the runs presented in this paper were repeated with the classic, unmodi-fied, Spitzer conductivity using more mesh points (typically 100,000 points). It is the case for the particular simulations that we choose to present in detail in Section 4 (see Figure 16) and some examples presented in Froment (2016). The overall pattern is not affected by this technique but some differences related to the precise nature of condensations can appear between the simulations using the Spitzer thermal conductivity and the ones using modified conductivity.
We choose a simple heating function that can be tuned with three free parameters. This heating function is the same as the one used in Paper I (see Equation 2): where g(s) = max(s − ∆, 0) and ∆ = 5 Mm is the thickness of the chromosphere, where the heating is constant.
H(s) is the volumetric heating rate, expressed in W m −3 . H 0 is the value to which H(s) tends at the apex and (H 0 + H 1 ) is the value of the heating in the chromosphere. λ 1 and λ 2 are the scale lengths for the energy deposition at the eastern and western leg of the loop, respectively.
For each loop geometry we test several values of H 1 , λ 1 and λ 2 . We choose to fix the H 0 value for each loop geometry, to a value that allows us to obtain a static loop (we set H 1 = 0 to use a uniform heating) with an apex temperature around 1 MK.
For each loop geometry, we explored values of H 1 in increments of a factor of 2. In that way we can easily compare the simulations. Note that the value of the factor itself is arbitrary. For each value of H 1 explored, we test a large set of combinations of λ 1 and λ 2 , specified in percentage of the total length L of each loop. The scan cube is then (H 1 , λ 1 , λ 2 ).
For each simulation, we can define the heat flux, i.e. the total heat that the loop receives over its length, normalized to the first loop footpoint cross-section area : with A(s) the cross-sectional area of the loop which is ∼ 1/B(s). As we will see in Section 2.2 (Figures 2 and 3), the magnetic field strength is similar at both end of each loop geometry chosen. The heat flux normalised at the second footpoint, Q 1 , is thus always similar to Q 0 . We thus only consider this latest value. The Q 0 heat flux value does not give information about the asymmetry of the heating. However, looking at Equations 1 and 2, we can see that the heating in a leg will be dominant when the scale height of energy deposition is larger, and the more the loop area expansion is important in this leg.
It is worth noting that the way we choose to explore the parameter space of the heating strength and stratification implies that Q 0 changes between two slices through the scan The orange contour delimits the area of the pulsations (9.0 hr of period) detected in the 335 Å passband of AIA (see Figure 4 in Froment et al. 2015), for a sequence of data between June 05, 2012 11:14 UT and June 08, 2012 11:16 UT. The red field lines match this contour. The white arrows indicate the loops B and C chosen for the parameter-space study. See Section 2.1 of Paper I for further details regarding these extrapolations. cube, i.e. between simulations with the same stratification (λ 1 , λ 2 ) but with different values of the heating at the footpoints (H 1 ). Indeed, Q 0 changes from one simulation to another. We could have chosen a different parametrization, fixing Q 0 instead of H 1 for each exploration of λ 1 and λ 2 . However, we confirmed that this different parametrization did not change our conclusions. In addition, whatever the parametrization, since Q 0 and H 1 are linked, the results can be explored at Q 0 constant. We tested as well an alternate way of creating asymmetry in the heating along the loop by applying a different amount of heating at each footpoint. We verified that this did not change the conclusions of this paper either.

Loop geometries
In Figure 2, we present the loop geometry of loop A, the test case semi-circular geometry we use for the first set of simulations. In Figure 3, we present the loop geometry of loops B and C which are from two field lines extracted from the LFFF extrapolations of the active regions presented in Figure 1. These loop geometries, even if from single fieldlines, are selected to model the behavior of the corresponding loop bundles. When we compare the simulations with observations, a simulated loop represents the average behav- ior of a loop bundle, whose individual threads all have similar properties.
In these figures, we show the quantities that are input to the simulations, i.e., the loop profile (the altitude of each point of the loop), the gravity projected along the loop (see Equation 1 of Paper I), the magnetic field strength along the loop and the loop area A(s), normalized to the first footpoint. For both loops, s = 0 corresponds to the eastern footpoint.
These three loop geometries have the following characteristics: Loop A -This loop is a semi-circular loop with the same length as loop B, i.e L = 367 Mm. The magnetic field along the loop is given by: Equation 4 of Mikić et al. 2013). The loop expansion factor reaches a value of 11 at the loop apex (i.e. at 58 Mm of height).
In input of the simulations we select the mesh spacing to be ∆s = 19 km in the chromosphere and transition region, increasing to ∆s = 190 km in the corona.
Loop B -This is the same loop that we studied in Paper I. It corresponds to the pulsating loop bundle detected in AIA data. It is quite a large and asymmetric loop with L = 367 Mm. As in Paper I, s = 0 corresponds to the east-  at s = 162 Mm, i.e. to the east of the loop apex. As seen in Paper I, the large low-lying portion at the eastern footpoint is due to the magnetic topology in this area, i.e. a low-lying null-point and many bald patches. As for loop A, the mesh spacing is ∆s = 19 km in the chromosphere and transition region, and ∆s = 190 km in the corona.
Loop C -This field-line corresponds to a non-pulsating loop bundle in the AIA data. It is located in the small active region at the West of NOAA AR 11499. We did not find any intensity pulsations in any of the loops of this region. The loop chosen is shorter than the previous ones, with L = 139 Mm. As for loop B, s = 0 corresponds to the eastern footpoint. The apex of this loop has an altitude of 42 Mm at s = 77 Mm, i.e. s = 0.55 L. The loop expansion is quite large, with a maximum value of 82 at s = 92 Mm, i.e. to the west of the loop apex. The magnetic field strength is about the same at each footpoint with respectively 1100 G and 1200 G. Using the same number of mesh points as for the two other loops, the mesh spacing is smaller for this loop: ∆s = 7 km in the chromosphere and transition region, increasing to ∆s = 70 km in the corona.

Exploration of heating parameter space
For each loop geometry, we present hundreds of simulations, which is still only a fraction of the parameter space we explored. We choose here to focus only on the area of the parameter space surrounding the simulations showing TNE cycles.
The results of the exploration for each loop geometry are displayed using grid plots, each individual plot showing the temperature evolution along the loop for three days of simulation time. Each figure represents a scan of the λ 1 and λ 2 values for a given value of H 1 , i.e. a cut through the scan cube. By this means we can analyze the global behavior of the simulations within the parameter space.
We choose arbitrarily to show only the temperature, although the density or the velocity would be also suitable to present a different view of the loop behaviors. Nevertheless, the density and the velocity profiles are shown for a few examples analyzed in detail (see Figure 16). Temperature, density, and velocity averaged around the apex are also shown in Section 2.3. Throughout the present paper the apex area is defined, as in Paper I, as the part of the loop above 90% of the apex height. Temperature, density, and velocity are averaged in this area, to avoid quoting values at a single location.

Physical limitations on the domains explored
As mentioned earlier, we explored a very large range of heating for each loop geometry, in particular in terms of scale heights: λ 1 and λ 2 . We did not limit the study to specific ranges that would be appropriate to the magnetic field configuration of each loop system. As a consequence, it is likely that the domain of the parameter space presented in this paper is larger than a "realistic" one that would be constrained by the magnetic field strength, for example. This was done deliberately, to separate the heating mechanism from the magnetic field, so as not to make any strong assumptions about the heating. Thus, we are maybe scanning scale heights that are not possible in reality for a given loop geometry.
The choice of the scanned H 1 values is justified by the temperature and density of the simulated loops, i.e. we aimed to produce loops with coronal temperature (close to 1 MK and at most 4 MK for the temperature and 10 14 m −3 to 10 15 m −3 for the density). We explored as well H 1 producing cooler loops but for these cuts through the scan cube, not enough density was injected, so no TNE was produced. We tested also higher H 1 but then the loops are extremely hot which does not correspond to our observations. Note that given the large loop expansion of loops B and C (respectively 38 and 82), some of the values of H 1 give really high Q 0 (up to 20 × 10 4 W m −2 for loop C, see Table 1). However, it can be delicate to interpret the absolute values of Q 0 . In our simulations the chromosphere and low-transition region are not well modeled. The part of the heating that is radiated away in these layers may be unrealistic and further studies would be needed to quantify it. Therefore, we are not interpreting the absolute values of the strength of the heating but rather the evolution of the size of the TNE domain when H 1 is doubled (relative heating strength).

Criteria to distinguish between the different behaviors
Within the parameter-space, we detect the TNE cases and determine the nature of the condensations, using only the temperature profiles.
TNE events -They are detected within the parameter-space using Fourier analysis. We look at periods between 2 and 16 hours as we did for the AIA observations. For each simulation, we look at the evolution of the temperature averaged around the loop apex. We do not consider the beginning of these temperatures curves, i.e. the first 10 hours of the simulations, to minimize effects related to the initialization of the simulation, and concentrate on the asymptotic behavior only. Simulations are labeled as TNE events when the Fourier power, for at least one frequency bin, is 20σ above an estimate of the average local power. Moreover we discard the simulations with amplitudes 1 smaller than 0.2 MK in the second half of the simulations to avoid having to much loops with damped cycles in the TNE domain.
Distinction between ICs and CCs -We look at the nature of the condensations not only around the apex but all along the coronal part of the loops. In fact if some CCs occur low 1 Temperature difference between the maximum and the minimum. enough in one of the loop legs, the evolution of temperature around the loop apex is not very different from an IC case (see e.g. the first, noted as IC and the third, noted as CC 2, simulations presented in Figure 16). Testing if the condensations are complete or incomplete only around the apex is then insufficient. The coronal part of the loop is defined as the parts of the loops above 10 Mm of altitude. The CC cases are then detected if the temperature drops locally under 0.5 MK. The other TNE cases are labeled as IC.

Loop A
For this group of simulations H 0 = 1 × 10 −7 W m −3 . We scan three values of H 1 : 640H 0 , 1280H 0 , and 2560H 0 . For each value of H 1 , λ 1 and λ 2 are scanned between 2% and 11% of L, i.e. we test 10 values between 7.3 Mm and 44.0 Mm. Every combination is tested, so we have eventually 3×10×10 different heating configurations, i.e. 300 simulations. All these simulations are presented in Figure 4.
Some of the simulations show cyclic cases of evaporation and condensation. Around a restricted domain of TNE cycles, the simulations produce either continuous siphon flows or loops reaching a static equilibrium. For each simulation, we indicate: • if it is a TNE case with a white dot, and either CC (for Complete Condensation) or IC (for Incomplete Condensation). CC cases can be visually identified by the dark blue and purple drops in the temperature evolution (temperature ≤ 0.5 MK); • SF is stated for continuous Siphon Flows; • and SE for Static Equilibrium.
In order to determine whether a simulation exhibits TNE cycles and what the nature of the condensations are, we use the criteria presented in Section 2.3.2. We see that TNE cycles are encountered only around the diagonal of each of these square grid plots, i.e. for simulations for which a symmetric heating function is applied. The upper limit for theses cycles is λ 1 = λ 2 = 33.0 Mm (i.e. 9% of L), i.e., the solutions with λ 1 or λ 2 larger than this value are stable. We also notice that the more heating is applied (H 1 high, and consequently Q 0 high), the more the TNE domain extends away from the diagonal.
In Figure 10, we show the maximum temperature and density averaged near the loop apex and the averaged velocity over the loop apex. This temperature is most of the time coronal (with very few simulations below ∼ 0.6 MK). It increases to 4 MK for large values of H 1 , λ 1 and λ 2 . We see a clear signature of heating for symmetric heating cases (with λ 1 ∼ λ 2 ). Indeed, comparing with Figure 4, we can identify that within the TNE domain (indicated by the field of white dots in Figure 4)  . Temperature evolution for 300 loops from 1D hydrodynamic simulations using the loop A geometry (semi-circular, see Figure 2). Each grid plot shows a cut through the heating parameter scan cube, i.e. each grid plot corresponds to a different value of H 1 (the heating imposed at the footpoints  Figure 5. Same as Figure 4 for the temperature evolution for 144 loops, using the loop B geometry (see Figure 3). Parameter scan with H 1 = 6.4 × 10 −5 W m −3 . λ 1 is scanned between 7% and 18% of L, i.e. 25.7 Mm and 66.  Figure 6. Same as Figure 5 with H 1 = 12.8 × 10 −5 W m −3 . The color scale is saturated at 4 MK. The red dots indicate TNE cases studied in detail in Section 4.1 and shown at higher resolution in Figure 16.  Figure 7. Same as Figure 4 for the temperature evolution for 144 loops, using the loop C geometry (see Figure 3). Parameter scan with H 1 = 4.0 × 10 −5 W m −3 . λ 1 and λ 2 are scanned between 4% and 15% of L, i.e. between 5.6 Mm and 20.9 Mm. The black dashed line indicate the symmetric heating simulations (λ 1 = λ 2 ), which corresponds to the diagonal of the grid plots for this loop geometry. The color scale is saturated at 2 MK.    easily. The maximum density plot shows also clearly a condensation pattern for the TNE simulations. We notice that this maximum density is up to ∼ 10 15 m −3 which is a reasonable value for a large coronal loop (Reale 2014). However these values are quite low for CCs cases but we have to bear in mind that we average these quantities around the apex before determining the maximum, and as we will see for other loop geometries, the density peak is not necessarily reached around the apex. We notice that the velocity around the loop apex is quite high (∼ 100 km s −1 ) when λ 1 +λ 2 < 40 Mm, when the heating is very stratified. The velocities are lower close to TNE cases.
We use the velocity maps 2 to analyze the loops without cycles. For each H 1 , in the region where λ 1 > λ 2 , we encounter loops whose evolution is dominated by siphon flows to the right footpoint (i.e. the less heated footpoint, see SF + in Figure 4). We witness the reversed siphon flows when λ 1 < λ 2 (see SF -in Figure 4). We have thus continuous siphon flows to the less heated footpoint when the heating is strongly asymmetric. The last main behavior encountered is static equilibrium (velocity close to zero along the loop, see SE in Figure 4), for the loops along the diagonal and with λ 1 , λ 2 > 33.0 Mm. Note that the TNE cases show periodic siphon flows (see Figure 16) but that the cases pointed out as SFs here are simulations dominated by continuous siphon flows for several days. For each values of H 1 , the TNE, SF +, SF -and SE domains do not overlap.
To conclude for this semi-circular loop geometry, we found that the majority of TNE cases produced CC cycles: between 82% and 89% of the TNE cases are CCs, depending on the H 1 used. This is consistent with the results of Mikić et al. (2013, see e.g. Case 7). A few ICs are encountered at the boundaries of the TNE domain (i.e. when the heating is asymmetric) when the total heating is increasing (i.e. for higher H 1 ). The domain within the parameter space in which loops are undergoing TNE cycles is rather restricted to symmetric (or close to) heating cases. We notice that there is more dispersion around the diagonal when the total heating is increased.
The TNE cycles are also located in a restricted domain but are now shifted to the region where λ 1 > λ 2 , i.e. asymmetric heating cases when the eastern footpoint (leg) is heated more than the western one. Compared to loop A, the sub-H 1 = 6.4 × 10 −5 W m −3 ; H 1 = 12.8 × 10 −5 W m −3 H 1 = 6.4 × 10 −5 W m −3 ; H 1 = 12.8 × 10 −5 W m −3 H 1 = 6.4 × 10 −5 W m −3 ; H 1 = 12.8 × 10 −5 W m −3 Figure 11. Same as Figure 10 for the simulations using the loop B geometry. The velocity is displayed between ±15 km s −1 . field of the parameter space presented here is thus not centered around the symmetric heating cases (indicated by the black dashed line). The scale heights for the energy deposition have to be larger than for loop A to reach TNE conditions.
Only one simulation shows TNE cycles with a symmetric heating function. For this simulation, λ 1 = λ 2 = 33.0 Mm (for H 1 = 1280H 0 , see Figure 6). However, this simulation is at the edge of the TNE domain. The envelope of the TNE domain is roughly restricted to the following values: 1.0 < λ 1 /λ 2 < 3.4, though the exact shape of the TNE domain is Figure 12. Same as Figure 10 for the simulations using the loop C geometry.The velocity is displayed between ±15 km s −1 . more complex. We notice that the range of λ 1 taken by TNE cases (between 7% and 17% of L) is wider than the one of λ 2 (between 2% and 9% of L). This is probably due to the asymmetry of the loop geometry, the field line from the LFFF extrapolations being skewed toward one footpoint.
As for Loop A, the more the loop is heated, the wider the TNE domain is. Looking now at the condensations in these TNE cases for the two H 1 values scanned, we notice that respectively 56% and 45% of them have cycles with incomplete condensations. Moreover, these IC cases tend to be at the edges of the TNE domain.
In the same way as for loop A, the maximum of the averaged apex temperature and density and the mean apex velocity are displayed in Figure 11. The values reached for both temperature and density are similar to the ones for loop A. We also see a larger apex temperature in the TNE domain, compared to the surrounding SF cases, as was the case for loop A. The velocities at the apex are much smaller than for loop A (maximum 15 km s −1 ). But we observe the same pattern of velocity evolution within the parameter space, i.e. higher velocities when the heating is highly stratified, λ 1 + λ 2 < 60 Mm in that case. Moreover, as for loop A, velocities at the apex are lower in the TNE domain.
At each side of the TNE domain, the simulations are dominated by siphon flows, the direction depending on the asymmetry of the heating. We also find a few SE cases (see location in Figures 5 and 6). Note that the white pixels in Figure 11 does not necessarily indicate SE cases as TNE cases can show zero velocities at the apex. SE cases point out no flows along the loop for most of the three days of the simulation.
Finally, it is worth noting the presence of high-frequency fluctuations (wavy pattern) in these simulations, especially around the eastern footpoint. We surmise that it is probably due to a combination of the thick chromosphere at this footpoint and the numerical treatment of the transition region. Loop B has a portion that is almost tangent to the photosphere at its eastern leg (i.e., with small projected gravity, see Figure 3 and Paper I). Indeed, this sawtooth pattern is not observed for the other loop geometries or for the highresolution simulations in Figure 16.

Loop C
For this last loop geometry, we parameterize the heating function with H 0 = 2 × 10 −6 W m −3 . We scan three values of H 1 : 20H 0 , as presented in Figure 7, 40H 0 , as presented in Figure 8, and 80H 0 , as presented in Figure 9. For each value of H 1 , λ 1 and λ 2 are scanned between 4% and 15% of L, i.e. we test 12 values between 5.6 Mm and 20.9 Mm. Every combination is tested, so we have in total 3 × 12 × 12 = 432 simulations.
With this loop geometry, we notice that TNE cycles appear first for symmetric or slightly asymmetric heating conditions (λ 1 > λ 2 ) when H 1 = 20H 0 . Then when we increase H 1 , more TNE appear for asymmetric heating conditions, specially for λ 1 > λ 2 . Finally, when H 1 = 80H 0 , we can notice that the TNE domain becomes much larger than for loop A and B. In particular, there does not appear to be a limit to TNE for large scale heights. However, the TNE domain remains limited as barely any TNE cases appear for a very high stratification of the heating (λ 1 or λ 2 smaller than 8.3 Mm).
We notice also that most of the TNE cases have CC cycles (respectively 0% of the TNE events for the two first value of H 1 and 1% for the last one). The siphon flow cases surround the TNE domain, with IC cycles starting to appear on the boundary of this domain. Figure 12 shows the maximum of the averaged apex temperature and density and the mean apex velocity. The temperature and density reached are similar to those of the other loops but we clearly see the large range covered by the CCs events (high temperature and high density). The velocities are close to the ones observed for loop B but the pattern we observed for loops A and B, i.e. higher velocities for small heating scale heights is not as clear here.

Conditions that favor TNE and constraints on the heating
Scanning the parameter space of heating configurations for different loop geometries, we have noticed that the distribution of the occurrence of TNE depends on the loop geometry. However, from this study, it seems that we are able to produce TNE-favorable conditions for any loop geometry. TNE will occur if the heating strength is sufficient to produce a loop dense enough to create a thermal runaway at high altitudes, and if this heating is deposited on specific scale heights.
From the heating parameter-space scan that we conducted with three different loop geometries, we can conclude that a stratified heating is a necessary condition to produce TNE, but it is not sufficient, as already found by many authors (e.g., Müller 2004;Susino et al. 2010;Mikić et al. 2013). For each loop geometry, the system undergoes TNE cycles for specific heating stratifications: • λ 1 λ 2 for loop A; • λ 1 > λ 2 for loop B; • For loop C, we observe two behaviors: λ 1 ∼ λ 2 when H 1 is small, λ 1 or λ 2 > 8.3 Mm when H 1 is large.
As stated before, in this paper we present only a sub-field of the parameter-space that we scanned. For each loop geometry, we also tested smaller and larger values of H 1 than the values presented here. However, when H 1 is too small, we do not reach typical coronal loops temperatures and densities and a large majority of the loops do not show any TNE cycles. When the H 1 are too large, the temperature of the loops is too high (> 4 MK) compared to the warm pulsating loops observed with AIA.
For all the loop geometries, the more heating (H 1 ) is applied at the footpoints the less the system requires heating symmetry to achieve TNE cycles. In other words, higher H 1 leads to more TNE within the parameter space. Higher H 1 induces more chromospheric evaporation, which results in denser plasma at coronal altitudes, favoring the thermal runaway. Indeed, the input heating H 1 has to be sufficient to inject the density required for triggering TNE events. The loops that we are modeling in this paper are quite large (367 Mm and 139 Mm) and thus a large chromospheric evaporation is needed to injected enough density, which can explain the large Q 0 values required to have TNE cycles (see values in Table 1).
It is worth remembering that the loop C geometry is extracted from a field-line corresponding to a loop bundle that is not undergoing cycles in the observations. Interestingly, it is for this loop geometry that the TNE cycles have the highest probability 4 to occur, according to our simulations results with the highest H 1 tested. If our model (quasi-steady stratified heating) is indeed correct, it would mean that we can constraint the heating for these particular loops. That would mean that loop C, which is not showing any pulsations in 4 With respect to the explored volume of the parameter space. the AIA data, is not heated enough at the footpoints to inject the excess of density needed in the loop bundle to trigger the thermal runaway and/or that the heating is very stratified. In the same way, loop B is showing pulsations in the AIA data so we can guess that for this loop bundle the heating is asymmetric, stratified and relatively important.

Exploration de-correlated from the magnetic field strength
For each of the geometries tested, not all the stratified heating configurations lead to TNE. The area where TNE occurs is limited to some range in the heating parameter space. This leads to the question as to whether the area explored within the parameter space is realistic. In particular, the heating is somewhat correlated to the magnetic field strength (see e.g. turbulent models, in e.g. Rappazzo et al. 2007) and therefore we may have explored heating parameters that are unrealistic. On the other hand, the strength of the magnetic field along the loops does not take into account the magnetic topology, which necessarily influences the heating as well (formation of separatrices, preferential reconnection sites, e.g. Aly & Amari 1997;Pariat et al. 2009;Parnell et al. 2010;Wyper et al. 2012). The heating parameter-space scan for loop B shows that we can produce TNE cycles for this loop geometry only with asymmetric heating profiles. This heating configuration was validated a posteriori in Paper I by the magnetic topology found around the eastern footpoints of this loop bundle, that can favor continuous reconnection and thus enhanced heating.

Common characteristics of TNE events
From the analysis of the flows, using in particular the averaged apex velocity (see Section 2.3), we noticed that the siphon flows are more intense for short heating scale heights. Moreover, they tend to be weaker close to the TNE conditions (see Figures 13 to 15 ).
We examine also some characteristics of the cycles of the TNE cases. Figures 13, 14, and 15 show, for each simulation, the periods of the cycles and time lags between the temperature and density averaged around the apex.  Figure 13 for the simulations using the loop C geometry.
Time lag between T e and n e -We compute the time lag between the temperature and density evolution. This delay is also a characteristic of TNE cycles; it is a signature of TNE when combined with the periodicity. It also explains the systematic cooling pattern observed between EUV channels, the intensity peaking first in the hotter channels and then in the cooler ones (e.g., Viall & Klimchuk 2012). In case of TNE events this observed cooling can be explained by a faster rise of the temperature than the temperature fall combined (or only if) the density is low during the heating phase compared to the cooling phase (see Section 3.2.3 in Paper I for more details).
This time lag is given here by the peak of the crosscorrelation between the average temperature and density curves around the loop apex. We choose to display them as a fraction of the period in order to compare the cases more easily. We explore systematically positive time delays between 0-60% of the TNE cycles period. Indeed, in our knowledge no TNE simulations have been reported to show an increasing of the density before the temperature and the signal being periodic we avoid in that way to detected spurious negative time lags between T e and n e . For loop A, we notice that there are very long delays when there are very strong CCs. In some cases this delay is close to the period (see the notes in Table 1). Moreover, for some CC cases the shape of the temperature and density curves are very different which leads to poor cross-correlations values and underestimated time lags (not catching the strongest density peak). For loops B and C the time lags tend to be maximum close to symmetric heating cases, otherwise becoming quite uniform within the TNE domain, i.e. about 20-30% of the period which is what was observed in Froment et al. (2015).  Figure 4 of Paper I). On the right of the 2D plots, we display the evolution of respectively T e , n e and v around the loop apex (mean value between the two dotted bars in the bottom panel). On the bottom of the 2D plots, we show two profiles (solid and dashed lines, corresponding respectively to the hot phase at t 1 and the cool phase at t 2 , indicated by the solid and dashed lines on the right plots. Note that t 1 and t 2 are different for each simulations. For the velocity, the red (positive) is for flows from the eastern footpoint to the western one and the opposite for blue. On the 2D plots and the loop profiles, n e is shown in logarithmic scale. However, the apex time series n e are in linear scale.   Figure 1. We trace as well the evolution of the intensity in a portion around the loop apex (looking at the profiles along the contour), and the intensity profiles in the same way as for the simulations. Note that the range of intensities displayed is not the same as for the profiles from the simulations. ration/condensation cycles (giving thus a period close to the one detected for event 1 in Froment et al. 2015). These simulations are extracted from Figure 6 and thus correspond to H 1 = 12.8 × 10 −5 W m −3 . They are indicated by a red dot in Figure 6. They all have similar heating conditions. However, and as discussed earlier, these simulations are repeated using the unmodified Spitzer conductivity and 100,000 mesh points, as in Paper I. The first simulation is an incomplete condensation case for which λ 1 = 40.4 Mm and λ 2 = 33.0 Mm. Note that this simulation is not the same as the one presented in Paper I (with λ 1 = 50 Mm and λ 2 = 20 Mm) that has the same loop geometry. The two other simulations present complete condensations, with different locations of the condensations. For the middle simulation: λ 1 = 44.0 Mm and λ 2 = 29.4 Mm, and the condensations form close to the apex. For the last simulation: λ 1 = 40.4 Mm and λ 2 = 29.4 Mm, and the condensations form closer to the eastern footpoint. At t 1 , we display the loop profiles when the temperature reaches a maximum at the apex for one of the cycles, and at t 2 the profiles when the temperature reaches a local minimum. These three loops have a maximum apex temperature around 3 MK. During the cooling phase, when the condensations are established, the temperature drops to ∼ 1 MK in the eastern leg of the IC simulation, while the density increases by a factor of ∼ 2. For the CC simulations, the temperature drops locally to 0.01 MK and the density increases by a factor of 10.
We also notice a larger velocity for CC (up to 140 km s −1 ) compared to the IC case (about 10 km s −1 at the footpoints), probably due to the increase of the density of the condensation that falls compared the density of the loop itself. The velocity is also higher when the CC start closer to the apex, due to the longer acceleration time. We notice periodic siphon flows for both the IC and the CCs cases. The ones for CCs being stronger.
As already indicated before, around the apex, the amplitudes of the temperature (and even the density) evolution are not dramatically different for the complete and incomplete cases. The complete condensation, which is triggered in the eastern leg, outside of the apex area, has a minimal effect on the temperature evolution around the apex. It is thus not possible to use the evolution of the parameters at the apex to distinguish between complete and incomplete cases.
In Figure 17, we trace the evolution of the EUV intensity as it would have been seen in the 171 Å, 193 Å, and 335 Å channels of AIA (see Section 3.2.1 of Paper I for computation details), for these three loop systems (noted as IC, CC 1 and CC 2). As in Paper I, we used the AIA response functions to isothermal plasma for each channels, calculated with CHI-ANTI version 8.0 (Del Zanna et al. 2015). For comparison, we add the same plot from the AIA observations. The intensity is given along the loop defined by a smoothed version of the orange contour 5 displayed in Figure 1. However, as we already pointed out in detail in Paper I, the synthetic intensity can only be examined in the coronal part of the loop (due to the limitations of the model, we exclude the chromosphere and the low-transition region of the intensity analysis). We will discuss in more detail the intensity variation and values along the observed contour in the next section.
The overall pulsating behavior is well reproduced in the three simulations. In both complete and incomplete cases we can find the same global cooling pattern, with the intensity peaking first in 335, then 193 and finally 171, following the order of the peak responses of the channels. Note that we choose these three simulations because they are showing condensations to the eastern footpoint and thus matching the intensity pattern (higher intensities close to the eastern footpoint) along the observed loop bundle. The simulation presented in Paper I was showing condensations to the western leg of the loop. This heating case was showing the most convincing intensity light curves at the apex (and directly comparable to the observed light curves) among the simulations of the parameter space explored. However the IC case chosen in the present paper is showing a more convincing pattern along the loop for the reason detailed above (asymmetry of the intensity between the two loop legs).
Between the IC and CCs simulations, the biggest differences 6 occur at the location of the triggering of the condensations, where we can see another 335 peak, corresponding to the 0.2 MK peak of that band 7 . However, it seems quite challenging to look for this smaller peak in AIA data, as it is probably hidden in the line-of-sight integration to distinguish between CC and IC cases.
The pulsations that we observed in Froment et al. (2015) may also include co-spatial and simultaneous coronal rain events. However, it is quite difficult to distinguish between complete and incomplete cases using only the coronal channel of AIA, as discussed in Froment et al. (2017). It remains also difficult to conclude firmly for on-disk observations. However, the time lag between the 171 and 131 channels can help to identify the nature of the condensations even if we only have access to the mean behavior of the loop bundles. In Auchère et al. (2018), the authors detect long-period EUV pulsations coincident with coronal rain, in a region ob-5 Note that the length of the observed loop is then a bit shorter that the length of the simulated loops (derived form LFFF extrapolations). It does not affect our analysis as we discard the synthetic intensities from the simulated footpoints. 6 We notice also smaller intensity values for the cooling phases of the CCs cases than for the ones of the IC case but with only one loop it can delicate be to focus on absolute values. 7 The peak at lower temperature (O III to O V lines) has been accounted in the AIA response function since February 2013 (version 4), following the measurement of Soufli et al. (2012) served off-limb. In this study the 131 intensity peaks after 171, which was not the case for the events studied on-disk in Froment et al. (2015). We found no time lag between these channels, which indicates that the temperature of the plasma did not decrease on average below the peak response of 171 (around 0.8 MK).

4.2.
Are all the results of these simulations realistic?
We have previously seen that some of the TNE cases can reproduce very well the average behavior observed with AIA, in case of long-period intensity pulsations (see also the results of Paper I). Looking beyond the TNE cases, we can ask ourselves if the non TNE cases produced in the parameter space are realistic. Only a few simulations are hydrostatic and most of the non-TNE simulations are dominated by continuous siphon flows lasting for most of the three days of the simulations. For these case the simulated intensity would not show any temporal variations, which is inconsistent with observed EUV loops. In this regard, we have to bear in mind that in this simplified analysis we have modeled the average behavior of loop systems, and have not included the variability of the heating that is likely to occur in the corona. Further studies, including an exploration of the temporal variations of the heating and/or loop geometry are needed to quantify the dynamics of these systems and their stability.
One way to check if the densities produced by our simulations are consistent with observations is to compare the observed AIA intensities with the simulated ones. In Figure 17, we display the intensity along loop B for three different simulations, as it would be seen with AIA, considering a radius of the cross-section of the loop bundle of 100 km at s = 0. The intensity values in DN s −1 are between 0.1 and 10, except during the cooling phases. Considering that the loop intensity is about 10% above the background emission (Del Zanna & Mason 2003;Viall & Klimchuk 2011), this is consistent with the intensity counts derived froms AIA observations (between about 1 and 100 DN s −1 ) displayed in the same Figure. The fact that we model only one loop explains as well why we can easily identify the condensations in the profiles. In contrast, in the AIA observations the difference of intensity between the heating and cooling phase profiles chosen at t 1 and t 2 is quite small. The signature of the condensations is probably hidden by the background and foreground emission.

SUMMARY
In this paper, we explored a large range of dynamics, scanning different regimes of thermal non-equilibrium (TNE) and other behaviors in coronal loops. Several parameter-space studies regarding TNE cycles have already been conducted (e.g., Müller 2004;Susino et al. 2010;Mikić et al. 2013). Our study takes into account the recent discovery that long-period intensity pulsations are commonly observed in coronal loops.
The 1D hydrodynamic description of loops we used allows us to rapidly scan the parameter space. The model presented is rather simple, but captures the highly nonlinear dynamics of coronal loops. Indeed, we are able to nicely summarize the thermodynamic evolution of loops, even though the transition region and chromospheric behavior cannot be examined in detail.
For this extensive study we chose to explore a broad range of heating configurations, without explicitly limiting the heating profiles to a function of the magnetic field strength. We present in this paper a subset of this study, showing the results of 1,020 simulations.
We found TNE events in specific regions of the parameter space explored. With the different loop geometries (one semi-circular and two from LFFF extrapolations) used for the heating parameter scan, we conclude that: • Any loop geometry seems suitable for a loop system to undergo TNE cycles.
• However for each loop geometry the heating requirements to obtain TNE cycles are not the same.
• A stratified heating is a necessary condition but it is not sufficient to produce TNE. For each loop geometry, the heating parameter domain where we obtain TNE is different.
• The domain where we find TNE in the heating parameter space is limited.
• The more the heating is important at the footpoints the more the loop is likely to undergo TNE cycles and in particular complete condensations (CCs), rather than incomplete condensations (ICs), due to the high density of the plasma injected in the loop from chromospheric evaporations.
These conclusions might at first sight imply that any loop system could undergo condensation and evaporation cycles. However, this is not the case. In reality, the geometry and heating conditions vary from point to point. For a given loop, only one set of heating parameters exists. TNE is only possible when there is a specific match between the loop geometry and the heating conditions. Indeed, the long-period intensity pulsations reported by Auchère et al. (2014) and identified as TNE cycle by Froment et al. (2015); Froment et al. (2017) are widely observed in the corona but not in every loop bundle. There are probably many more cases of such cycles in loops, with heating conditions that change too much over time, producing more limited and irregular cycles. The Auchère et al. (2014) technique was designed to detect regular intensity pulsations, and is thus most sensitive to TNE events with stable pulsations.
The detection of the probably more frequent cases in which coronal conditions evolve with time would require a different method.
Our work present several limitations, in particular: simple input heating, poor treatment in the chromosphere and the transition region, simulation of a single loop, no time dependence of the heating. However, it aims to be a first step towards the exploration of the complex parameter space we only merely touched upon. More parameters could play an important role to trigger and maintain evaporation and condensation cycles. Eventually, elaborate simulations, possibly multi-dimensional, with a proper forward modeling could help to constrain the heating of the loop observed by comparing their behavior (cycles or not, period, time lag between the temperature and the density evolution, ...) with the results of such simulations.
This extensive parameter space study allowed us as well to explore some characteristics of the TNE events. These characteristics are summarized in Table 1 and Figures 13 to 15. We found that the period (from 2.4 to 15.5 hrs) is increasing with the length of the loop and with the maximum temperature reached. These periods tend also to be longer for CC compared to IC for the same loop geometry. The time delay between the temperature and density evolution, characteristic of TNE events when combined with the periodicity, is constant to within 20-30% of the period for most of the simulations (strong CC cases are an exception). We found as well that some loop geometries are more favorable to CC cases (see loop C). Moreover, looking at IC and CC cases in more details we show that both are showing siphon flows during the cooling phases. For CCs these flows are stronger. This is consistent with 2.5D magneto-hydrodynamic simulations of Fang et al. (2013Fang et al. ( , 2015. To conclude, we presented a unified picture of numerical simulations of cooling/heating in loops. We reaffirm in particular that coronal rain and long-period intensity pulsations are two manifestations of the same phenomenon, as demonstrated observationally by Auchère et al. (2018).
This work is an outgrowth of the work presented during the VII Coronal Loop Workshop and at Hinode 9. The authors acknowledge useful comments from attendees of these conferences. The authors would like to thank Jim Klimchuk and Patrick Antolin for fruitful discussions on thermal nonequilibrium and long-period pulsations in loops. The SDO/AIA and SDO/HMI data are available by courtesy of NASA/SDO and the AIA and HMI science teams. This work used data provided by the MEDOC data and operations centre (CNES / CNRS / Univ. Paris-Sud), http://medoc.ias.u-psud.fr/. Zoran Mikić was supported by NASA Heliophysics Supporting Research Grant NNX16AH03G. This research was supported by the Research Council of Norway, project number 250810, and through its Centres of Excellence scheme, project number 262622.