© The Authors 2026

1. Introduction

Solar eruptions are extreme manifestations of the Sun’s magnetic activity, explosively releasing parts of the vast amount of energy contained within the solar interior out into the heliosphere, including to Earth (Schrijver et al. 2012). Strong eruptions, exemplified by the 1859 Carrington event, can severely affect technology and life on Earth and in its orbit, making space weather an important field of study for the sake of mankind. The most powerful events display bright flares, storms of solar energetic particles, and coronal mass ejections (CMEs). They can erupt from both active regions and long filament or prominence channels, due to a sudden destabilization of force-free magnetic fields in the corona.

To understand and eventually predict solar eruptions, Patsourakos et al. (2020) stress the need for a “clear understanding of the nature of the preeruptive magnetic field configurations of [CMEs]”. However, the mechanisms that bring the magnetic fields and plasma to such a preeruptive state are not fully understood, in spite of a plethora of theories (Aulanier 2014; Patsourakos et al. 2020). Many of these theories include a flux rope in the preeruptive magnetic field, but even that is contended by the magnetic breakout model (Antiochos et al. 1999). Flux ropes are also a common magnetic geometry in filaments that do not erupt (Mackay et al. 2010). In the present study, the focus is on analyzing in detail the mechanism in which flux cancellation leads untwisted arcades to form a twisted coronal flux rope. In doing so, we aim to increase the understanding of preeruptive magnetic fields, which in turn can lead to better prediction of eruptions.

There are several competing theories for the onset of eruptions related to the destabilization of prominence carrying magnetic flux ropes. Furthermore, there are multiple theories for the generation of these flux ropes. These theories are difficult to prove or disprove through observations alone, since flux ropes and their formation are not directly observable. They are only indirectly observable, through, for example, the corresponding cold and denser prominences or through field extrapolations from photospheric magnetograms (Nakagawa & Raadu 1972; Wiegelmann 2008), and through coronal cavities (Tandberg-Hanssen 1995; Gibson 2018). Therefore, it is important to supplement theories and observations with numerical simulations.

The suggested mechanisms for flux rope formation can in short be split in two: (i) Pre-twisted flux ropes emerging from the convection zone through the photosphere into the solar atmosphere, as was observed by Tanaka (1991), Lites et al. (1995), Leka et al. (1996); (ii) Flux ropes forming above the photosphere, due to flux cancellation in the photosphere, driven by plasma convection of the field line footpoints, motivated by observations (Harvey & Martin 1973; Martin & Harvey 1979; Wang et al. 1989; Livi et al. 1989; Green et al. 2011). Both theories have been studied extensively through numerical simulations.

The emergence of pre-twisted flux ropes has been simulated with various simplifications. Some simulations model the photosphere as the upper boundary, studying the rise of a pre-twisted buoyant flux rope inside an initially stratified convection zone, first in 1D and then in 2D (Emonet & Moreno-Insertis 1998; Jouve et al. 2013). Then kinematic emergence simulations brought the flux rope into the corona by modeling the rise of the flux rope through the convection zone kinematically, whereupon the physics in the atmosphere was described by a full magnetohydrodynamic (MHD) model (Fan & Gibson 2004; Amari et al. 2004, 2005). Also, the convection zone was modeled by the full MHD equations, both in 2D (Archontis et al. 2007) and 3D (Fan 2001; Archontis & Hood 2008), although the convection zone was initially stratified in these simulations. The flux ropes do not fully emerge in these simulations, but are tied to the photosphere due to the heavy plasma that they do not manage to carry into the atmosphere. Additionally, as the convection zone is initially stratified but not convective, they do not fully model the mechanisms of the Sun. Hotta & Iijima (2020) have recently run simulations of an emerging flux rope in a convectively unstable and relaxed convection zone, showing that it can be done, achieving the goal of forming sunspots.

How flux ropes can form in the atmosphere by flux cancellation was theorized in the seminal paper by van Ballegooijen & Martens (1989). A continued shearing of overlying coronal arcades, combined with a priori convectively driven converging flows toward a straight polarity inversion line (PIL), would cause reconnection as the opposing fluxes would cancel, gradually creating a twisted flux rope. This mechanism has typically been studied numerically in 3D by considering the photosphere as the lower boundary, where the magnetic field line footpoints are line-tied. The flux rope formation and dynamics is then driven by prescribed convection and/or diffusion of the field line footpoints at the boundary, in the manner described by van Ballegooijen & Martens (1989). This approach has been successfully repeated several times over the last two decades, using smooth surface flux distributions and flux cancellation all along well-defined PILs (Amari et al. 2003a,b; Mackay & van Ballegooijen 2006a,b; Aulanier et al. 2010; Jiang et al. 2024; Xing et al. 2024a,b).

A more critical problem of the line-tied flux-cancellation models, where it breaks physics, is that the magnetic flux diffuses as a scalar at the photospheric boundary. That leads to a diminishing flux until it disappears by diffusion (Aulanier et al. 2010, Fig. 5). In reality, these field lines reconnect and continue down below the photosphere where they are redistributed and reenter the magnetic structure. By not properly modeling the reconnection of field lines, which is at the root of the mechanism for forming flux ropes, there are valid concerns whether the flux ropes generated are physical or simply artifacts of the model and the prescribed driving at the boundary. Therefore, we need a proper treatment of flux cancellation to assess whether it can produce flux ropes in the corona or not.

While observations alone cannot explain all details of how flux ropes are formed, they do show several attributes of the real process that the models should replicate. Real flux cancellation does not occur along a straight PIL, as in several of the simulations described above. In observations, the flux cancellation is fragmented between patchy flux concentrations along the PIL, even in δ spots, between opposite polarity sunspots isolated from other flux, and in decaying active regions (Schrijver et al. 2011; van Driel-Gesztelyi et al. 2003; Wang & Muglach 2013). Flux ropes also typically do not form in a single cancellation event, but are rather built up through several events that contribute to the same rope or filament (Wang & Muglach 2007, 2013). Both of these are also true for flux ropes that erupt. In the observation by Green et al. (2011), a bipole with an initially smooth PIL connected by slightly sheared arcades turned into an eruptive sigmoid through several cancellation episodes between increasingly patchy polarities.

In this study, we simulate flux rope formation through flux cancellation of sheared coronal arcades, driven by stochastic convection from a convectively unstable convection zone. This counters many of the artificial approximations of previous flux-cancellation simulations. We do so by adding a linear force-free field (LFFF) in the form of sheared coronal arcades to a relaxed quiet Sun (QS) simulation in Bifrost (Gudiksen et al. 2011), which combines into an enhanced network (EN) with a patchy PIL. Then, we let the stochastic convection drive the evolution of the atmospheric magnetic fields. We find that the result is a stable twisted flux rope formed in the solar atmosphere, but that unlike in the smooth flux cancellation model, the flux rope formation results from a variety of local events, including previously undocumented convectively driven thick-photosphere tether-cutting (TPTC) reconnection.

2. Numerical method

2.1. The Bifrost code

The simulations discussed in this work were run with the parallel radiative-MHD code Bifrost (Gudiksen et al. 2011). The code is designed to study a box in the Sun, starting in the upper convection zone and extending through the photosphere into the corona. We refer to these vertical boundaries as the lower and coronal boundaries, respectively. The complex dynamics within this box in the Sun are modeled as realistically as currently possible by modeling several different physical processes in a modular approach. The implemented physics used here include radiative transfer, Spitzer thermal conductivity along magnetic field lines using the wave method discussed by Rempel (2016), hydrogen kept in local thermal equilibrium (LTE), and numerical hyper-diffusivity that recently was proven to model the resistivity realistically, even at poor resolutions (Færder et al. 2023, 2024).

There are various vertical boundary conditions (BCs) implemented in Bifrost, including constant, symmetric, antisymmetric, hydrostatic, and divergence-free magnetic field extrapolations. These again can be treated either as standard time-derivative equations or by the method of characteristics. The divergence-free magnetic field extrapolations suppress numerical growth of nonphysical ∇⋅B, which reduces the need for iteratively canceling it to suppress nonphysical magnetic monopoles. The horizontal BCs are typically periodic in Bifrost simulations. Other types of horizontal BCs are possible in Bifrost, but are typically not used. Open BCs lead to loss of convective power and hard wall BCs lead to strong constraints on the magnetic field and plasma motion close to the boundary. To get a compact simulation box with minimal artifacts, periodic horizontal BCs are preferred.

The staggered grid used here was 3D. It is uniform in the horizontal (x, y) plane and nonuniform in the vertical z direction, in such a way that it is better resolved closer to the photosphere and chromosphere than high up in the corona. In Bifrost, z increases along the line of sight (LOS), meaning it is positive in the convection zone and negative in the corona. The spatial derivatives are sixth-order and the time stepping is done with the explicit third-order Hyman predictor-corrector scheme (Hyman 1979). The state of the simulation was written to file every 10 s (solar time), hereafter referred to as “snapshots”. That includes for each grid cell the density, three momentum components, internal energy, and three magnetic field components: eight variables in total. One can additionally store other auxiliary variables, such as the temperature, that alternatively can be calculated a posteriori from the eight independent variables.

2.2. Initial conditions: A modified QS simulation

The simulation presented here started at a copied state from a QS simulation that had been relaxed and run for 8000 s = 133m20s, twice as long as necessary to typically reach a QS state in a Bifrost simulation. That the simulation had indeed reached a QS state, was determined by considering the average signed and unsigned flux at the photosphere, reported below, and checking for the absence of organization of the flux into poles, similar to the state in Fig. 3a. The solar times referred to in this paper are relative to this state, i.e., the copy was made and the solar time was reset to 0 s. This original QS simulation has a grid consisting of 768³ grid cells. The horizontal cross-section of the simulation was square with side lengths L_x = L_y = 24 Mm, divided uniformly such that the horizontal resolution was 31.25 km. The vertical direction extended 2.5 Mm into the convection zone, referred to as the lower boundary in this article, and –14.3 Mm into the corona, referred to as the coronal boundary. The vertical grid was not split uniformly, the resolution was higher close to the photosphere than in the corona.

Up until the copy was made, the QS simulation was subject to magnetic flux feeding from the convection zone when the plasma velocity through the lower boundary was directed into the simulation box. This was done to replace the magnetic flux that was convected out and to compensate for the lack of small scale dynamo. The injected flux was horizontal, first with a strength of B_y = 45 G for 4000 s, followed by B_y = 110 G for another 4000 s, until the time of the copy. As a result, at the time of the copy of the QS simulation, the middle signed horizontal flux was measured to be ⟨B_y(z = 2.5 Mm)⟩ ≈ 100 G at the lower boundary and about 10 G at the photosphere. This was due to the flux feeding. The vertical field was balanced with a middle signed flux of ⟨B_z⟩≈0 G, as it should be for a QS simulation, and a middle unsigned flux at the photosphere of ⟨|B_z(z = 0)|⟩ ≈ 20 G.

The copied state of the QS simulation was modified before restarting it. The resolution of the original simulation was first halved in each spatial dimension to conserve numerical resources, using the backstaff¹ tool (Frogner & Gudiksen 2024). This tool was designed to “provide a fast, reliable, and flexible framework for computations on Bifrost simulation data”, such as interpolations to different grids, tracing of field lines, and synthesizing of optically thin spectral lines. Then, the coronal boundary was extended from –14.3 Mm to –28.0 Mm above the photosphere by a hydrostatic extrapolation and adding of extra vertical grid cells. This successfully alleviated various numerical problems close to the coronal boundary, further elaborated in Appendix A. The new staggered grid consisted of m_x × m_y × m_z = 384 × 384 × 472 cells. That gave a horizontal resolution Δx = Δy = 62.5 km, while the vertical resolution varied between a minimum Δz = 24 km in the chromosphere and transition region and a maximum Δz = 155 km in the corona. Since restarting it, there has been no magnetic flux fed through the vertical boundaries.

The horizontal BCs were periodic. The vertical boundary ghost zones were filled by extrapolating from inside the main domain. The magnetic field BCs were set to be symmetric around the boundary in the horizontal components, while the vertical component was used to make the field divergence free. The density was extrapolated via the hydrostatic condition in the coronal ghost zones and logarithmically extrapolated in the lower boundary ghost zones. The internal energy was constant in the coronal ghost zones and extrapolated in the lower boundary ghost zones. To keep the lower boundary at a certain energy level, the simulation continuously edged the average density and internal energy toward some constant values. The momentum was kept symmetric around the border in both the coronal and lower boundary ghost zones.

Some trial and error was necessary to make the best possible numerical model for this physical effect. More discussion of this experimentation is presented in Appendix A.

2.3. Linear force-free magnetic fields (LFFFs)

According to van Ballegooijen & Martens (1989), flux ropes can be generated by first shearing potential magnetic arcades along the PIL whereupon the shear is increased by convective flows toward the PIL from both polarities. Here, we inserted already sheared arcades in the form of an LFFF. Such magnetic fields, B, and corresponding currents, j, have zero Lorentz force, thereby satisfying

$$\begin{array}{cc}\hfill \mathbf{\nabla }\mathbf{·}\mathbf{B}& =0,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill {\mu }_0\mathbf{j}& =\mathbf{\nabla }\mathbf{\times }\mathbf{B}=\alpha \mathbf{B},\hfill \end{array}$$(2)

where α is a shearing scalar, being constant for an LFFF, and μ₀ is the magnetic vacuum permeability. By taking the curl of Eq. (2), and using Eq. (1), one gets the Helmholtz differential equation

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\mathbf{B}=-{\alpha }^2\mathbf{B}.\end{array}$$(3)

In Cartesian coordinates, the solution of Eq. (3) in a box can be expressed as periodic harmonics in the x and y directions that are exponentially decaying along z into the atmosphere, (here toward negative z) as expressed by Nakagawa & Raadu (1972), Demoulin et al. (1989). When removing the dependence on the y coordinate, the magnetic field components become

$$\begin{array}{cc}\hfill B_x& =-{\displaystyle \sum _{n=1}^{n_h}}\frac{B_n{l}_n}{k_{\mathit{nx}}}cos(k_{\mathit{nx}}x)exp[l_n(z-z_{\mathrm{CZ}})],\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill B_y& =+{\displaystyle \sum _{n=1}^{n_h}}\frac{B_n\alpha }{k_{\mathit{nx}}}cos(k_{\mathit{nx}}x)exp[l_n(z-z_{\mathrm{CZ}})],\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill B_z& =-{\displaystyle \sum _{n=1}^{n_h}}{B}_{n}sin(k_{\mathit{nx}}x)exp[l_n(z-z_{\mathrm{CZ}})],\hfill \end{array}$$(6)

where z_CZ is a constant reference depth, here set to the depth of the lower boundary in the convection zone, B_n is the amplitude of the nth harmonic, and

$$\begin{array}{cc}\hfill k_{\mathit{nx}}& =\frac{2\pi }{L_x}n,\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill l_n& =\sqrt{k_{\mathit{nx}}^2-{\alpha }^2}=\frac{2\pi }{L_x}\sqrt{n^2-{\alpha }_{\mathrm{norm}}^2},\hfill \end{array}$$(8)

are the corresponding wave number and exponential decay rate, respectively. Furthermore, α_norm is the normalized shearing parameter, which has a maximum real value of 1. It has been assumed that the horizontal interval is x ∈ [0, L_x]. In the limit of n_h → 1 or z ≪ z_CZ (far above the lower boundary) the shearing angle between B and the x axis at the PIL is given by

$$\begin{array}{c}\hfill {\displaystyle \underset{n_h\to 1}{lim}}tan(\theta )=\frac{B_y}{B_x}=\frac{\alpha }{l_1}=\frac{1}{\sqrt{1-{\alpha }_{\mathrm{norm}}^2}},\end{array}$$(9)

such that, for example, α_norm² = 0.5 corresponds to θ = 45deg, which just happens to be the shearing constant used in the simulation presented in this paper. With additional harmonics (n_h > 1) the shearing angle is lower closer to z_CZ. That is caused by the higher harmonics, for which l_n/α is higher, before their field strength fall off quicker with height due to the factor exp[l_n(z − z_CZ)] in Eqs. (4)–(6), as discussed by Demoulin et al. (1989).

Due to the orthogonality of the different harmonics, one can generate different LFFFs by appropriately selecting the amplitudes B_n. The goal here is to make it dipolar and focused around the central PIL at x_PIL = L_x/2 as opposed to around the PIL at the x = 0 boundary. To achieve the second goal, one can set two requirements. First, assuming that there are more than one harmonic, require the horizontal derivative of B_z of ever higher orders to be 0 at z = z_CZ and x = 0, making B_z as small as possible in the close vicinity,

$$\begin{array}{c}\hfill {\frac{{\partial }^p{B}_z}{\partial x^p}|}_{x=0}\propto \{\begin{array}{cc}{\displaystyle \sum _n^{n_h}}{n}^p{B}_n=0\hfill & ,p\text{odd}\text{and}{n}_h>1\hfill \\ sin(0)=0\hfill & ,p\text{even}.\hfill \end{array}\end{array}$$(10)

Secondly, make the field as strongly dipolar as possible at x = L_x/2, by making the horizontal derivatives of the different vertical harmonics have the same sign,

$$\begin{array}{c}\hfill {\frac{\partial B_{\mathit{zn}}}{\partial x}|}_{x=\frac{L_x}{2}}=B_{n}cos(\pi n)=B_n{(-1)}^n>0,\end{array}$$(11)

which requires every other amplitude to be of opposite sign. It turned out that, in all the cases we tested, this was achieved automatically by solving Eq. (10).

To find n_h magnetic amplitudes, one requires as many linearly independent equations. The first degree of freedom is the overall amplitude of B, which can be set by B₁ = 1. All the amplitudes can later be scaled by a constant to get a different overall amplitude. To get the remaining n_h − 1 equations, one must solve Eq. (10) up to order p = 2(n_h − 1)−1; for example, for n_h = 3, one must go to the derivative of order p = 3 to get two additional equations. Coefficients for the first numbers of harmonics are given in Table 1. Except for B₁, the other B_n depend on n_h, and are of opposite sign, as required by Eq. (11).

The plan is to add the LFFF to an existing QS Bifrost simulation. In general, one introduces new Lorentz forces by adding a magnetic field to another, even if the added field is force-free. To illustrate this, imagine a configuration with a background magnetic field, B₀, with a force density, f₀. When you add an LFFF denoted B₁, this can introduce new force terms since

$$\begin{array}{cc}\hfill \mathbf{f}& =\mathbf{j}\times \mathbf{B}\hfill \\ \hfill & =({\mathbf{j}}_0+{\mathbf{j}}_1)\times ({\mathbf{B}}_0+{\mathbf{B}}_1)\hfill \\ \hfill & ={\displaystyle \underset{{\mathbf{f}}_{\mathbf{0}}}{\underbrace{{\mathbf{j}}_{\mathbf{0}}\times {\mathbf{B}}_{\mathbf{0}}}}}+{\displaystyle \underset{\text{new force terms,}{\mathbf{f}}_{01}+{\mathbf{f}}_{10}}{\underbrace{({\mathbf{j}}_0\times {\mathbf{B}}_1+{\mathbf{j}}_1\times {\mathbf{B}}_0)}}}+{\displaystyle \underset{{\mathbf{f}}_{\mathbf{1}}=\mathbf{0}}{\underbrace{{\mathbf{j}}_{\mathbf{1}}\times {\mathbf{B}}_{\mathbf{1}}}}}.\hfill \end{array}$$(12)

Note that the additional force terms, f₀₁ + f₁₀, are proportional to the field strength of the added LFFF. To minimize the consequential adjustment, the LFFF should be added incrementally. We detail how we did this in our experiment in Sect. 3.2.

2.4. Tracking field lines with corks in Vapor

In addition to the physics modules that make the simulations of Bifrost realistic, there are also modules implemented to aid the analysis. One such module of great relevance to this study is the corks module (Zacharias et al. 2018; Druett et al. 2022). This module was employed in a similar fashion by Robinson et al. (2022, 2023). In short, corks are Lagrangian point markers that are convected by the velocity-field in the simulation without affecting the physics. These markers track the plasma motion at the numerical time resolution of the simulation, not at the often much larger time resolution between output snapshots. Hence, the corks give an accurate description of the simulation evolution, otherwise lost by limiting the output cadence, due to disk storage concerns, to once every 10 s of solar time in this case.

The corks were inserted at a specific time, one cork per grid cell. To avoid over- or under-representation of any areas, the corks where pruned (removed) if there were more than two per grid-cell and added if there were zero in a grid-cell. For each snapshot, an identification number and the three spatial coordinates were stored per cork. Simple arithmetic shows that N grid cells will approach 2N corks, making the storage requirement comparable to that of the eight physical parameters of the plasma itself.

To accurately visualize the magnetic field lines in 3D, we used the Visualization and Analysis Platform for Ocean, Atmosphere, and Solar Researchers (Vapor) (Li et al. 2019). This software can visualize 3D parameters in a multitude of ways, such as showing the magnetic field lines integrated in 3D on top of a 2D magnetogram illustrating B_z at z = 0. In particular interest for this study, Vapor can select random seed points from which to start the integration of magnetic field lines in 3D, biased by another parameter. We often selected such seed points at a specific time by biasing the algorithm to find large values of

$$\begin{array}{c}\hfill \frac{|\mathbf{j}|}{{\mu }_0|\mathbf{B}|}=\frac{|\mathbf{\nabla }\mathbf{\times }\mathbf{B}|}{|\mathbf{B}|}·\end{array}$$(13)

Equation (13) is a normalized measure of the current density or complexity of the magnetic field, almost identical to the absolute value of the shearing parameter α in Eq. (2). We selected corks close to these seed points and found where those corks had been or would be at different times to track and analyze the evolution of the magnetic field lines and connectivity in particular.

Tracking individual field lines is a key component of the analysis presented in Sect. 4. While the corks module proved utterly necessary in this endeavor, it has its limitations, some of which are further elaborated in Appendix B.

3. Numerical experiment

The Bifrost simulation presented and analyzed in this paper was run as following: (i) Relax a modified QS simulation for 1000 snapshots (166m40s), starting at a reference time of t₀ = 0 s; (ii) Ramp up the LFFF for 120 snapshots until solar time 186m40s; (iii) Initialize the corks at solar time 203m30s; (iv) End the simulation at solar time 368m00s.

3.1. Relaxing the modified QS simulation

We started the simulation as a relaxed QS simulation so as to not have a smooth idealized photosphere with a straight PIL, but rather one that was more comparable to real-Sun observations. After the resolution was reduced and the coronal boundary was extended, as is described in Sect. 2.2, it was necessary to first relax the modified initial condition. As is discussed in Appendix A, it was found more than sufficient to relax it for 10 000 s = 166m40s. In Fig. 1, we display the horizontal averages of the temperature, density, plasma pressure, and magnetic pressure in snapshot 1000 after this relaxation.

Fig. 1. Horizontal averages of the Bifrost simulation. The magnetic pressures (PB) are given in dyn/cm2 for the QS simulation before the ramp at snapshot 1000 (dashed green), for the time-independent LFFF (dashed red), and for the simulation after the ramp at snapshot 1120 (dashed purple). The pressure (orange), the density (light blue), and the temperature (dark blue) of the atmosphere are plotted at snapshot 1000. The density is given in g/cm3, multiplied with a factor 1012 to fit on the same axis as the magnetic pressures, while the temperature is given on the right axis. The shaded gray area (z ∈ [0, −1] Mm) represents the “thick” photosphere, up to the altitude of the temperature minimum, i.e., the height where the density stratification is the strongest (see Sect. 3.4 for more details).

A vertical 2D cut of the magnetic field is shown in Fig. 2a and a horizontal 2D cut of the magnetogram at the photosphere is shown in Fig. 3a. After this relaxation without flux feeding through the bottom boundary, the vertical field was still balanced, while the mean signed flux in the y direction had dropped to below 1 G, both at the lower boundary and at the photosphere. Hence, the simulation was ready as a blank slate for being injected with an LFFF.

Fig. 2. 2D cross sections of the magnetic field at y = 10 Mm. The streamlines show field lines, but their density does not signify the field strength. The field strength is shown by the color. (a) is the QS field right after the relaxation, (b) is the LFFF, (c) is the sum of the two proceeding, while (d) shows the EN right after the 20 min long ramp of the LFFF into the QS simulation.

Fig. 3. Magnetograms at the photosphere (z = 0) in the Bifrost simulation. (a) is just after the relaxation of the QS. (b) is just after the ramp of the LFFF. The white polarity (north pole) signifies flux coming out of the solar surface, being negative here as z is positive along the LOS. The black polarity (south pole) signifies flux entering the solar surface. The associated movie is available online

3.2. Ramping of the LFFF

The LFFF added to the simulation was chosen to fulfill several requirements. These requirements represent our hypothesis for a field with a high likelihood of eventually forming a flux rope through several episodes of flux cancellation and magnetic reconnection as in van Ballegooijen & Martens (1989), but driven by the stochastic convection. The shape was to be dipolar in B_z at the lower boundary, sheared along the PIL and be stronger across the PIL in the middle of the simulation box than at the x boundary. The strength was to dominate the preexisting magnetic pressure and plasma pressure in the atmosphere, but not dominate the plasma pressure in the convection zone. Through a thorough analysis, we decided to define the shape of our LFFF by α_norm² = 0.5, to make the shearing angle θ ≤ 45°, and n_h = 30, to make the field sufficiently focused around the central PIL, as illustrated in Fig. 2b. After finding the relative strengths of the coefficients, the overall scaling of the field was decided by setting the strongest vertical field at the photosphere to max[B_z(z = 0)] = 50 G. This value was chosen to ensure that there would be coronal arcades from which the flux rope could be built after the insertion. Simultaneously, the arcades should not be present in the convective layer, because the convection zone was to self-consistently drive the motion. This is achieved, as can be seen in Fig. 1, where the magnetic pressure of the LFFF does dominate the QS magnetic and plasma pressure in most of the atmosphere, but not in the convection zone.

Immediately adding the designed LFFF to the QS simulation gave the magnetic field illustrated in Fig. 2c. The LFFF dominates the original magnetic field in most of the box, except for in the convection zone and highest up in the corona. To add the field this abruptly perturbed the QS simulation too much, which often crashed after significant plasma flows. The crashes typically occurred at the coronal boundary when calculating the Spitzer conductivity. Even without the crashes, the induced plasma flows were deemed an artifact of the sudden change of magnetic field and therefore undesirable in a realistic modeling of this mechanism.

To reduce the reactions to the injected field, the LFFF in the presented simulation was instead slowly introduced in 120 equal parts every 10 s, totaling up to 20 min of solar time. This ramp of the LFFF was conducted from snapshot 1000 to snapshot 1120, including a 10 s cool-down period after the last injection. The hyper-diffusivity was kept slightly larger during the ramp and an additional 1000 s, to diffuse small-scale effects that caused crashes after the injection. This was acceptable as we were interested in long-term large-scale effects. The magnetic pressure of this EN after the ramp is quite close to that of the LFFF in most of the atmosphere, as seen in Fig. 1. A vertical 2D cut of this magnetic field after the ramp is shown in Fig. 2d, which is appreciably different from Fig. 2c, particularly higher up in the corona.

The change in the simulation caused by adding the LFFF is clearly seen from the QS magnetogram in Fig. 3a to the EN magnetogram in Fig. 3b. After the ramp, there is a clear separation in white and black poles. The stochastic and self-consistent plasma convection has moved the field lines to the intergranular lanes, where they have concentrated and consequently increased the field strength locally. While the maximum field strength of the smooth LFFF was 50 G, these poles have much higher maximum field strengths. One thing lacking from this plot, combined with Fig. 2d, is that the field lines are also sheared. In Fig. 3b, the field lines would cross the PIL at approximately 45° upward to the right. This is exemplified by the field lines in Fig. 4a.

Fig. 4. Gradual formation and evolution of a flux rope in the Bifrost simulation. Each subfigure contains two panels, an angled view above a vertical top view of the same set of field lines. The angled views are identical in all subfigures, oriented such that the x axis increases down to the right and the y axis increases down to the left. The numbers on the axes are in units of Mm. (a) is right after the end of ramp of the LFFF, and differs from the rest by only showing field lines seeded at x = 12 Mm, y ∈ {5, 10, 15} Mm, and z ∈ (0, −5) Mm in gradually brighter colors. (b)–(f) show field lines representative for the flux rope in yellow-red colors, darker colors are shown for larger values of |J|/|B|, and overlying arcades in cyan. (d) also shows the four poles W1, W2, B1, and B2 referred to in the text.

3.3. Stochastic driving by self-consistent convection

During and after the ramp of the LFFF, the self-consistent Bifrost simulation controlled what happened to the magnetic field. In other words, there was no prescribed convection pattern, and the surface flows developed freely and self-consistently, following the laws of physics. Since the plasma beta β_p = P_p/P_B ≫ 1 in the convection zone, the plasma dictated the motion of the “photospheric footpoints” of the coronal parts of the magnetic field lines. During the ramp, the convection moved the magnetic flux into the granular lanes, as indicated by the end of the ramp in Fig. 3b. After the ramp, the self-consistent convection continued to stochastically make positive and negative polarities converge and cancel in the photosphere. This can be seen as the mutual cancellation of black and white polarities in Fig. 3 (see additional online made available online for long-term evolution of the magnetogram). According to van Ballegooijen & Martens (1989), these cancellations can cause reconnection and flux rope formation.

3.4. Reconnection in the thick photosphere

Before describing what happened in the simulation, we find it instrumental to stress that the photosphere is not an infinitesimal boundary or plane. We define the photosphere as the extended region located between z = 0 and the height of the temperature minimum, z(T_min), where the density stratification is the strongest with the smallest pressure scale-height. The width of this layer is typically approximated in the literature as 0.5 Mm (Priest 1984). However, the thickness of the photosphere is not a constant, it varies in space and time.

In the simulation presented in this article, the altitude of the temperature minimum is z(T_min) = 0.73 ± 0.24 Mm after ramping of the LFFF. There are even transient excursions up to 2 Mm. To be conservative, we define the average thick photosphere to extend up to 1 Mm above z = 0. We illustrate the width of this extended thick photosphere as the gray area in Fig. 1. It is within this layer that we have observed much of the atmospheric reconnection stochastically driven by the self-consistent convection.

4. Flux rope generation by flux cancellation

The primary goal of this work is to know whether or not a twisted flux rope could form in the atmosphere of a stochastic self-consistent realistic radiative MHD simulation. We actually find that it can.

4.1. High-level summary of the flux rope generation

The gradual generation and evolution of such a flux rope, over 2.5 h of solar time after the end of the ramp of the LFFF, is shown in Fig. 4. Fig. 4a shows the shearing angle of the LFFF being 45° high up in the atmosphere, but slightly less closer to the photosphere. The arcades are also less smooth closer to the photosphere, where the plasma pressure dominates the magnetic pressure. No appreciable twisting is visible in Fig. 4b, 30 min later. A twisted flux rope is gradually formed through Figs. 4c and 4d, with footpoints in a white polarity (W1) and a black polarity (B1). In Fig. 4d, there are also cyan field lines going from a white pole (W2) at (x, y) = (9.5 Mm, 15.5 Mm) to a black pole (B2) at (x, y) = (13.5 Mm, 9.5 Mm). The flux rope W1B1 later reconnects with these field lines W2B2 to form a longer flux rope in Fig. 4e, before the black pole B1 starts to disentangle toward Fig. 4f. The resulting flux rope in Fig. 4e has a maximum footpoint separation of approximately 12 Mm. At this time, the top of the flux rope extends vertically up to z = −2 Mm into the solar atmosphere.

The secondary goal of this work is to understand how such a flux rope is formed. In particular, how the forming and twisting of the flux rope are related to the flux cancellations at the photosphere. After the ramp of the LFFF, described in Sect. 3.2, the simulation was left to its own devices. Various behaviors were observed, of which we in the following analyze and discuss four distinct processes in the order in which they occur during the cancellation process: (1) Assumed slipping reconnection, aiding in forming the flux rope between Figs. 4b and 4c; (2) Flux cancellation linked to U-loop emergence, increasing the twist of the flux rope between Figs. 4c and 4d; (3) Flux cancellation linked to Ω-loop submergence (also between Figs. 4c and 4d), located right under the flux rope, but not increasing the twist of the flux rope appreciably; (4) Thick-photosphere tether-cutting (TPTC) reconnection and flux cancellation between a preexisting flux rope and comparably untwisted arcades between Figs. 4d and 4e.

4.2. Slipping reconnection

Possible examples of slipping reconnection are shown in Fig. 5. The colored field lines were selected at the time shown in Fig. 5e, as it was close to the end of the slipping process we were observing. That made it ideal for finding field lines that were or had been slipping, whereupon we studied them backward in time to see how they arrived at such a configuration. The chosen field lines had large values of the expression in Eq. (13), i.e., they were related to strong and narrow current layers, and their footpoints were in the canceling white polarity within the red circle. We chose three separate groups of field lines that passed through the cross section of the resulting flux rope in physically separated areas, to see how they wrapped around each other. The field lines were tracked with corks chosen between the middle of the field lines and their footpoints in the black polarity at this time.

Fig. 5. Slipping reconnection in the Bifrost simulation. Combined, the subfigures labeled (a)–(f) show the evolution until the time ta = 245m40s when this slipping process is complete. Each subfigure contains a top panel with an angled view of three families of 3D field lines (cyan, green, magenta), tracked by corks, above a 2D magnetogram. The angled view is angled from the bottom right corner of the 2D magnetogram toward the top left corner. On the bottom of each subfigure is a panel with only the 2D magnetogram, seen from directly above. The red circles in the bottom panels of each subfigure all highlight the same area of interest.

All the field lines remain rooted in the black pole (B1) at approximately (x, y) = (15 Mm, 6 Mm) throughout this process and beyond. The footpoints in the white polarity, however, are gradually slipping to the same white pole (W1) at (x, y) = (12.1 Mm, 10.6 Mm) within the red circle in Fig. 5f. First the green field lines slip, then the magenta, then the cyan. Each family spend 6–10 min on this slipping, in which the photospheric footpoints in the white polarity move approximately 4 Mm. Hence, the average field line footpoint slipping velocities are 7–11 km/s. The maximum slipping velocity observed in this Fig., which has 2 min cadence or more between the panels, is approximately 25 km/s. The peak velocities are higher. This matches slipping velocities observed during a solar flare by Dudík et al. (2014), which ranged from 4.5 km/s to 136 km/s. The slipping of the footpoints in Fig. 5 stops after the last shown snapshot. Later, these field lines remain twisted around each other, in approximate this position, as a flux rope that later becomes more twisted.

During the slipping reconnection, there is convergence and eventually cancellation between the black pole and W1 within the red circle in Fig. 5a. This convergence can be forcing the slipping reconnection through a quasi-separatrix layer that resides higher up in the corona. Since the reconnection does not occur in the photospheric layer, this process is different from the process proposed by van Ballegooijen & Martens (1989). Nevertheless, it still contributes to the generation and twisting of the flux rope.

4.3. U-loop emergence

An example of U-loop emergence is shown in Fig. 6. This is detectable as flux cancellation in the magnetogram. The green field lines were selected to go down into the intermediate black pole. They were tracked by corks selected along the longest visible part of the field lines pre-emergence, 50 s before t_b, close to the intermediate white pole. The pale yellow field lines represent approximately the core of the flux rope.

Fig. 6. U-loop emergence in the Bifrost simulation. From left to right, there are three different times relative to tb, a representative time for the emergence: (a) 100 s earlier, (b) tb, (c) 200 s later, after the end of the relevant flux cancellation. The upper panels of each subfigure show the 3D field lines of interest in green, tracked by corks, and pale yellow field lines representing the core of the flux rope, above a 2D magnetogram. The view in the upper panels is angled toward lower y. The lower panels show the same 2D magnetogram, from directly above. The cancellation of interest occurs within the red circle.

The entire area of interest within the red circle is convected to the left, toward the white polarity footpoints of the flux rope. The black polarity cancels mostly against the white polarity to the left of it in Fig. 6a. As the U loop starts to emerge in Fig. 6b, and the field lines to a lesser extent are weighed down by the denser plasma, the field lines curl around the flux rope in Fig. 6c, which increases the overall twist of the flux rope. That this U-loop emergence led to increased twist was dependent on the field line connectivity and orientation relative to the flux rope.

Judging from the magnetogram alone, this was a candidate for reconnection à la van Ballegooijen & Martens (1989). That could have been the case if the long green field lines ending in B1 instead would have been emerging from the white polarity to the left of the black pole in focus, instead of the white polarity to the right of it. It was necessary to visualize the 3D field lines and track them with corks to deduce what happened in this event. Even though it is not a reconnection event, this flux emergence is still a flux cancellation event that adds twist to the flux rope.

4.4. Ω-loop submergence independent of the flux rope

An example of Ω-loop submergence is shown in Fig. 7. This is detectable as flux cancellation in the magnetogram. The three sets of colored field lines were tracked by corks selected at time t_c, corresponding to Fig. 7b. The green and cyan field lines were selected to go down through different areas of their shared black pole. The corresponding corks were selected closer to their footpoints in this black pole, about 20–40% along the part of the field lines visible above the photosphere. The magenta field lines were selected to come up through the white pole that later cancels. The corresponding corks were selected close to the middle of the part of the field lines visible above the photosphere. A few of the cyan field lines follow the magenta lines into the Ω loop at t_c. Similarly, beyond the Ω loop (meaning at lower x), most of the magenta field lines, a few above and many below the photosphere, go in the same direction as the cyan field lines toward their footpoints in the white polarity.

Fig. 7. Ω-loop submergence in the Bifrost simulation. From left to right, there are three different times relative to tc, a representative time for the submergence: (a) 300 s earlier, (b) tc, (c) 250 s later, after the end of the flux cancellation. The top panel of each subfigure shows three families of 3D field lines of interest in different colors (cyan, green, magenta), tracked by corks, and pale yellow field lines representing the core of the flux rope, above a 2D magnetogram. The middle panel shows the same, but with a partially transparent photospheric magnetogram, to reveal the field lines below the surface. The bottom panel shows only the 2D magnetogram, where the cancellation of interest occurs within the red circle. We stress that the two upper rows are viewed from the right, toward lower x, to better visualize the dynamics of interest.

Most of the area of interest, within the red circle, is convected horizontally toward larger y (downward in the 2D magnetograms in the bottom row). The white pole within the red circle in Fig. 7b cancels evenly against the black pole at larger y and the one at larger x, before reaching the state in Fig. 7c. The magenta field lines are convected downward along z, into the convection zone, between these two snapshots. The previously photospheric Ω loop is submerged in Fig. 7c. This submerging and related flux cancellation is not linked to reconnection for most of the magenta field lines. This can be seen by how the magenta field lines are connected beyond the, now, submerged Ω loop (meaning at lower x). Most of the magenta field lines are still going in the same direction as in Fig. 7b, meaning along the cyan field lines that are visible above the photosphere toward their footpoints in the white polarity. A single magenta field line follows the green field lines above the photosphere. This reconnection happens below the photosphere, in the high plasma β regime. In other words, it is not the kind of coronal reconnection leading to a coronal flux rope that we are looking for. This cancellation does not contribute to the twist of the flux rope.

Judging from the magnetogram alone, this cancellation between the footpoints of the flux rope was also a candidate for reconnection à la van Ballegooijen & Martens (1989). However, as for the event in Sect. 4.3, the analysis of how the field lines evolved in 3D clearly revealed that it was not. Since the green and cyan field lines both enter the black pole in the same direction, they cannot reconnect to make a longer field line twisted around the pale yellow flux rope. We note not this event partially overlaps in time with the U-loop emergence discussed in Sect. 4.3, without them seeming to strongly affect each other.

The cancellation discussed in this section is an Ω-loop submergence, but it happens just after an Ω-loop emergence. The white pole is not visible in Fig. 7a, because the Ω loop is still below the photosphere. The loop is convected upward, through the photosphere, until t_c, before it is going down again. Such a peak-a-boo behavior of flux appearance and cancellation is an example showing that cancellation is not always related to reconnection or twisting of a flux rope. In comparison to other types of simulations described in the introduction, this mechanism could not have happened naturally on its own, meaning without prescription, without modeling the convection zone.

4.5. Thick-photosphere tether-cutting (TPTC)

An example of TPTC reconnection is shown in Fig. 8. It is detectable as flux cancellation in the magnetogram. There is no clear PIL between the white and black polarities, a trait it has in common with cancellation events often observed. That is also part of what separates this simulation from the other cancellation simulations that were mentioned in the introduction. Furthermore, the cancellation is more extended in both spatial and temporal dimensions than the events previously discussed, making an equivalent event on the Sun more easily detectable. The overall geometry of the field lines involved is similar to that in the tether-cutting model for solar eruptions (Moore et al. 2001). However, the altitude of the reconnection site is located neither in the corona as during a flare, nor in the transition region as in models of eruptive flux ropes formed by flux emergence (Manchester et al. 2004; Archontis & Hood 2008), but instead the reconnection occurs within the photosphere, as in the original flux-cancellation model by van Ballegooijen & Martens (1989).

Fig. 8. TPTC reconnection in the Bifrost simulation. In this figure, time advances downward, unlike in the previous figures, while all three panels per subfigure row show different views of the same snapshot. The left panels show four families of 3D field lines of interest in different colors (cyan, green, magenta, yellow), tracked by corks, above a 2D magnetogram. The magnetogram in the left panel in (c) is transparent by choice to show the magenta field line below the surface. The middle panels show the same, but from directly above. The right panels show only the 2D magnetogram, from directly above, where the cancellation of interest occurs within the red circle. There is no red circle in (d) as the cancellation then has completed.

All four families of colored field lines, displayed in Fig. 8, were chosen to have large values of the expression in Eq. (13) at t = 291m40s, corresponding to Fig. 8a. The yellow field lines correspond to the flux rope we have been presenting the gradual twisting of so far. Its footpoints are W1 and B1, and therefore we denote it W1B1. The footpoints are marked in Fig. 8a. The corks were selected between W1 and the middle of the field lines. As is mentioned in Sect. 4, here we show how this flux rope W1B1 reconnects with W2B2, the green untwisted field lines in Fig. 8a. The corresponding corks were selected closer to B2 than W2. The corks of the magenta field lines, W1B2, were selected in the middle between the two footpoints. The corks of the cyan field lines, W2B1, were selected approximately above B2.

The photospheric area of interest, within the gradually shrinking red circles, display significant flux cancellation over a much more extended time than the events previously discussed. The first and last snapshots displayed are separated by over 40 min. Within this time, W1 and B2 gradually approach each other through Figs. 8b and 8c, before having canceled in Fig. 8d. However, this non-line-tied simulation highlights that this flux cancellation is associated with a TPTC reconnection, including an X-crossing. The reconnection does not take place exactly in the z = 0 plane where we visualize the flux cancellation, but over an extended volume within the photosphere, mainly within z ∈ [0, −1] Mm, which we call the thick photosphere. There, the plasma β = P_p/P_B is still high, but decreasing, and the temperature is close to a minimum, as seen from Fig. 1.

Following the flux cancellation, W1B2, the magenta Ω loop, is submerged and transported by solar convection down and out through the bottom of the simulation box, which is at 2.5 Mm below the photosphere. In Fig. 8c, there is only one cork still within the simulation box; hence, the single magenta field line below the photosphere. As W1B2 (magenta) is submerged, W1B1 (yellow) and W2B2 (green) reconnect in such a way that the flux rope becomes longer, stretching from W2 to B1. In Fig. 8c, the green and yellow field lines still provide the generated flux rope with feet in W1 and B2. In Fig. 8d, after the complete cancellation of these poles, the feet are gone and the green and yellow field lines are twisted around each other, with a maximum photospheric footpoint separation of 12 Mm. The flux rope extends about 2 Mm into the solar atmosphere. The black pole B1, which have been quite steady so far, now starts to break up.

As for the previously discussed events, one could not have concluded from the magnetograms alone that this cancellation was reconnection à la van Ballegooijen & Martens (1989). The 3D field lines were necessary to analyze the snapshots individually, while the corks used to track the field lines over time were necessary to analyze the dynamics. This analysis confirmed the picture from van Ballegooijen & Martens (1989) for a realistic treatment of convection, including the occurrence of the submergence of the small reconnected Ω loop. Additionally, our analysis revealed that there is a small but nonzero difference in altitude within the photosphere between where the magnetic fields are measured and where the reconnection occurs (i.e., up to 1 Mm). Hence, the reconnection is driven by the converging motions of the polarities toward one another, not the cancellation itself. Therefore, the traditional picture of reconnection has to be expanded to allow reconnection to start before the cancellation between opposite flux concentrations.

5. Discussion

A key motivation for this study is to see if flux ropes actually can form above the photosphere from flux cancellation, as previously has been shown in line-tied simulations with prescribed driving at the photospheric boundary. Such prescribed driving has the benefit that one can control everything, but one has to know how the Sun works for that control to represent reality. In other words, all ropes that have come out of such simulations could have been artifacts of the prescribed driving. This is just one of the many logical arguments against these simulations. However, everyone does it – one of the authors has run many – but it does have its assumptions that have not been treated before. A key contribution of the present study is that it shows that flux ropes in fact can be generated also without these assumptions. Hence, the simplifying assumption of line tying is not a prerequisite for qualitatively forming flux ropes in the atmosphere. However, the line tying assumption still has other problems such as the decaying signed flux discussed in the introduction.

The present study took on an exploratory approach in simulating the formation of flux ropes through flux cancellation without the assumption of line-tied footpoints at the photosphere. This methodology for forming a flux rope has been proven to work with the following set of parameters for the LFFF: n_h = 30, α_norm² = 0.5, and max[B_z(z = 0)] = 50 G. Only one set of parameters has been tested because of the numerical cost and complexity of running such a simulation, as is expanded upon in Appendix A. The simulation did produce a flux rope as hypothesized, through multiple distinguishable flux cancellation events along a non-smooth PIL, driven by self-consistent convection.

We gradually ramped up the LFFF in the QS simulation. As a result, sheared arcades appear in the Bifrost simulation. That is similar to the configuration after steps (a) and (b) in the theory by van Ballegooijen & Martens (1989, see Fig. 1). The next step of the theory, and several other simulations discussed in the introduction, was to further shear the arcades by a prescribed convection of the plasma and field line footpoints toward or along the PIL, until the arcades would reconnect. In Bifrost, as in the Sun, the convection is not prescribed, but stochastic and self-consistently driven by the convection zone. This causes the stochastic patterns exemplified by the magnetograms in Fig. 3. That is the most significant difference of the idealized theory and simulations to the more realistic simulation presented in this paper.

The formation of the flux rope in the simulation bears some resemblance to observations referred to in the introduction. The flux rope forms along an irregular PIL that separates two areas of patchy polarities, as in Schrijver et al. (2011), Wang & Muglach (2013). Throughout the simulation, multiple flux cancellation events, separated in time and space, contribute to the formation of the flux rope, as in Wang & Muglach (2007, 2013). The resulting flux rope becomes gradually more aligned along the irregular PIL. However, the resulting flux rope is rather short and low, possibly due to the limited shear of the initial arcades or size of the simulation box. Hence, it does not warrant quantitative comparison to observed flux ropes.

Now that the process of forming flux ropes through cancellation of sheared coronal arcades has been established with this simulation, it would of interest to test other field configurations. To use an LFFF with a high number of harmonics (n_h) caused the field to be centered around x = L_x/2 as intended, but it also caused the shearing angle to be smaller close to the photosphere than higher up in the corona, as visible in Fig. 4a. This is opposite to the configuration in van Ballegooijen & Martens (1989), where there were several sheared arcades in parallel over an even more sheared field line. It is also opposite to the observations by Schmieder et al. (1996), which showed that the shear is concentrated at lower heights, close to the PIL. Similar setups to the one presented, with different parameters for the LFFF, are expected to lead to qualitatively similar results, albeit with quantitative differences; for example, a stronger shear would likely cause a longer flux rope. Hence, instead of inserting an LFFF, one could add a nonlinear force-free field with a shear distribution more in line with observations. Such a simulation would be more quantitatively comparable to observations.

The work presented in this paper focus on the cancellation process in a dispersed flux concentration environment, which is the case for decaying active regions. Flux cancellation can also happen between emerging sunspots. That was modeled in Rempel et al. (2023), in which they formed a flux rope and an eruption. This simulation was run with MURaM, which also can simulate atmospheric dynamics driven by stochastic and self-consistent flows in the convection zone. In order to build the full picture of flux rope formation in both situations, it may be interesting to also perform a similar kind of detailed analysis to the one presented here, on the reconnection process in the setting of colliding sunspots.

6. Conclusion

In this paper, we have shown and analyzed the stepwise formation and twisting of a flux rope in a self-consistent radiative MHD simulation including a free convection zone. Other simulation studies on flux rope formation often limit themselves to simulating the atmosphere, keeping the field lines line-tied to the photospheric boundary where there is prescribed convection or diffusion to drive the cancellation, or driving it by flux-emergence through a stratified convection zone and atmosphere with reconnection happening in the transition region, often leading to eruptions. A key aspect of the present study, in comparison, is the modeling of the stochastic convection of an unstable convection zone. This extension was necessary to include the complexities of the real Sun, as well as being able to understand what happens below the photosphere post-cancellation and post-submergence of Ω loops. The simulation was run with the state-of-the-art parallel code, Bifrost. To enable the tracking of field lines without following the prescribed convection of their footpoints at the boundary, we used Lagrangian markers, denoted corks, to accurately track the motion of selected field lines with time.

Within the simulation, which extends 28 Mm into the corona, an LFFF with sheared arcades was slowly introduced into a relaxed QS simulation. After the careful introduction of the LFFF, the magnetic field at the photosphere was mainly bipolar on the large scale, but patchy and without a straight PIL on the small scale. This is qualitatively the same as observed PILs on the Sun. Following this initial setup, the stochastic self-consistent plasma flows in the convection zone formed a flux rope above the photosphere, through several events connected to observable flux cancellation in the photosphere. Over the course of about 2.5 h of solar time, a flux rope was formed with a maximum horizontal footpoint separation of approximately 12 Mm and extending 2 Mm into the solar atmosphere (Fig. 8d).

The gradual formation of this flux rope included multiple events, separated in time and space. A few representative events with different atmospheric responses have been studied in more detail. A slipping reconnection event forms the initial flux rope. This event might be driven by observed flux cancellation at the z = 0 plane, but the reconnection occurs higher up in the atmosphere. A U-loop emergence later contributes to the twist of this flux rope. However, the flux cancellation observed as the loop emerges is not linked to reconnection. An Ω-loop submergence occurs shortly after, between the photospheric footpoints of the flux rope. This submergence, observable as flux cancellation, was found to occur shortly after the emergence of the same field lines. Nevertheless, this cancellation event was found to not directly contribute to the twist of the flux rope, even though it occurred in close vicinity to the flux rope. This analysis raises awareness of the need to take caution so as to not mistake correlation for causation when analyzing observations.

Thick-photosphere tether-cutting (TPTC) is the most significant event analyzed in this study. This event combined the already formed flux rope with additional untwisted arcades (see Sect. 4.5). Together, they formed an even longer flux rope, as well as some shorter field lines that were later submerged. This event is the realistic counterpart to the more ideal and symmetrical reconnection proposed by van Ballegooijen & Martens (1989). It is also similar to the observed reconnection of two loops to a longer loop in Wang & Muglach (2013). However, it is the first time that such an event has been simulated in a thick photosphere in a convectively driven MHD simulation with a gradually forming stable flux-rope. We find that the reconnection occurs within an extended volume of the thick photosphere; hence, the TPTC name. The reconnection also occurs over an extended period of time during the convergence of two eventually canceling polarities. In other words, the flux cancellation is clearly related to the reconnection and does not only appear in more artificial models. By simulating the convection zone, we also follow the shorter field lines, which were formed during the reconnection, down below the photosphere after the submergence. This is how such a TPTC event occurs on the Sun.

A flux rope did form in this simulation, in spite of the fragmented nature of the large-scale PIL that we introduced in the middle of the simulation box. That an atmospheric flux rope was formed is the same end result as in the planar line-tied simulations, the non-convective flux-emergence simulations, and the theoretical illustration by van Ballegooijen & Martens (1989). However, these past models all have smooth elongated PILs, along which reconnection happens all at once, to form the flux rope. Many observed solar prominences are not formed above smooth PILs. Hence, a question that we wanted to answer with this work was whether the Sun could form flux ropes above patchy canceling PILs. It can.

Data availability

Movie associated with Fig. 3 is available at https://www.aanda.org

Acknowledgments

We would like to thank M. Szydlarski and M. Carlsson for help setting up and running Bifrost simulations, M. Druett for help tracking the field lines with corks in Bifrost, as well as A. Prasad for visualization with Vapor. This research has been supported by the European Research Council through the Synergy Grant number 810218 (“The Whole Sun”, ERC-2018-SyG) and by the Research Council of Norway through its Centres of Excellence scheme, project number 262622 (Rosseland Centre for Solar Physics – RoCS). The simulations were performed on resources provided by Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway. This work was supported by the Action Thématique Soleil-Terre (ATST) of CNRS/INSU PN Astro, co-funded by CNES, CEA and ONERA.

References

Appendix A: Relaxing a reconnecting simulation

To prepare, run, and analyze the simulation presented in this paper caused unforeseen challenges. Numerical experimentation was therefore necessary to get the simulation to relax to a realistic model of the process under study. In this appendix, we add an explanation for some of the modifications done to the simulation, described in Sect. 2.2, to combat these unforeseen challenges. While not critical for the physical conclusions drawn in this paper, it could be relevant for someone attempting to run a similar simulation in the future.

At first, the flux feeding through the lower boundary was turned off, to not affect the magnetic field topology and reconnection we later wanted to study. After this modification, the simulation needed time to relax to the new BCs. The spatial resolution was reduced at the same time to conserve numerical resources, using the backstaff tool. Unexpectedly, the modified resolution caused enhanced p-mode waves through the box originating at the lower boundary. Furthermore, these waves were not fully canceled at the original coronal boundary at −14.3 Mm, but were in part reflected back into the box after about 10 min, which stirred the interior of the simulation box.

To minimize the reflections from the coronal boundary, it was attempted – mid-simulation, in the copied state of the QS simulation – to change the vertical BCs. Instead of extrapolating into the ghost zones and using the standard time-derivative equations, we tested characteristic BCs. However, that led to nonphysical point heating in the coronal boundary, likely as a consequence of the abrupt change of the boundary equations. The combined impact of the waves and their reflection is believed to have caused heating, possibly due to reconnection, high up in the corona in early tests of this simulation.

To rather minimize the push causing the wave at the lower boundary in the first place, we ended up with a modified lower BC instead. The simulation drove the average mass and energy densities over time to keep them close to a set of desired values. The reduction of the resolution shifted the depth of the lower boundary, which made the set of desired values inappropriate at the new depth. Hence, the desired values for these parameters at the lower boundary were reset to their new calculated averages just after modifying the resolution. This minimized the kick which drove the problematic waves, alleviating the problem.

When we later introduced the LFFF, a process described in Sect. 3.2, more waves occurred within the box as the simulation was again trying to relax to a self-consistent configuration including this new field. When we inserted the LFFF all at once, it created strong waves, which stirred and mixed the interior of the simulation box. When we rather ramped the extra field in 120 smaller insertions over 20 min, the simulation box had time to readjust, which greatly reduced this problem.

When the waves generated by ramping the LFFF reached the coronal boundary, which in part reflected them, it would sometimes cause numerical instabilities. It was found that the magnetic field of the added LFFF, when convected into the coronal boundary, again caused isolated heating, which the wave method in the Spitzer module was unable to convect away correctly, even at very short time steps. This problem caused further studies into the Spitzer module by Cherry et al. (2024), Furuseth et al. (2024). These studies found that the Courant condition for the wave method in fact gives a lower limit on the time step, below which, the method becomes unstable. This is opposite most explicit methods, which typically has an upper limit for the time step.

The challenges near the coronal boundary, both when relaxing the QS simulation and when ramping the LFFF, were alleviated by extending the corona almost twice as far away from the photosphere. This was done by adding horizontally stratified layers. However, as a consequence of adding stratified layers, the simulation required to be relaxed for longer before we could ramp the LFFF. These simulations with an extended corona, typically reached an approximately steady state after about 700 snapshots, or about 2 h. We relaxed the simulation for almost one more hour after this before we started to ramp the LFFF.

With the presented openness and transparency of this effort, it may seem that the gradual formation and twist of a flux rope presented in this paper was a stroke of luck, which finally occurred after hundreds of attempts. That is not the case. Similar flux ropes occurred in several earlier simulations that had, for example, nonphysical strong waves through the domain due to the change of the resolution, nonphysical strong reconnection at the coronal boundary due to the BC, and a lower resolution of Δx = 125 km.

Appendix B: Tracking an evolving flux rope in 3D

A flux rope is easy to visualize in single snapshots with both home-made numerical analysis tools and more professional tools such as Vapor (Li et al. 2019). However, to display how a flux rope evolves over time in 3D is challenging, particularly in this kind of simulation where the motion of the field lines are not controlled by the convection at the boundary (photosphere), but rather by self-consistent plasma convection.

We tracked the motion of the field lines with corks, point particles that were convected with the plasma. They worked very well to track the field lines over shorter periods of time, but over longer periods, they could be convected far away from the area of interest. That made it impossible to, for example, track the full history of the field lines over the 2.5 h of solar time shown in Fig. 4. Therefore, those figures were not made with corks. The slightly more than 40 min shown in Fig. 8 seems to be borderline, in our simulation, indicated by the magenta field lines convected down and out of the simulation box. The other colored families of field lines in this figure were also challenging to illustrate, as the corks were convected far along the field lines during that time interval. This challenge should be kept in mind if one were to plan a similar numerical experiment in the future.

All Tables

All Figures

Fig. 1. Horizontal averages of the Bifrost simulation. The magnetic pressures (PB) are given in dyn/cm2 for the QS simulation before the ramp at snapshot 1000 (dashed green), for the time-independent LFFF (dashed red), and for the simulation after the ramp at snapshot 1120 (dashed purple). The pressure (orange), the density (light blue), and the temperature (dark blue) of the atmosphere are plotted at snapshot 1000. The density is given in g/cm3, multiplied with a factor 1012 to fit on the same axis as the magnetic pressures, while the temperature is given on the right axis. The shaded gray area (z ∈ [0, −1] Mm) represents the “thick” photosphere, up to the altitude of the temperature minimum, i.e., the height where the density stratification is the strongest (see Sect. 3.4 for more details).

In the text

	Fig. 2. 2D cross sections of the magnetic field at y = 10 Mm. The streamlines show field lines, but their density does not signify the field strength. The field strength is shown by the color. (a) is the QS field right after the relaxation, (b) is the LFFF, (c) is the sum of the two proceeding, while (d) shows the EN right after the 20 min long ramp of the LFFF into the QS simulation.
In the text

Fig. 3. Magnetograms at the photosphere (z = 0) in the Bifrost simulation. (a) is just after the relaxation of the QS. (b) is just after the ramp of the LFFF. The white polarity (north pole) signifies flux coming out of the solar surface, being negative here as z is positive along the LOS. The black polarity (south pole) signifies flux entering the solar surface. The associated movie is available online

In the text

Fig. 4. Gradual formation and evolution of a flux rope in the Bifrost simulation. Each subfigure contains two panels, an angled view above a vertical top view of the same set of field lines. The angled views are identical in all subfigures, oriented such that the x axis increases down to the right and the y axis increases down to the left. The numbers on the axes are in units of Mm. (a) is right after the end of ramp of the LFFF, and differs from the rest by only showing field lines seeded at x = 12 Mm, y ∈ {5, 10, 15} Mm, and z ∈ (0, −5) Mm in gradually brighter colors. (b)–(f) show field lines representative for the flux rope in yellow-red colors, darker colors are shown for larger values of |J|/|B|, and overlying arcades in cyan. (d) also shows the four poles W1, W2, B1, and B2 referred to in the text.

In the text

Fig. 5. Slipping reconnection in the Bifrost simulation. Combined, the subfigures labeled (a)–(f) show the evolution until the time ta = 245m40s when this slipping process is complete. Each subfigure contains a top panel with an angled view of three families of 3D field lines (cyan, green, magenta), tracked by corks, above a 2D magnetogram. The angled view is angled from the bottom right corner of the 2D magnetogram toward the top left corner. On the bottom of each subfigure is a panel with only the 2D magnetogram, seen from directly above. The red circles in the bottom panels of each subfigure all highlight the same area of interest.

In the text

Fig. 6. U-loop emergence in the Bifrost simulation. From left to right, there are three different times relative to tb, a representative time for the emergence: (a) 100 s earlier, (b) tb, (c) 200 s later, after the end of the relevant flux cancellation. The upper panels of each subfigure show the 3D field lines of interest in green, tracked by corks, and pale yellow field lines representing the core of the flux rope, above a 2D magnetogram. The view in the upper panels is angled toward lower y. The lower panels show the same 2D magnetogram, from directly above. The cancellation of interest occurs within the red circle.

In the text

Fig. 7. Ω-loop submergence in the Bifrost simulation. From left to right, there are three different times relative to tc, a representative time for the submergence: (a) 300 s earlier, (b) tc, (c) 250 s later, after the end of the flux cancellation. The top panel of each subfigure shows three families of 3D field lines of interest in different colors (cyan, green, magenta), tracked by corks, and pale yellow field lines representing the core of the flux rope, above a 2D magnetogram. The middle panel shows the same, but with a partially transparent photospheric magnetogram, to reveal the field lines below the surface. The bottom panel shows only the 2D magnetogram, where the cancellation of interest occurs within the red circle. We stress that the two upper rows are viewed from the right, toward lower x, to better visualize the dynamics of interest.

In the text

Fig. 8. TPTC reconnection in the Bifrost simulation. In this figure, time advances downward, unlike in the previous figures, while all three panels per subfigure row show different views of the same snapshot. The left panels show four families of 3D field lines of interest in different colors (cyan, green, magenta, yellow), tracked by corks, above a 2D magnetogram. The magnetogram in the left panel in (c) is transparent by choice to show the magenta field line below the surface. The middle panels show the same, but from directly above. The right panels show only the 2D magnetogram, from directly above, where the cancellation of interest occurs within the red circle. There is no red circle in (d) as the cancellation then has completed.

In the text