© The Authors 2023

1. Introduction

Solar prominences or filaments are cool and dense structures in the solar corona (Tandberg-Hanssen 1995; Vial & Engvold 2015). As the surrounding corona is tenuous, they can easily be observed in multiple wavelengths, such as Hα and extreme ultraviolet (EUV). This is convenient for detecting waves or oscillations in the solar corona, for instance, as well as for examining how they are manifested in prominences (Arregui et al. 2018; Chen et al. 2020).

Due to the magnetic buoyancy, heavy prominences can be suspended high up in the solar corona. There have been several mechanisms proposed to explain how such dense and cool plasma can enter the solar corona (see Mackay et al. 2010, for example). Meanwhile, many of these mechanisms have been demonstrated either by works of observation (e.g., Chae 2003; Berger et al. 2011; Zou et al. 2016) or simulation (e.g., An et al. 1988; Kaneko & Yokoyama 2015; Fan 2018). The evaporation-condensation mechanism, which has been studied in 1D geometry in detail in a number of previous works (see e.g., Antiochos et al. 1999, 2000; Xia et al. 2011), is among these mechanisms. Later on, this approach was extended to multi-dimensional magnetohydrodynamic (MHD) models (see e.g., Xia et al. 2012; Xia & Keppens 2016; Li et al. 2022).

This evaporation-condensation mechanism works in the following way. Cool plasma in the chromosphere, mainly from the higher chromosphere to lower transition region (TR) heights, enters the solar corona through evaporation. The heating that causes such evaporation is usually modeled as a parametrized localized heating source. It is denoted as H_loc so that it may be distinguished from the background heating, H_bgr, which is used to balance the coronal radiative cooling and maintain a hot corona. In the simulations, taking Zhou et al. (2021) as an example, the following expression is used for the background heating, H_bgr:

$$\begin{array}{c}\hfill H_{\mathrm{bgr}}=H_{0}exp(-(y-y_0)/{\lambda }_0),\end{array}$$(1)

where y is the vertical coordinate perpendicular to the solar surface and λ₀ is the heating scale height. Control parameters y₀ and H₀ serve to set the magnitude of the background heating. The expression in Eq. (1) is typically used in simulations which lack detailed lower-lying convection aspects and/or neglect adequate treatments of the full radiative transfer equation. Its justification is based on the observational fact that the total energy transported into the corona (per unit area) should be in the order of ∼2 × 10⁵ erg cm⁻² (Withbroe & Noyes 1977; Aschwanden 2001). Recently, Brughmans et al. (2022) studied the effect of other observationally motivated forms of background heating (see Mandrini et al. 2000, for example), especially on the in situ thermal-instability-driven formation process of prominence structures in flux ropes (while ignoring condensation-evaporation aspects).

The localized heating term H_loc usually takes the following form:

$$\begin{array}{c}\hfill H_{\mathrm{loc}}=H_{1}exp(-{|y-y_1|}^s/{\lambda }_1^s),\end{array}$$(2)

where s = 1 or 2. Although this expression is similar to that of H_bgr, in this case, H_loc behaves quite differently. From previous parameter surveys (Xia et al. 2011; Johnston et al. 2019; Pelouze et al. 2022), successful simulations of evaporation-condensation have taken H₀ to be on the order of ∼10⁻⁵ − 10⁻⁴ erg cm⁻³, while H₁ is typically larger by two to three orders of magnitude (H₁ ∼ 10⁻² − 10⁻¹ erg cm⁻³). The scale height, λ₀, has a typical value of ∼10² Mm or even higher, making the heating in the corona nearly a constant. While λ₁ is usually much smaller, making the heating decay very quickly when the height y deviates from y₁. That is why H_loc is called localized heating. In previous 1D simulations, λ₁ could be taken from a range 1%−10% of the total length of the loop, which is typically 2−20 Mm for a 200-Mm-long loop (Müller et al. 2004; Klimchuk et al. 2010; Xia et al. 2011; Pelouze et al. 2022). However, in multi-dimensional simulations, a narrower range of values ∼2 − 3 Mm is usually adopted (Keppens & Xia 2014; Xia & Keppens 2016; Zhou et al. 2020; Li et al. 2022; Jerčić & Keppens 2023). It should be noted that although the magnitude parameter H₁ is larger than H₀ by orders of magnitude, the total energy of H₁ integrated throughout the chromosphere and corona could be even one order smaller than H₀. This could be the reason why we still have no strong observational evidence for this localized heating.

The observable manifestation of the H_loc occurrence is, thus, still under debate. For example, random footpoint motions can lead to small-scale reconnections (nanoflares, as suggested by Parker 1988) and/or Ohmic dissipation (Pontin & Hornig 2020), more recently evidenced as “campfires” (Berghmans et al. 2021). Alfvén waves generated at the photosphere can also cause wave energy leakage into the chromosphere or transition region, through a variety of mechanisms (see, e.g., De Groof & Goossens 2002; van Ballegooijen et al. 2014; Howson et al. 2019). Torsional motions or tornado-like features may be responsible for energy transfer from the lower atmosphere (Su et al. 2012; Wedemeyer-Böhm et al. 2012). Small-scale energy events should in any case be ubiquitous, but are hard to detect directly in observations in short wavelength bands. Recent observations in non-thermal meter-wave radio channels provide some supporting evidence (Sharma et al. 2022).

Apart from studying formation of prominences, their oscillations have also widely received attention in recent years (Ballester 2005; Tripathi et al. 2009; Arregui et al. 2018; Chen et al. 2020). Prominence oscillations are ubiquitous, and they have been observed since the 1930s (Dyson 1930). At that time, filaments were found appearing and disappearing periodically in the Hα line center and wings when their line-of-sight (LOS) velocity is high enough. As a result, this phenomenon has also been dubbed “winking filaments” and it has, in fact, been frequently observed during an earlier time when spectroscopic observations were more popular, due to the lack of spatial resolution for imaging observations (e.g., Hyder 1966; Ramsey & Smith 1966). From these observations, it has been deduced that the period of prominence oscillations should be a sort of intrinsic property.

Starting from the beginning of this century, with the development of instruments, more and more winking filaments have been reported. Some of the winking filaments were reported as by-products of Moreton waves or EUV waves. For example, Eto et al. (2002) observed a winking filament caused by a Moreton wave using Hα line center and ±0.8 Å wings and, similarly, Okamoto et al. (2004) found a winking filament triggered by an EUV wave. Winking filaments could also result from lower atmosphere reconnections or coronal shocks (Isobe & Tripathi 2006; Grechnev et al. 2014).

Gilbert et al. (2008) proposed that the occurrence of the winking filament phenomenon requires LOS velocities of at least 30 − 40 km s⁻¹, which is related to the intensity contrast between the filament and the background. With only a few selected wavelengths available, instead of having high spectral resolution data, the LOS velocity could only be estimated using approximations. For example, Isobe & Tripathi (2006) observed a filament oscillation with a period of 2 h. Although this was mainly a horizontal oscillation, they used the method mentioned in Morimoto & Kurokawa (2003) obtaining LOS velocities of 20−30 km s⁻¹, which are much higher than the value of 4 km s⁻¹ for the horizontal direction. Similarly, Gilbert et al. (2008) observed a winking filament with a period of 29 min and a maximum LOS velocity of 41 km s⁻¹. However, Jackiewicz & Balasubramaniam (2013) revisited this event with another method and considered that the maximum LOS velocity is only few km s⁻¹. In Table 1, we list several parameters of winking filament obtained from selected observations.

We point out that according to the definition of a winking filament, only the LOS component of velocity is responsible for the winking. For a filament near the center of the solar disk, this would be an oscillation basically perpendicular to the solar disk, which is in line with intuitive expectation. However, for a filament near the solar limb, it should be predominantly an oscillation parallel to the solar disk. As demonstrated by our previous 3D MHD simulations of oscillating prominences (Zhou et al. 2018), the restoring force (and the associated period) for these two kinds of oscillations would be rather different. Therefore, in the present work, the term “winking filament” only refers to an oscillation perpendicular to the solar disk. In reality, however, solar filaments are expected to be in a rather dynamic equilibrium state over most of their lifetimes (see, e.g., Berger et al. 2010, 2017; Xia & Keppens 2016; Zhou et al. 2020; Jenkins & Keppens 2022; Jerčić & Keppens 2023). Thus, all kinds of dynamics (perpendicular and parallel) are likely to be coupled with each other.

In this work, we use the simulation to show that the formation process of solar filaments by evaporation-condensation mechanism can actually be responsible for the curious phenomenon of the winking filament. Our novel idea invokes an observationally justifiable temporal variability in the localized heating, and we show how the magnetic stretching in a 2D setup then leads to the winking in synthetic spectroscopic views. The paper is organized as follows. In Sect. 2, we discuss the setup of the simulation and the phenomenon of magnetic stretching; Sect. 3 presents a simulation of a winking filament. Section 4 gives our conclusion and a discussion.

2. Evaporation-condensation and magnetic stretching

In this section, we use a 2D model to demonstrate that a magnetic stretching phenomenon is embedded in the traditional evaporation-condensation mechanism. This is a purely multi-dimensional effect that can never be studied with the restricted 1D fixed-field assumption that is often made in evaporation-condensation scenarios.

2.1. Numerical setup

The 2D magnetohydrodynamic (MHD) model used here is similar in setup with previous work (Keppens & Xia 2014). The simulation box ranges from −50 Mm < x < 50 Mm and 0 < y < 80 Mm. We used a 768 × 768 uniform grid so that the resolution was 130 km × 104 km.

The atmosphere used in our simulation is composed of pure hydrogen. For the region below a certain height y_c = 2543 km, we adopted the temperature profile T(y) from the traditional VAL-C model (Vernazza et al. 1981). However, for the region above y_c, the following expression is used to extrapolate the semi-empirical VAL-C model so that the vertical heat flux is constant:

$$\begin{array}{c}\hfill T(y)={(7F_{\mathrm{c}}(y-y_{\mathrm{c}})/\left(2{\kappa }_{\Vert }\right)+{T_{\mathrm{tr}}}^{7/2})}^{2/7}.\end{array}$$(3)

We set the temperature at y_c, which is the typical temperature at the TR, T_tr = 4.47 × 10⁵ K so that the two profiles are continuously connected at y_c. The constant vertical thermal conduction F_c is 2 × 10⁵ erg cm⁻² s⁻¹ and the Spitzer-type conductivity κ_∥ is 8 × 10⁻⁷T^5/2 erg cm⁻² s⁻¹ K⁻¹ in this work. The atmosphere is partially ionized. The ionization degree of hydrogen is approximated by a single function of temperature (Heinzel et al. 2015). The technical details will be given in the appendix. Then, together with hydrogen number density at the bottom of the computational domain n_Hb = 9.45 × 10¹³ cm⁻³, we can setup the hydrostatic equilibrium atmosphere. This specifies both the density and the pressure profile (in 1D along the height).

After setting up the atmosphere, we added a potential, quadrupolar field that has a magnetic dip located at the horizontal center of our 2D domain:

$$\begin{array}{cc}& B_x=+B_{\mathrm{p}0}cos\left(\frac{\pi x}{2L_0}\right)e^{-\frac{\pi y}{2L_0}}-B_{\mathrm{p}0}cos\left(\frac{3\pi x}{2L_0}\right)e^{-\frac{3\pi y}{2L_0}},\hfill \end{array}$$(4)

$$\begin{array}{cc}& B_y=-B_{\mathrm{p}0}sin\left(\frac{\pi x}{2L_0}\right)e^{-\frac{\pi y}{2L_0}}+B_{\mathrm{p}0}sin\left(\frac{3\pi x}{2L_0}\right)e^{-\frac{3\pi y}{2L_0}}.\hfill \end{array}$$(5)

We set L₀ = 50 Mm while B_p0 is chosen to be 2 G in the demo case so that the plasma β near the TR is approximately unity and the magnetic field strength near the magnetic dip is approximately 4−5 G.

The governing equations used in this purely 2D (all vectors have only x, y components) simulation are as follows:

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \mathit{v})=0,\hfill \end{array}$$(6)

$$\begin{array}{cc}& \frac{\partial (\rho \mathit{v})}{\partial t}+\mathrm{\nabla }·(\rho \mathit{vv}+p_{\mathrm{tot}}\mathit{I}-\frac{\mathit{BB}}{{\mu }_0})=\rho \mathit{g},\hfill \end{array}$$(7)

$$\begin{array}{cc}& \frac{\partial e}{\partial t}+\mathrm{\nabla }·\left(e\mathit{v}+p_{\mathrm{tot}}\mathit{v}-\mathit{BB}·\mathit{v}\right)=\rho \mathit{g}·\mathit{v}+\mathrm{\nabla }·\left(\mathit{\kappa }·\mathrm{\nabla }T\right)\hfill \\ \hfill & -n_{\mathrm{e}}{n}_{\mathrm{H}}\mathrm{\Lambda }(T)+H,\hfill \end{array}$$(8)

$$\begin{array}{cc}& \frac{\partial \mathit{B}}{\partial t}+\mathrm{\nabla }·(\mathit{vB}-\mathit{Bv})=\mathbf{0}.\hfill \end{array}$$(9)

Here, $\mathit{g}=-g_{\odot }{r_{\odot }}^2/{(r_{\odot }+\mathit{y})}^2{\widehat{\mathit{e}}}_{\mathit{y}}$ is the gravitational acceleration and we have g_⊙ = 274 m s⁻², r_⊙ = 691 Mm. The expression ∇ ⋅ (κ ⋅ ∇T) is the field-aligned Spitzer-type anisotropic thermal conduction where κ is tensor defined as $\kappa ={\kappa }_{\Vert }\widehat{\mathit{b}}\widehat{\mathit{b}}$ and n_en_HΛ(T) is the optically thin radiative cooling term taken from the tables from Dalgarno & McCray (1972, for T < 14 000 K) and Colgan & Abdallah (2008, for T > 14 000 K). For the region below y = 2.5 Mm, this radiative cooling term is turned off because, in the lower layers, we would need a more consistent radiative transfer treatment for the chromosphere. The heating term, H, is given later. All the other symbols in the equations have their usual meanings.

For boundary conditions, symmetric settings are adopted for ρ, v_y, p, and B_y, while v_x and B_x are antisymmetric, on side boundaries. All the variables are fixed on the bottom boundary except that v_x and v_y are asymmetric. For the top boundary, B_x and B_y are extrapolated, v_x and v_y are asymmetric, while ρ and p are calculated according to a hydrostatic equilibrium assumption based on an extrapolated T. Other numerical settings, including the scheme, adopted limiter and divergence cleaning method, are the same as in our previous work (Zhou et al. 2021), except the transition region adaptive conduction (TRAC) method, which is not used in this work. Simulations are done with our open source code MPI-AMRVAC¹ (Xia et al. 2018).

2.2. Magnetic stretching

We first relaxed the system till t_relax = 214.7 min to a thermodynamic equilibrium state. During the relaxation, only a background heating term, H = H_bgr, was imposed to balance the radiative cooling. A time-dependent version of Eq. (1) is:

$$\begin{array}{c}\hfill H_{\mathrm{bgr}}(y,t)=R_{\mathrm{bgr}}(t)H_{0}exp(y/{\lambda }_0),\end{array}$$(10)

where R_bgr(t) = max(5−12t/t_relax,1) is used to setup a steady background heating after t = t_relax/3. We chose H₀ to be 6 × 10⁻⁵ erg cm⁻³ and λ₀ is 100 Mm.

The black solid lines in Fig. 1a show the topology of the magnetic field with some selected field lines, at the end of the relaxation stage t = 214.7 min. The color scale shows the temperature distribution at this moment. To show these distributions more clearly, a vertical slice is taken along the y-axis at x = 40 Mm. Temperature and number density (of hydrogen) distributions along this slice are shown in Fig. 1b, with blue and red solid lines, respectively. We can see that the TR is located a little bit higher than y = 2 Mm.

Fig. 1. Initial condition and magnetic stretching. (a) Temperature distribution and magnetic configuration after relaxation. (b) Temperature and number density distributions along x = 40 Mm, as indicated by the dashed line in panel a. (c) Selected magnetic field line at t = 100.2 min and after relaxation. (d) Selected magnetic field line before and after we impose the localized heating.

Actually, the system is already in a relatively stable state after about 100 min, as shown in the later analysis. However, we continue to relax the system for around 100 min to make a clear comparison between the relaxation stage and the further localized heating stage. To show how a magnetic field line indeed no longer changes in the latter half of the relaxation phase, we select the field line which starts from (x, y) = (45, 0) Mm (so that it ends at (x, y) = (−45, 0) Mm due to its symmetry) and plot it in Fig. 1c. In that panel, we compare the trace of this field line at t = 100.2 min (blue dashed line) and at t = 214.7 min (red solid line). Clearly, we can see that during the extended relaxation stage (without localized heating, H_loc), this field line can stay at the same position for more than 100 min.

Then we start to impose the localized heating term H_loc. The localized heating term is a 2D version of Eq. (2) together with a modulating ramp function, R_ramp, so that the localized heating can increase smoothly from zero after t_relax:

$$\begin{array}{cc}& H_{\mathrm{loc}}(x,y,t)=H_1{R}_{\mathrm{ramp}}(t)H_x(x)H_y(y),\hfill \end{array}$$(11)

$$\begin{array}{cc}& H_x(x)=exp\left(-\frac{{(x-x_r)}^2}{{\sigma }^2}\right)+exp\left(-\frac{{(x-x_l)}^2}{{\sigma }^2}\right),\hfill \end{array}$$(12)

$$\begin{array}{cc}& H_y(y)=exp\left(-\frac{{(y-y_1)}^2}{{\lambda }_1^2}\right)·\hfill \end{array}$$(13)

Here, H₁ = 2 × 10⁻² erg cm⁻³ s⁻¹, y₁ = 4 Mm and λ₁ = 3.16 Mm. In the x-direction, we use the parameters x_l = −41.5 Mm, x_r = 41.5 Mm, and σ = 5.48 Mm to ensure that the heating is concentrated at footpoints. We simply take a periodic sine function for the ramp function R_ramp:

$$\begin{array}{c}\hfill R_{\mathrm{ramp}}=max(sin(2\pi (t-t_{\mathrm{relax}})/P),0).\end{array}$$(14)

The period P is taken as 20 min, a typical value listed in Table 1 and also a typical value from the our previous simulation work on vertical oscillations (Zhou et al. 2018). Then the first heating pulse will end in 10 min.

After the first 10 min, when t = 224.7 min, as shown in Fig. 1d, we can see that the field line we showed earlier gets stretched significantly (black dashed line), compared to its shape from 10 min before (red solid line). Quantitatively, the apex of the loop rises from 3.15 Mm to 3.64 Mm, with an increment of 16%. The inferred velocity of such a rise during these 10 min is on average 8.2 km s⁻¹.

In previous works, such a vertical stretching of the magnetic field line is usually considered to be associated with eruptive events, typically erupting solar filaments and/or magnetic flux ropes (see, for example, Chen et al. 2002; Amari et al. 2014; Jiang et al. 2021), in which plenty of energy is released. However, in our numerical experiment shown above, even a small portion of gradual energy injection into the system can lead to a significant stretching. We will look into the details of the physics here in Sect. 4.

Our 2D MHD simulation shows that this kind of evaporation-induced field line stretching, where the increased thermal pressure gives rise to an overall inflated magnetic topology, should be ubiquitous on the Sun. With such a natural multi-dimensional mechanism, we can now propose an alternative explanation for the winking filament phenomenon, that would be entirely driven by the popular evaporation-condensation scenario for prominence formation.

3. Winking filament

In Sect. 2, we have shown that the magnetic field lines in an arcade setup would be stretched during an evaporation-condensation process. However, we cannot actually see the configuration of magnetic field lines in observations. Luckily, there are at least two ways of how to “see” this stretching indirectly by observing the motion of a filament in the arcade that is supported by the field line or by observing a possible wave generated by this stretching. We go on to demonstrate that this can indeed induce the vertical oscillation of the filament, thereby leading to the winking filament.

As we did not yet have a filament in our setup at t = 224.7 min, we continue our simulation from the point where we stopped in Sect. 2. We still adopt the periodic heating from Eq. (11) with the period of P = 20 min. As shown in Johnston et al. (2019), adopting different periods and durations of the added impulsive heating will lead to different scenarios involving thermal non-equilibrium evolutions and/or condensations driven by thermal instabilities. For the typical period range listed in Table 1, this influence should be minor and we expect to be able to form a large-scale prominence.

Considering the condensation-evaporation mechanism, it is expected to cause a condensation, when localized heating is added for a limited period of time. This first condensation grows larger and gradually forms a filament. Figures 2a–c show the temperature evolution during the filament formation process in our simulation. In panel a, at t = 224.7 min, the average coronal temperature inside the black rectangle (−5 Mm < x < 5 Mm and 5 Mm < y < 35 Mm) reaches a maximum value of 2.23 MK, increased from the temperature of 1.55 MK before we imposed the localized heating. After that, the temperature (inside the chosen rectangle) starts to drop gradually. At around t = 430.8 min, catastrophic cooling due to thermal instability happens so that a clear condensation starts to appear in the simulation domain (panel b). Here, we define a condensation as cool material or filament as plasma below 14 000 K, because the values of the ionization degree listed in Table 1 in Heinzel et al. (2015) were calculated for temperatures from 6000 K to 14 000 K. Then, the area and mass of the cool filament increase gradually. Panel c shows the temperature distribution at t = 480.9 min, when we can clearly see a vertical sheet-like filament that formed. From the attached animation, we can already clearly see the oscillation in the vertical direction. At the same time, we can also see a periodic shrinking and expansion of the filament in the horizontal direction.

Fig. 2. Formation of the filament. (a)–(c) Temperature distributions during the formation of the filament at the times t = 224.7 min, t = 430.8 min and t = 480.9 min. An animation showing their evolution is available online. (d) Ionization fraction distribution at t = 480.9 min. (e) AIA 304 Å synthetic image at t = 480.9 min. (f) Approximate Hα synthetic image at t = 480.9 min. The units of the last two panels are arbitrary. See text for details.

Figure 2d gives the distribution of the ionization fraction for the time as is for panel c. Because the ionization fraction in Heinzel et al. (2015) only ranges from 0.17 to 0.94, we call here the plasma fully ionized when its ionization degree is near 0.94, instead of 1. Most of the hydrogen within the condensation region is only partially ionized, as expected. We can also see that most of the region within the outmost heated arcade is not fully ionized, though the temperature there is typically above 10⁵ K.

To make a direct comparisons with observations, we have made synthetic images in He II 304 Å and Hα, which are shown in Figs. 2e and f. It is noted that neither of these two lines is usually considered optically thin. However, for the typical environment of a quiet corona and filament, the optically thin assumption still works well in practice for He II 304 Å. This is shown in, for example, Chen et al. (2015) and Xia & Keppens (2016). Thus, we can simply calculate the synthetic 304 Å radiation using:

$$\begin{array}{c}\hfill I_{304}(x,y)=G_{304}(T(x,y))n_{\mathrm{e}}(x,y)n_{\mathrm{H}}(x,y)\mathrm{\Delta }z.\end{array}$$(15)

The temperature-dependent response function G₃₀₄(T) is obtained from the CHIANTI atomic database (Dere et al. 1997; Del Zanna et al. 2021), where the lower temperature limit is 10⁴ K. Therefore, for regions cooler than 10⁴ K in our simulations, we used the fixed value G₃₀₄(T = 10⁴ K) as an approximation. In addition, a thickness of Δz = 10 Mm is assumed along the LOS direction. From Fig. 2e, we can see that the heated arcade section is slightly brighter than the background, while the prominence–corona transition region (PCTR) that envelopes the filament is the brightest part in the corona. We note that in this 2D simulation, we cannot see the whole PCTR covering the filament, especially the part in the LOS (z) direction. That explains why the filament appears as an artificially dark structure. In fact, in a three-dimensional (3D) simulation, where we can see the whole PCTR, the filament would be a bright structure, as seen in observations.

The Hα image in Fig. 2f is synthesised using the approximate method in Heinzel et al. (2015), which is found to be a good approximation compared to a one-dimensional (1D) full radiative transfer approach (Jenkins et al. 2023). This method works well only for the filament in the corona. Thus, synthesised radiation in the lower region, where y < 5 Mm, may deviate from actual values. Again, a thickness of Δz = 10 Mm in the third direction is assumed. In this panel, we can clearly see a vertical sheet-like bright structure, representing the bright prominence.

In order to study the formation process in detail and to better show the oscillations, we picked out a particular region, namely, the region within the black rectangle in Fig. 2a and we studied the property changes inside it. Figures 3b–d show the time evolution of the average thermal pressure, $\overline{p}$, the average number density of hydrogen, $\overline{n_{\mathrm{H}}}$, and the average temperature, $\overline{T}$, over this region, together with the time evolution of R_ramp, as described in Eq. (14) shown in panel a. Dashed lines in these panels show the starting point of the localized heating and the starting point of the condensation, namely, t = 214.7 min and t = 430.8 min, respectively.

Fig. 3. Time evolution of (a) Rramp, (b) average pressure, (c) average number density, (d) average temperature in the rectangle domain shown in Fig. 2a. (e) and (f) Time evolution of the y-position and vertical velocity of the apex of the selected magnetic field line. (g) and (h) Time evolution of the y-position and vertical velocity of the mass center of the selected rectangle domain. Dashed lines show the starting time of localized heating (red) and starting time of condensation (blue).

With the evaporation proceeding, the density in the corona increases gradually. For the temperature and thermal pressure, these will first increase due to the evaporation, but then decrease due to the strong radiative cooling. Ten periodic oscillations of various amplitudes could be found in Figs. 3b–d, between the two dashed lines, corresponding to the parametrically imposed heating periods in panel a.

After condensation occurs, the number density increases in steps, with no clear period to be found (see Fig. 3, panel c). Meanwhile, the time evolution of thermal pressure and temperature (panels b and d) still show clear periods, but they are decreasing overall. The thermal pressure will drop to even lower than the value before we impose the localized heating. These results are similar to our previous 1D or 2D simulations (Xia et al. 2011; Keppens & Xia 2014).

According to the magnetic stretching mechanism mentioned in Sect. 2, this kind of periodic heating will cause a correspondingly periodic shrinking and expansion of the magnetic field lines. Thus, based on the MHD frozen-in theory, the coronal plasma as well as the filament will move together with the field line. Then, we would expect to detect the filament oscillation in the vertical direction, with the same period of the periodic impulsive heating.

To show the oscillation more clearly, we first pick out the selected field line again (as described in Sect. 2). The time evolution of the y-position of the apex of the selected field line, namely, y_pk(t), is shown in Fig. 3e. Panel f shows the corresponding velocity v_pk(t) in the y-direction at this point. The oscillation of the field line is slightly different from the physical parameters in panels b–d. Its oscillation period seems to be only half of the given P. And every two periods are composed of a larger peak and a smaller one. The oscillation of the selected field line has a typical amplitude of about 5 Mm and a velocity amplitude of about 20 km s⁻¹. And since the plasma around this field line does not experience condensation, the oscillation pattern does not change significantly after condensation occurs on the inner field lines.

The oscillation of the filament is described by the y-position of the mass center y_mc(t) of the rectangle region (see Fig. 2, panel a) as well as its corresponding density-weighted average velocity v_mc(t) in the y-direction. Their time evolution could be found in Figs. 3g and h, respectively. Similarly to the oscillation of the selected field line, the oscillation of the filament also has a period of only 1/2 P. It can be seen (especially from the y_mc curve) that the oscillation becomes stronger after condensation occurs. The amplitude of v_mc also increases slightly, with a typical value slightly higher than 10 km s⁻¹.

To demonstrate the oscillation of the filament more clearly, we took a vertical slice along x = 0, and made a time-distance plot in Fig. 4. The evolution of the number density of the filament, synthesised 304 Å radiation and Hα radiation along this slice are shown in panels a–c, respectively. From these panels, we can clearly see the oscillation during the formation and growth of the filament. In the number density panel, we also drew the y-position of the mass center using dashed lines. In this way, we can clearly see that the period of the oscillation is not P, but approximately 1/2 P. A similar dashed line, which is weighted by the Hα radiation, is shown in panel c, displaying a similar oscillation. In the case of the 304 Å panel, the bright PCTR structure shows the oscillation with period of approximately 1/2 P, similarly to the density and Hα plots.

Fig. 4. Time-distance plots along the y-axis for (a) number density, (b) 304 Å synthetic radiation, and (c) Hα synthetic radiation. Dashed lines in panels a and c trace the density-weighted and radiation-weighted center position, respectively, after the condensation.

We refer to the filament oscillating in the vertical direction as “winking” because it will appear and disappear periodically in the Hα line center and wings (especially). To confirm this, in addition to the oscillation in the prominence view (face-on view), we must check the filament view (top-down view). The synthetic Hα radiation in the filament view is also calculated following the method in Heinzel et al. (2015), for both Hα line center and wings at wavelengths ±0.8 Å. Observations of prominences and filaments by Eto et al. (2002) and Okamoto et al. (2004) were also made in these wavelengths. However, in that view, the chromosphere (the region below y = 2.5 Mm) is not included in this upward integration since the emission from the lower atmosphere is treated as background emission. Our simulation is a 2D simulation, and when observed from the top, the data cube will collapse into an integrated emission sheet. We show the time evolution of this Hα sheet in Fig. 5a.

Fig. 5. Time evolution of the synthetic (a) Hα line center (b) Hα red wing (c) Hα blue wing radiation from the filament view (top-down view). Vertical dashed lines are drawn to help to see the periods.

Before condensation occurs, nothing special could be observed. After the condensation, we can see a “black cloud” appear in the figure. The area of the Hα filament changes periodically, becoming sometimes small and faint while sometimes large and clear. The maximum area is typically 3−4 times bigger than the minimum area. Typically, in this case, the filament always exists in the Hα band. In some of our tests (see Sect. 4), it will be very faint and nearly invisible. In observations, both situations could occur (see also the references listed in Table 1).

Synthetic filament views observed in the red and blue wings are presented in Figs. 5b and c, respectively. Similarly to the previous oscillation patterns, the period here is typically 1/2 P. It is clear in these two panels that the filament is visible in the red wing and blue wing alternately, which is consistent with the observations. Thus, the simulated filament can be called a winking filament.

4. Discussion and conclusion

4.1. Details of the magnetic stretching

In this part of the discussion, we look to find an answer to the question of why magnetic field lines become stretched during the evaporation. The physical scenario should be as follows. When the lower atmosphere is heated by the localized heating, the local pressure will increase so that it can produce a pressure gradient, by which plasma will be pushed to move upward, generating upward flow. This is the so-called evaporation. This pressure gradient is typically directed vertically upwards (i.e., along the positive y-direction). However, the magnetic field lines are inclined, having an angle (x-varying) with the vertical y-direction. According to the MHD frozen-in theorem, the upward moving plasma will rise together with the magnetic field lines, dragging it into a more stretched state. We can illustrate the analogy with an inflated balloon: as pressure rises below, the line-tied magnetic field lines stretch outwards. We go on to analyze this in detail for a simple, but representative, ramped-up heating phase.

Therefore, instead of the periodic R_ramp as used in the previous sections, here we use a linear R_ramp function instead, whereby:

$$\begin{array}{c}\hfill R_{\mathrm{ramp}}(t)=\{\begin{array}{cc}4\left(t-t_{\mathrm{relax}}\right)/P& \mathrm{if}{t}_{\mathrm{relax}}<t<t_{\mathrm{relax}}+P/4,\\ 1& \mathrm{if}t\ge t_{\mathrm{relax}}+P/4.\end{array}\end{array}$$(16)

This new R_ramp(t) function is shown in Fig. 6a, and grows from 0 to 1 linearly in 5 min, that is, during t = 214.7 min and t = 219.7 min. We again selected the field line starting from (x, y) = (45, 0) Mm and analyzed the forces at a particular height on this line, for instance, at y = 4 Mm on the right side. During the simulation, the magnetic field line deforms and stretches, so that its position moves in the horizontal direction. Therefore, it cannot always be the “same” point. Considering that the deformation of the field line at the lower atmosphere is minor, we treated the selected point as the same point approximately, as indicated by the blue circle in Fig. 1d.

Fig. 6. Time evolution of Rramp shown in panel a. Time evolution of (b) |By/Bx|, (d) gas pressure gradient in the y-direction, (f) Lorentz force in the y-direction, (h) gravity force, (j) total force in the y-direction before the localized heating is imposed are shown on the left column. Panels in the right column show the time evolution of these parameters after the localized heating is on, correspondingly.

In the second row, we show the time evolution of the inclination of the magnetic field line at this point in the two steady heating phases, namely, before (during relaxation) and after the ramping up of the localized heating. In panel b, before we introduce the localized heating, namely, prior to t = 214.7 min, |B_y/B_x|, which indicates the inclination of the magnetic field line at the selected point, is typically 2.4. This increases a little bit after the heating is fully on (panel c), which means that the magnetic field line becomes more vertical. Hence, it gets stretched in the y-direction. In panel d the y-component of pressure gradient, namely, −(∇p)_y, is typically −3 × 10⁻¹⁰ in cgs units at the end of the relaxation phase, but increases to about 5 × 10⁻¹⁰ (panel e), showing a significant increase. Although it is not shown here, the x-component of the pressure gradient shows the same tendency. That means during the relaxation phase, the pressure gradient is slightly inwards (negative −(∇p)_x, y values) for the selected loop, but its direction turns outward in the evaporation phase. Correspondingly, the Lorentz force J × B changes in the opposite direction to ensure net overall force balance, from slightly directing outward to pointing inward (panels f and g). Meanwhile, the evolution of the gravity force (ρg)_y shows a similar response to the localized heating as the Lorentz force (panels h and i). This means that before the localized heating is imposed, the gravity together with the gas pressure gradient exerts a force to collapse the loop. Meanwhile, the Lorentz force plays a supporting role. In panel j, we can clearly see that the total force F_tot (in the y-direction) is nearly balanced during the relaxation. Then, after the extra localized heating is on, the gas pressure gradient pushes the plasma moving upward, trying to expand the loop, while the Lorentz force will suppress such a motion, and changes its direction to act downwards. Such a clear change in the interaction of all forces eventually causes the change in the shape (i.e., the local direction) of the magnetic field. We note that the evolution of the total force (panel k) is very close to the evolution of pressure gradient (panel e), indicating the most important role played by the pressure gradient here.

4.2. Parameter survey

Since the stretching of the magnetic field line is triggered by the evaporation flow caused by the upward pressure gradient, plasma β will definitely be important in this physical process. To see how different plasma βs would influence the result of our simulation, we can choose to use different magnetic field strength, namely, B₀, or alternatively, to put the localized heating at different positions. These two methods are not totally equivalent, but both of them should work. Here, for convenience, we chose to put the localized heating at different positions, to simulate heating at different heights. For all these different runs, we simply compare the heights of the apex of the selected field line. The results are shown in Table 2.

The localized heating typically takes place at the chromosphere and TR. Thus, we chose for y₁ in Eq. (12) to change from 2 Mm to 6 Mm, while the x-position (e.g., x_r) is changed correspondingly along the field line (x_l = −x_r). The corresponding plasma β at (x_r, y₁) ranges from 0.1188 to 0.1581, thus, it changes only slightly. However, we can see that in different runs, the vertical velocity v_y at y = 10 Mm (on the same field line connecting (x_r, y₁)) increases from 37.6 km s⁻¹ to 69.11 km s⁻¹. In addition, the apex of the selected field line y_pk (t = 243.3 min) also rises significantly. This means that when the same magnitude of localized heating is acting at a lower position in the atmosphere, the stretching effect becomes progressively minor. When it takes place in the higher chromosphere or TR, the stretching can work very effectively.

Besides plasma β, another important ingredient of this mechanism is the strength of the heating, H₁. It represents the energy release rate throughout the lower atmosphere. The results of different runs with different H₁ are listed in Table 3, similar with the values used in previous works (Xia et al. 2011; Pelouze et al. 2022). For localized heating injected into the system at the same location, stronger heating will generate higher upflow velocity and the magnetic field lines will get more stretched. However, while stronger heating will induce a more effective magnetic stretching, it will prohibit condensation from occurring, as demonstrated by Xia et al. (2011).

Notice should be taken that, the test runs here quantify the oscillation of field lines before the formation of filaments. Therefore, these runs cannot yet be compared directly to the observations, but if we were to continue these runs using similar periodic cycling, the same winking effects on the filaments (once formed) can be expected.

In this work, we focus on the evaporation-condensation mechanism. A similar mechanism for filament formation is the injection model, in which thermal instability is not necessary, and energetic events in the lower atmosphere will push directly the cool material into the solar corona. From the “unified model” for prominence formation (Huang et al. 2021), the difference between the evaporation-condensation model and the injection model should typically be achieved by varying the values of y₁ and H₁ here. The injection model usually adopts a lower y₁ and stronger H₁. Thus, in principle, the injection model could also produce such kind of vertical oscillations.

4.3. Driven period

The period of the impulsive heating P is chosen to be 20 min in this work. However, the period of the winking filament seems to be approximately only half of this value – and this is not a surprise. We should distinguish three types of periods, or frequencies here: the frequency of the impulsive heating or driven force, ω_dirv, the frequency, ω_wink, of the observed oscillation of the winking filament, and the eigenfrequency, ω₀, of the magnetic loop. In the research of coronal seismology, especially the kink mode, ω_dirv is usually equal to ω₀ (see recent reviews in e.g., Van Doorsselaere et al. 2020; Nakariakov et al. 2021). This is because the loop system will naturally select out a narrow band frequency close to the eigenfrequency, so that the kink mode can work efficiently.

However, in our simulation, the frequency of the additional heating ω_dirv is more likely to be inconsistent with ω₀. A similar study was recently carried out by Ni et al. (2022) where the authors use periodic jets to drive the oscillation of filaments. In a more general sense, this scenario is similar to the forced oscillators in the classic mechanics, where the solution usually takes the following form:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill s(t)& =\frac{f}{{\omega }_0^2-{{\omega }_{\mathrm{dirv}}}^2}\left(cos{\omega }_{\mathrm{dirv}}t-cos{\omega }_{0}t\right),\hfill \\ \hfill & =\frac{2f}{{\omega }_0^2-{{\omega }_{\mathrm{dirv}}}^2}sin\frac{\left({\omega }_0-{\omega }_{\mathrm{dirv}}\right)t}{2}sin\frac{\left({\omega }_0+{\omega }_{\mathrm{dirv}}\right)t}{2},\hfill \end{array}\end{array}$$(17)

where s(t) is the displacement of the oscillator and f is the driven force per mass. According to Eq. (17), the observed frequency ω_wink is composed of two frequencies: (ω₀−ω_dirv)/2 and (ω₀+ω_dirv)/2. We note that the amplitude of the oscillation is inversely proportional to $\left({\omega }_0^2-{\omega }_{\text{dirv}}{}^2\right)$. Therefore, ω_dirv cannot deviate from ω₀ too much. With this assumption, then, considering that we would usually recognize the higher frequency in observations, we have:

$$\begin{array}{c}\hfill {\omega }_{\mathrm{wink}}=({\omega }_0+{\omega }_{\mathrm{dirv}})/2.\end{array}$$(18)

We can take ω₀ = 2π/P₀ = v_Aπ/L, which is usually used in the study of coronal seismology (Roberts et al. 1984), as an approximation. Then, P₀ is the corresponding eigenperiod, v_A is the Alfvén speed, and L is the length of the loop. In our simulations, ω_dirv is fixed while ω₀ could be different from field line to field line and from time to time. However, before the formation of filament (i.e., between t = 224.7 min and t = 430.8 min), the eigenperiod P₀ within the region where our filament later forms is typically around 450 s and gradually drops to around 330 s after condensation. Using Eq. (17), we then get the observed period P_wink of the winking from the range 542 s–655 s, which is approximately 10 min or 1/2 P. This value is similar to our results.

Therefore, it is also not necessary to have the period of the impulsive heating of exactly 20 min, and even the impulsive heating does not necessarily need to be periodic. A quasi-periodic heating should be enough to get the winking phenomenon.

4.4. Conclusion

Within the one-fluid ideal MHD description, magnetic field lines move together with the plasma. Small energy events in the lower atmosphere can produce localized heating that increases the local pressure gradient, causing the evaporation of the plasma. Since the magnetic field is coupled with the plasma, the plasma can move freely along the field line, and at the same time, the field line will be dragged by any upward moving plasma. We have shown that at least a 2D (but completely analogous in 3D) model is required to show that line-tied stretching of field line is inevitable whenever a localized heating turns on, beyond the background coronal heating. This process has been overlooked by all restricted 1D models, which reduce coronal dynamics and prominence formation scenarios in strong field settings to 1D hydrodynamic along a fixed, prechosen field line shape.

In this work, we demonstrate that the evaporation-condensation scenario of prominence formation, during which a small portion of energy is injected into the low atmosphere, can produce oscillations of the filament in the vertical direction. When the heating source then acts (quasi-) periodically, the winking filament phenomenon can be observed. However, this also depends on the position and strength of the heating source.

In our simulation, we found that the periods of impulsive heating and the oscillations of the filament are equal to each other, and period of the oscillations is typically reduced by a factor of two. This must relate to the eigenfrequency of the perturbed loop being different from the frequency of the impulsive heating. The interplay and beating of two frequencies result in the altered period of oscillations that we see in our simulations.

In previous observations, a winking filament is usually interpreted as being caused by waves or energetic events. Here, we demonstrate through numerical simulations that winking filaments may also be observed due to the evaporation during the formation phase of filaments. As far as we know, this has not yet been reported in observations, but we expect that future observations will be able to fully confirm this scenario.

Movie

Movie 1 associated with Fig. 2 (movie) Access here

Acknowledgments

We thank A. Hillier, S. Gunár, M. Guo for valuable suggestions. We thank the referee for detailed suggestions on writing standard. Y.Z. acknowledges funding from Research Foundation – Flanders FWO under the project number 1256423N. X.L. and R.K. acknowledge the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No. 833251 PROMINENT ERC-ADG 2018). J.H. acknowledges funding from NSFC under grant 11903020. R.K. acknowledges support by Internal funds KU Leuven, project C14/19/089 TRACESpace and FWO project G0B4521N. Visualisations used the open source software Python (https://www.python.org/). Resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation – Flanders (FWO) and the Flemish Government.

References

Appendix A: Ionization Degree in MPI-AMRVAC code

In the original MPI-AMRVAC 2.2 code (Xia et al. 2018), the atmosphere is assumed to be fully ionized. To treat the Hα emission in a more self-consistent way, we include ionization degree (of hydrogen) in our simulation. Since the method is not described either in our previous work (Zhou et al. 2020) or the MPI-AMRVAC paper (Xia et al. 2018), we give the details in this Appendix.

In the MPI-AMRVAC code, conservative variables, for instance, density ρ, momentum ρv, energy (volume) density e and magnetic induction B, are solved. The quantity e is expressed as: e = e_int + ρv²/2 + B²/2μ₀. It is composed of internal energy e_int, kinetic energy ρv²/2, and magnetic energy B²/2.

When ionization is considered, we have the ionization fraction of i = n_e/n_H, where n_H = ρ/m_H is the total number density of hydrogen, m_H is the mass of a hydrogen atom, and n_e is the number density of electron or ionized hydrogen. The ionization fraction i is a function of vertical altitude y, thermal pressure p and temperature T; for instance, we have i = i(y, p, T), as tabulated in Heinzel et al. (2015).

Then, the internal energy e_int and gas pressure p could be written as

$$\begin{array}{cc}\hfill e_{\mathrm{int}}& =\frac{p}{\gamma -1}+i{n}_{\mathrm{H}}\chi ,\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill p& =(n_{\mathrm{e}}+n_{\mathrm{H}})k_{B}T=(1+i)n_{\mathrm{H}}{k}_{\mathrm{B}}T,\hfill \end{array}$$(A.2)

respectively. γ = 5/3 is the adiabatic index, k_B is the Boltzmann constant and χ = 13.6 eV is the ionization energy of hydrogen.

In the future 3.1 version of MPI-AMRVAC, the ionization degree, i, is will be function of temperature only. Such i(T) function has been tabulated, for instance, in Carlsson & Leenaarts (2012) and Hong et al. (2022). Then, we can rewrite Eq. A.1 and Eq. A.2 as:

$$\begin{array}{c}\hfill \frac{e_{\mathrm{int}}}{n_{\mathrm{H}}}=\frac{(1+i)k_{\mathrm{B}}T}{\gamma -1}+i\chi .\end{array}$$(A.3)

Considering that i(T) is monotonically increasing with T, the right-hand side of Eq. A.3 is monotonic. Therefore, T as well as i could be determined from e_int/n_H, according to the dependence of i(T). Therefore, after the update of e and ρ in the code, the variables T and i are then determined from e_int/n_H.

However, in the current work, i is not a single function of T, but instead it is dependent on y, p, and T. We cannot use the same simple way to calculate i and T from e_int and n_H. For convenience, approximately, we can drop the ionization energy term n_eχ on the right hand side of Eq. A.1. This approximation will overestimate the temperature a little bit especially when temperature is low, but still acceptable. Then, at a certain grid point, with known p and y, the ionization fraction, i, is again a single function of temperature, T. In this way, again, i and T could be obtained from n_H and e_int simultaneously.

Since the table given in Heinzel et al. (2015) is sparse, with only 3-6 points on y, p, and T, we used linear interpolation in between in our code. For values beyond the table, we simply use the nearest value. For example, the altitude y is from 10 Mm to 30 Mm in the table. Thus, for the region above 30 Mm, we used the value i(y = 30 Mm, p, T) directly, instead of extrapolating on the basis of the table. As a result, the ionization fraction i has a range from 0.17 to 0.94 in our simulation, instead of 0 to 1.

All Tables

All Figures

	Fig. 1. Initial condition and magnetic stretching. (a) Temperature distribution and magnetic configuration after relaxation. (b) Temperature and number density distributions along x = 40 Mm, as indicated by the dashed line in panel a. (c) Selected magnetic field line at t = 100.2 min and after relaxation. (d) Selected magnetic field line before and after we impose the localized heating.
In the text

Fig. 2. Formation of the filament. (a)–(c) Temperature distributions during the formation of the filament at the times t = 224.7 min, t = 430.8 min and t = 480.9 min. An animation showing their evolution is available online. (d) Ionization fraction distribution at t = 480.9 min. (e) AIA 304 Å synthetic image at t = 480.9 min. (f) Approximate Hα synthetic image at t = 480.9 min. The units of the last two panels are arbitrary. See text for details.

In the text

Fig. 3. Time evolution of (a) Rramp, (b) average pressure, (c) average number density, (d) average temperature in the rectangle domain shown in Fig. 2a. (e) and (f) Time evolution of the y-position and vertical velocity of the apex of the selected magnetic field line. (g) and (h) Time evolution of the y-position and vertical velocity of the mass center of the selected rectangle domain. Dashed lines show the starting time of localized heating (red) and starting time of condensation (blue).

In the text

	Fig. 4. Time-distance plots along the y-axis for (a) number density, (b) 304 Å synthetic radiation, and (c) Hα synthetic radiation. Dashed lines in panels a and c trace the density-weighted and radiation-weighted center position, respectively, after the condensation.
In the text

	Fig. 5. Time evolution of the synthetic (a) Hα line center (b) Hα red wing (c) Hα blue wing radiation from the filament view (top-down view). Vertical dashed lines are drawn to help to see the periods.
In the text

Fig. 6. Time evolution of Rramp shown in panel a. Time evolution of (b) |By/Bx|, (d) gas pressure gradient in the y-direction, (f) Lorentz force in the y-direction, (h) gravity force, (j) total force in the y-direction before the localized heating is imposed are shown on the left column. Panels in the right column show the time evolution of these parameters after the localized heating is on, correspondingly.

In the text