Tearing Instability and Periodic Density Perturbations in the Slow Solar Wind

In recent years, the usual dichotomy between a fast (>600 km s⁻¹) and a slow (∼350–400 km s⁻¹) wind has been challenged as the main discriminating factor to identify the source of the solar wind plasma. A more accurate picture could be described as follows: a rather steady component (often fast but not always), coming from open regions, or coronal holes, is associated with a well-developed turbulent spectrum of Alfvénic fluctuations that actively contribute to the acceleration of the solar wind through an extended wave pressure and dissipation at kinetic scales (Belcher 1971; Leer et al. 1982; Velli et al. 1991; Verdini & Velli 2007; Verdini et al. 2009). As a consequence of the extended energy deposition, the wind will reach higher speeds than a purely thermally driven outflow. Yet, the terminal velocity is fundamentally dependent on the processes occurring in the low corona, and a dense, slow wind can still emerge from a coronal hole if a strong expansion and/or heating is present below the sonic point (see, e.g., Velli 2010).

A second component, more intermittent and dynamic, originates from or nearby streamers, i.e., closed coronal structures. This characterization of the "slow" wind is backed by composition analyses, which show a strong enrichment in low first ionization potential (FIP) elements (see, e.g., Laming 2004). Some dynamical processes must then allow for exchange between the closed structures and the outflow (Antiochos et al. 2012; Higginson & Lynch 2018). Moreover, the wind surrounding the heliospheric current sheet (HCS) reveals density perturbations that are not observed in the fast wind. These density perturbations have been shown to be periodic, with high power concentrated around a 90 minute period, coherently over many different instruments at various places in the inner heliosphere (Viall et al. 2010; Viall & Vourlidas 2015). Sanchez-Diaz et al. (2017a, 2019) have further identified longer periods of 10–20 hr driven by the release of extended structures along with the aforementioned substructures of 10 times smaller periods. Recent solar wind measurements by the Parker Solar Probe (PSP), combined with remote-sensing observations taken near 1 au, provide additional support to the idea that the slow wind comes in at least two states with different bulk properties and levels of variability (Rouillard et al. 2020a).

The suggestion that reconnecting plasmoids at the tips of helmet streamers were a slow wind source had already been made by Sheeley et al. (1997) and Wang et al. (1998). The study of Sanchez-Diaz et al. (2017b), using observations of the Solar and Heliospheric Observatory (SOHO) and STEREO spacecraft, has recently shown that the formation of density structures in the corona were associated with inflows that support formation by magnetic reconnection. Endeve et al. (2003, 2004), while studying simple dipolar configurations, unveiled an instability that leads to magnetic reconnection and the periodic release of flux ropes from the tips of the streamers. These works identified thermal processes as the origin of the streamer's instability, the thermal conduction coefficient being the critical parameter to make the instability vanish. The periods obtained in these works are of the order of 15 hr, and thus can only account for the longer timescales observed in the inner heliosphere, which suggests that another mechanism might be at play.

The HCS is pinched at the tips of streamers, which can lead to a tearing instability (Furth et al. 1963; Biskamp 1986; Velli & Hood 1989; Loureiro et al. 2007; Pucci & Velli 2014; Tenerani et al. 2015). Einaudi et al. (1999) modeled the region above the cusp of a helmet streamer as a current sheet embedded in a broader wake flow. The combined tearing Kelvin–Helmholtz instability was shown to induce the acceleration of density-enhanced magnetic islands. The model was further developed by Rappazzo et al. (2005), who showed how the diamagnetic plasmoid expulsion in a spherical geometry would lead to a rapid plasmoid acceleration profile. However, these works started with a finite thickness current sheet and did not take into account either the specific geometry of the helmet streamer cusp, nor the natural thinning of the sheet arising from the converging plasma flow, although they did comment on the role such convergence might play in the evolution.

In this Letter, we revisit the instability of streamers using a resistive magnetohydrodynamic (MHD) model of the solar corona and wind. We investigate the effect of the Lundquist number S on the streamers' stability and the periodicity of the structures released in the slow solar wind. We show that, for sufficiently high S, a tearing instability is launched in the HCS and that its scaling properties with S are consistent with the ideal tearing scenario (Pucci & Velli 2014; Tenerani et al. 2015). We observe two main periodicities in the simulations: a small one related to the fastest growing tearing mode, and a longer one associated with cycles of linear onsets of the instability. Using the ideal tearing scalings to very high S typical of the solar corona, we show that these periodicities are fully compatible with the timescales of density structures observed in the slow solar wind.

We use the Alfvén-wave turbulence-driven MHD model presented in Réville et al. (2020), itself based on the PLUTO code (Mignone et al. 2007), and refer the reader to this paper for all the details on the numerical schemes and equations solved. For this study, we perform 2.5D simulations in spherical geometry of a dipolar configuration, corresponding to a solar minimum case. We focus on the properties of the HCS, located around θ = π/2. For the need of this particular study, some improvements have been brought to the numerics of the code. First, we use a HLLD Riemann solver (Miyoshi & Kusano 2005), which handles discontinuities and sharp gradients better than the previously used HLL solver (Einfeldt 1988). In particular, our studies of reconnecting current sheets using both solvers have shown that only the HLLD was able to correctly capture the growth of a tearing instability. Second, we use a constrained transport method (Balsara & Spicer 1999) to maintain ∇ · B = 0 at machine accuracy. Finally, explicit resistivity is introduced to ensure a good control of the Lundquist number close to the current sheet. The simulations are integrated on a non-uniform grid, strongly refined around the current sheet with Δθ = 10⁻³ for θ ∈ [π/2 − 0.1, π/2 + 0.1]. The radial grid is stretched, although we maintain an almost constant Δr = 10⁻² up to 20R_⊙, while the solution is extended up to 50R_⊙ with a coarser grid. The aspect ratio of the cells rΔθ/Δr is close to unity between 1 and 15R_⊙, which encompasses the growth region of the tearing instability that we will characterize in the following section. Although the resolution might seem coarse for this kind of problem, we use a fourth-order spatial scheme obtained by a parabolic reconstruction, which effectively ensures results that are consistent with linear theory up to S ∼ 10⁶ (see Section 3).

The model has already shown a very good agreement with bulk solar wind properties at the PSP first perihelion (see Réville et al. 2020). There are essentially four input parameters in the model: the amplitude of Alfvén waves at the coronal base δv_⊙ = 30 km s⁻¹, the base density ρ_⊙ = 5 × 10⁻¹⁶ g cm⁻³, the correlation length of the turbulence at the base of the corona λ_⊙ = 0.025R_⊙, and the input magnetic field. These input parameters are very close to the ones used in the study of the first PSP perihelion except for the input magnetic field, which was then given by a solar magnetogram. Instead, we use here a purely dipolar field of 5 G at the equator, which reproduces typical fast and slow wind bulk properties at solar minimum.

In Figure 1, we compare the wind speed profile obtained in our simulations with the speed of density structures observed in the inner heliosphere with the white-light instruments Solar and Heliospheric Observatory/Large Angle and Spectrometric Coronagraph (SOHO/LASCO) and Solar Terrestrial Relations Observatory/Coronagraph (STEREO/COR; Sheeley et al. 1997; Abbo et al. 2016). The fast wind speed in orange is obtained at high latitudes, while the slow wind speed profile is extracted at the HCS in our model. The slow wind speed is very close to the fit made to the observations in Rouillard et al. (2020b) and fully consistent with the speeds of the slow wind. The whole system of flows and density fields around the HCS is obtained self-consistently with the Alfvén-wave driven heating model and ensures that realistic physical conditions are met at HCS, where a tearing instability is triggered, as described in the following sections.

Zoom In Zoom Out Reset image size

This Letter presents six simulations in which the Lundquist number S = LV_A/η varies around the HCS. The value of L = 4R_⊙, the characteristic length, and v_A = v_kep ∼ 437 km s⁻¹ are the same in all simulations and corresponds to the growth region of the tearing instability and the Alfvén speed away from the current sheet. We only change η so that S varies between 10⁴ and 10⁶. At the beginning of the simulations, the magnetic field is a pure dipole and there is no current sheet. The first phase hence consists in the creation and the thinning of the HCS around θ = π/2 within a solar wind flow. In addition to the flow, latitudinal, equatorward Lorentz forces act to thin further the HCS (see, e.g., Réville & Brun 2017). The thinning process takes roughly 90 hr in all our simulations.

Studies of current sheet reconnection usually start with a given equilibrium of magnetic field and velocities. The Harris current sheet, for example, sets the field along the sheet, B_R in a cylindrical coordinate system (R, Z, Φ), to be

$$\begin{eqnarray}&&{B}_{R}(Z)=A\tanh \left(\displaystyle \frac{Z}{a}\right),\end{eqnarray} \tag{ 1 }$$

where a is the thickness of the current sheet. Recent works (Pucci & Velli 2014; Tenerani et al. 2015) have tried to assess what aspect ratios a/L could be physically constructed and lead to fast reconnection. In particular, Pucci & Velli (2014) showed that Sweet–Parker current sheets would generally be impossible to create because of the diverging (infinite) growth rates for the tearing instability as $S\to +\infty$. They concluded that current sheets would disrupt once a limiting inverse aspect ratio scaling as a/L ∝ S^−1/3 would be reached, at which point the tearing mode growth rate becomes independent of S and of the same order as that of ideal MHD dynamics.

In the top panel of Figure 2, we show the structure of the current sheet obtained in the simulations at the end of the thinning phase. The thickness of the current sheet is computed using Equation (1) and the value $| A| =| {B}_{Z}(R=5,Z=2)|$, the amplitude of the field away from the current sheet. As the Lundquist number increases the current sheet thins and we give the obtained scaling in the bottom panel of Figure 2. After the thinning process, different evolutions are observed depending on the value of S. For the lowest S = 1.2 × 10⁴, the current sheet is stable. For S = 2.4 × 10⁴, a regime of steady reconnection sets in and for all higher S, a tearing instability is triggered. The bottom panel shows the aspect ratio a/L and compares with theoretical scalings. The plain blue line is obtained with L = 4R_⊙, while the shaded area represents the interval L ∈ [2R_⊙, 6R_⊙]. Interestingly, the aspect ratios a/L that we get are always smaller than S^−1/3. The reason for this smaller ratio might be found in the presence of a small normal component of the magnetic field, whose stabilizing effect can lead to smaller inverse aspect ratios for fast ideal growth rates shown by Pucci et al. (2019). Nonetheless, at large S, the aspect ratio a/L approaches the S^−1/3 scaling, diverging clearly from the well-known Sweet–Parker scaling S^−1/2.

Zoom In Zoom Out Reset image size

Figure 3 shows the growth of the tearing instability for two unstable cases. The quantity in the background color is the latitudinal Alfvén speed v_A,θ, which is characteristic of the tearing mode. Using the Alfvén speed instead of the magnetic field component across the sheet is a way of accounting for the expansion of the wind. Analysis of the toroidal electric field yields similar results. The snapshots are taken at the end of the linear growth phase of the fastest growing mode k_m in the simulations. As S increases, we can see that k_m increases (or that the characteristic wavelength decreases). The growth region is found for all simulations to be between 5 and 9R_⊙, which explains a posteriori our choice of L = 4R_⊙.

Zoom In Zoom Out Reset image size

According to the ideal tearing scenario (Pucci & Velli 2014; Tenerani et al. 2015), when the current sheet aspect ratio a/L ∼ S^−1/3, we get the following asymptotic scaling for the fastest growing mode:

$$\begin{eqnarray}&&{k}_{m}L\propto {S}^{1/6},\,\gamma {t}_{{\rm{A}}}=C,\end{eqnarray} \tag{ 2 }$$

with the constant C ≃ 0.62 for the Harris equilibrium profile and taking precisely a/L = S^−1/3. Figure 4 shows the computed scaling of k_m for the fastest growing mode and the associated rate for all unstable cases. As discussed earlier, at our values of S the fastest growing mode has a growth rate that increases with S, and although our simulations only reach S ∼ 10⁶, the slope of k_m (S) is seen to decrease toward the predicted limiting scaling, illustrated by the black dashed line. The growth rate shown here also tends to converge toward a fixed value γt_A ∼ 2.5, or at least does not explode with increasing S. This value, higher than expected for Harris equilibria, could be the result of gradients in velocity and density, which are subsequent to the realistic wind conditions met around the current sheet in the simulations.

Zoom In Zoom Out Reset image size

Magnetic islands, or flux ropes in three dimensions, created by the tearing instability described in the previous section, do create a local enhancement of the density and their propagation in the solar wind could lead to the in situ and white-light observations of periodic density perturbations. In Figure 5, we show the propagation of the structures created by the tearing instability and the reconnection at the tip of the helmet streamer. The transverse Alfvén speed v_A,θ and the expansion compensated density nr² are taken at the HCS and stacked over time. The substructures of alternate signs in the Alfvén speed are the flux ropes created by the tearing mode. Density enhancements are associated with the edges of these flux ropes and reach about two to three times the local density, which is consistent with observations (see, e.g., Sanchez-Diaz et al. 2017a). The scales predicted by the tearing mode analysis of Section 3 can be computed as follows: because the simulation with the largest Lundquist number is close to the asymptotic regime, we can scale up the fastest growing mode flux ropes at S ∼ 10¹², which is the expected value in the solar corona:

$$\begin{eqnarray}&&{\lambda }_{m}(S={10}^{12})=1/{k}_{m}(S={10}^{6})\times {\left({10}^{-6}\right)}^{1/6}\sim 0.15{R}_{\odot }.\end{eqnarray} \tag{ 3 }$$

This length scale does not account for the expansion and acceleration of the wind. Figure 5 shows the characteristic curve of the propagation of a passive scalar advected by the wind speed near the HCS. This curve matches closely the kinematics of the islands and periodic density perturbations shown in the colorplots. When released, likely by a combination of flow dragging and melon-seed forces (see Rappazzo et al. 2005), density structures are advected by the slow solar wind. They gradually accelerate out to a heliocentric distance of 10–15 solar radii and then propagate at nearly constant speed. Moreover, we can see in Figure 5 that the structures grow as they propagate in the slow solar wind. This can be understood as a consequence of the solar wind acceleration in the growth region of the tearing mode. If these structures are pressure balanced and advected by the wind, two points initially separated by dλ = (x₁ − x₀) will eventually be separated by

$$\begin{eqnarray}&&d\lambda ^{\prime} (t)=d\lambda +{\int }_{t}v({x}_{1}(t))-v({x}_{0}(t)){dt},\end{eqnarray} \tag{ 4 }$$

after a time t (see also Bulanov et al. 1978; Shi et al. 2018). We find using the slow wind speed profile used in Figure 1 that between 3R_⊙ and 20R_⊙ a structure of characteristic scale λ_m = 0.15R_⊙ can expand by a factor 50 along the sheet. A straightforward integration of Equation (4) for x₀(0) = 4R_⊙ and dλ(0) = λ_m = 0.15R_⊙ gives dλ = 2R_⊙, out of the solar wind acceleration region. Provided that these structures propagate at some 300 km s⁻¹, we get a typical period of 80 minutes out of the solar wind acceleration region, which is very close to the periodicity of density structures reported in Viall et al. (2010), Viall & Vourlidas (2015). Depending on the precise birth location of the flux rope, however, the growth will vary, as the expansion is an increasing function of the time spent in the acceleration region.

Zoom In Zoom Out Reset image size

Finally, we also observe in Figure 5 a second periodicity that could be the longer one reported by the observations. The bottom panel of Figure 5 displays three episodes of linear onsets followed—for the first two—by a nonlinear phase of the tearing mode, that are periodically evacuated by the solar wind outflow. During the nonlinear phase, flux and plasma is transferred from the streamer to the HCS, as shown in the supplementary movie in the Appendix. The composition of the density structures thus should be close to the composition of the helmet streamer, which is consistent with observations. These episodes are separated by roughly 30 hr, which is slightly more than observations at activity maximum, yet in the right order of magnitude (Sanchez-Diaz et al. 2017a, 2019). This periodicity is of course controlled by the strength and duration of the nonlinear stage. We can expect that as S increases, and the size of the structures decreases, the solar wind flow will be able to evacuate more easily the perturbations, hence reducing the time interval between linear onsets. However, further studies are necessary to confirm this behavior of the nonlinear regime at high S.

In this work, we have shown that a tearing instability was naturally occurring at finite resistivity in simulations of the HCS with realistic outflow velocity and density gradients. At low Lundquist numbers, the current sheet is stable and no reconnection is observed. Increasing S, the current sheet thins, and reconnection is triggered via a tearing mode. The current sheet before the onset of the instability is thinner than the ideal tearing scenario for low S but the inverse aspect ratio approaches the scaling a/L ∝ S^−1/3 as S increases. The fastest growing mode and associated growth rate also reach an asymptotic regime, which is consistent with the ideal tearing scenario, with k_m ∝ S^1/6 and γt_A ∼ 2.5. Gradients and inhomogeneities in the simulations could explain the difference with the theoretical values predicted by Pucci & Velli (2014). Another reason might be the previously mentioned finite normal component in the neighborhood of the helmet streamer cusp. This has a stabilizing influence, requiring sheets to become thinner for tearing modes to occur on fast ideal timescales (Pucci et al. 2019).

The convergence toward an ideal regime allows some extrapolation at much larger S, characteristic of the conditions in the solar corona. We can expect that such a tearing mode for the actual HCS could give typical length-scales of about 0.15R_⊙ at the tip of the streamers. These flux ropes are then ejected and grow as pressure balanced structures inside the solar wind acceleration region. Accounting for this expansion, the length and timescales of flux ropes as well as the density structures that separate them are fully consistent with periods of a few hours observed in heliospheric images (Viall et al. 2010; Viall & Vourlidas 2015) and in situ data (Kepko et al. 2016; di Matteo et al. 2019).

This Letter is mostly focused on the linear regime of the tearing instability. Our simulations show an interesting behavior, where several linear growth phases of the tearing mode are evacuated with the solar wind. Such a recurring process could explain why the linear fastest growing mode has a peculiar status in the observations, as well as the longer periods of flux ropes events of 10–20 hr. However, the precise duration of these cycles is controlled by the properties of the nonlinear regime, and needs to be studied in further detail. Nonetheless, our results suggest that in the nonlinear phase plasma initially trapped on coronal loops below the helmet streamer and enriched in low FIP elements could be released in the slow solar wind (see the supplementary video in the Appendix).

Finally, it is important to stress here that our results differ quite fundamentally from the streamer's instability mechanism proposed by Endeve et al. (2003, 2004). Interestingly, while we do get siphon flows (see, e.g., Cargill & Priest 1980) for resolutions equivalent to that of these latter works, we do not observe any thermal instabilities at the relatively high resolution used for our present dipolar study. Hence, although thermodynamic processes related to coronal heating are likely to be responsible for reconnection in the low corona, they may not be necessary to explain the observations of periodic density perturbations. Moreover, a tearing-induced reconnection is likely much less dependent on the scale and properties of the streamers, which is consistent with regular periodicities of solar wind density structures over the solar cycle (see Viall et al. 2009). Analysis of early FIELDS and SWEAP data during the first PSP orbit seem to suggest that the whole heliospheric plasma sheet may be generated by the ensemble of reconnection exhausts originating from the tips of helmet streamers (Lavraud et al. 2020). In the near future the WISPR imager on board PSP may be close enough to the streamer cusp region to provide global observations of the reconnection process.

This research was funded by the ERC SLOW_ SOURCE project (SLOW_ SOURCE—DLV-819189) and the NASA Parker Solar Probe Observatory Scientist grant NNX15AF34G. V.R. thanks Fulvia Pucci for useful discussions. The authors are grateful to A. Mignone and the PLUTO development team. Simulations were performed on the Extreme Science and Engineering Discovery Environment (XSEDE; Towns et al. 2014) SDSC based resource Comet with allocation number AST180027, and GENCI supercomputers (grant 20410133). XSEDE is supported by National Science Foundation grant number ACI-1548562. This study has made use of the NASA Astrophysics Data System.

Figure A1 shows the onset of the tearing instability for the S = 2.4 × 10⁵ simulation. The time used for the still frame is t = 90.7 hr and corresponds to the first onset of the instability. The animation corresponding to Figure A1 displays the evolution of the density perturbations following the tearing mode and the multiple onsets of the instability.