MHD Simulation in the Solar Corona to Obtain Conditions for the Acceleration of Cosmic Rays during Solar Flares

1. Introduction

Solar cosmic rays are generated during the primordial energy release in solar flares. This explosive process takes place in the solar corona above the active region at altitudes of 15 000 to 70 000 km (~1/40–1/10 of the solar radius). This has been proven by direct measurements of the thermal X-ray emission of flares on the limb [1]. There is a number of other evidence for the energy release of a flare in the corona, based on an analysis of observations (see, for example, [2,3]). The possible positions of the flares found from the results of numerical simulations in the corona above the active region confirm the appearance of flares at these altitudes.

The slow accumulation of magnetic energy that occurs in the solar corona, and then its fast release during a flare, is explained by the mechanism according to which, during a flare, the energy accumulated in the magnetic field of the current sheet is released [4,5]. The current sheet, which is formed in the vicinity of a singular line of the X-type magnetic field as a result of the accumulation of disturbances, is first in a stable state. Further, in the process of quasi-stationary evolution, during which the plasma flows fast out of the sheet under the action of the magnetic tension force, the current sheet transfers into an unstable state. The instability in the current sheet causes a flare release of energy with all the observed manifestations, which are explained by the electrodynamic model of the flare proposed by I.M. Podgorny [6]. The model was developed based on the results of observations and numerical MHD simulation and uses analogies with the substorm electrodynamic model proposed earlier by its author [7].

Solar cosmic rays appear as a result of acceleration of charged particles, mainly protons, by an inductive electric field in the current sheet, equal to the field E=-V×B/c near the current sheet [8]. This inductive electric field is caused by a rapid change in the magnetic field during the current sheet instability. I. M. Podgorny noted an important property of this electric field for a wider class of physical phenomena: it is equal to E=-V×B/c in the field B of the powerful current formed near its location and in the plasma accelerated by magnetic pressure towards the location of the current. Besides the current sheet such an electric field arises in the pinch [9] created in an installation for solving the problem of controlled thermonuclear fusion. In both processes, charged particles are accelerated by an inductive electric field. For solar flares, they (mainly protons) under favorable conditions can leave the magnetic field above the active region, which will lead to the appearance of solar cosmic rays (SCR).

This mechanism can accelerate galactic cosmic rays during superflares on stars to energies five orders of magnitude higher than the SCR energy [10].

To study the mechanism of solar flares, the process of formation of a current sheet and the accumulation of energy in its magnetic field for each specific flare, and to obtain the possibility of improving the prognosis of solar flares based on an understanding of their physical mechanism, it is necessary to carry out magnetohydrodynamic (MHD) simulation of a flare situation in the solar corona above a real active region. The electric and magnetic fields obtained by MHD simulation above a real active region at the site of the primordial flare energy release and near it are necessary for studying the generation of solar cosmic rays. The acceleration of charged particles and the possibility of their exit from the acceleration region must be studied by calculating the trajectories of charged particles in electric and magnetic fields obtained by MHD simulation. The first results of studying the acceleration of charged particles were obtained by calculating trajectories in electric and magnetic fields, found by MHD simulation above the active region under simplified conditions [11].

Here, MHD simulation is continued above the real active region in order to more accurately determine the configurations of the magnetic and electric fields for studying SCR generation by calculating the trajectories of charged particles. More complex field configurations are obtained and analyzed near singular lines at the places where current sheets appear. A more detailed comparison of the MHD simulation results with observations has led to some problems, the solution of which requires further development of the numerical simulation technique.

When setting the conditions for MHD simulation, no assumptions were made about the physical mechanism of the solar flare; the boundary conditions were taken from observations. For MHD simulation in the real scale of time, a significant computational speed is required, which is achieved only through the use of parallel computing on supercomputers. Parallelization of calculations was carried out by computational threads on graphic cards (GPU). Experience has shown that the physical mechanism of a solar flare can be studied by MHD simulation only if the calculation starts a few days before the appearance of flares, when magnetic energy for the flare has not yet accumulated in the corona. Simulation in this setting of the problem showed the appearance of a current sheet above the active region [12]. Also, the appearance of a current sheet in such a setting of the problem was demonstrated by MHD simulations carried out by C. Jiang, S.T. Wu et al [13,14]. MHD simulations in the solar corona in other problem settings were performed in [15,16,17,18].

2. Methods of MHD Simulation

MHD simulation is carried out above the active region of AR 10365. The computational domain in the corona is a rectangular parallelepiped (0≤х≤1, 0≤y≤0.3, 0≤z≤1) (the length unit was chosen L₀=4×10¹⁰ cm). The lower boundary of the computational domain y=0 (XZ) is located on the surface of the Sun (photosphere) and contains the active region; the Y axis is directed from the Sun perpendicular to the photosphere. The three-dimensional system of MHD equations is solved numerically in a dimensionless form:

A typical magnetic-field strength in the active region, B₀ = 300 G, was used as the unit of the magnetic field. The units of the plasma density and temperature were taken to be their typical values in the corona above the active region, n₀=ρ₀/m_i=10⁸ cm^-3, T₀=10^{6 o}K, which were assumed to be initially constant (m_i = 1.67×10^-24 g – ion mass). The units of plasma velocity and time are, respectively, the Alfvén velocity and Alfvén time $V_0=V_A=B_0\sqrt{4\pi {\rho }_0}$ = 0.65×10¹⁰ cm/s, t₀=L₀/V₀= 6.1 s. In equations (1) - (4) γ=5/3 is adiabatic constant, Re_m=V₀L₀/ν_m0 is the magnetic Reynolds number, ν_m0=с²/4πσ₀ is the magnetic viscosity for the conductivity σ₀ at the temperature T₀, σ is conductivity, σ/σ₀=T^3/2, β=8πn₀kT₀/B₀2 (n₀=ρ₀/m_i, m_i is the ion mass). Re=ρ₀L₀V₀/η is the Reynolds number, η is the viscosity. G_q=L(T₀)ρ₀t₀/T₀, L'(T)=L(T₀T)/L(T₀) is the dimensionless radiation function. κ_dl=κ/(Πκ₀) is the dimensionless thermal conductivity along the magnetic field lines, Π=ρ₀L₀V₀/κ₀ is the Peclet number, κ₀ is the thermal conductivity for the temperature T₀, κ is the thermal conductivity, κ/κ₀=T^5/2.

κ_⊥dl=[(κκ₀^-1Π^-1)(κ_Bκ_0B^-1Π_B^-1)]/[(κκ₀^-1Π^-1)+(κ_Bκ_0B^-1Π_B^-1)] is the dimensionless thermal conductivity perpendicular to the magnetic field lines, Π_B=ρ₀L₀V₀/κ_0B is the Peclet number for the thermal conductivity across the strong magnetic field (when the cyclotron radius is much smaller than the free path length). This thermal conductivity is denoted as κ_B, and κ_0B is its value for the temperature T₀, plasma density ρ₀, and field B₀; κ_B/κ_0B=ρ²B^-2T ^-1/2.

To select the parameters, the principle of limited simulation [19] was used, according to which, much larger and much smaller units of dimensionless parameters are set in calculations much larger and much smaller than unity without accurately preserving their values. The main parameters are the magnetic and ordinary Reynolds numbers Re_m and Re. Dimensionless values of artificial viscosities (ν_{m_Art_Ph} и ν_{Art_Ph}) were set to stabilize the numerical instabilities.

For the numerical solution, an upwind absolutely implicit finite-difference scheme was developed, which is conservative with respect to the magnetic flux [20]. The scheme is solved by the method of iterations. The methods were applied with the aim of constructing a scheme that remains stable for the maximum possible time step. The scheme is realized in the computer program PERESVET.

A stable finite difference scheme (high stability is achieved if this scheme is implicit) should approximate the magnetic induction equation (1), in which the diffusion term is taken as ν_mΔB rather than -rot(ν_mrotB). For this purpose, it is necessary to preserve the magnetic flux through the boundary of each cell of the difference scheme, which means that the finite-difference analog [divB] is equal to zero with high accuracy, if the values are correctly specified on the grid elements of the difference scheme. For a difference scheme conservative with respect to the magnetic flux, the finite-difference analogs of [ν_mΔB] and [-rot(ν_mrotB)] coincide. If in the finite difference scheme the magnetic flux is conserved with high accuracy, then the scheme with finite difference approximation of the diffusion term in the form [-rot(ν_mrotB)] is also stable, since this diffusion term is equal to [ν_mΔB] with high accuracy. In addition, when solving a finite difference scheme by iteration method, the values at the central point of the template can be taken both at the current iteration and at the next iteration. It becomes possible to use various variants of the difference scheme. The above considerations for the constant magnetic viscosity ν_m extend to the case where ν_m depends on spatial coordinates. To do this, it is necessary to add finite-difference analogs of the products of the spatial derivatives of ν_m and the spatial derivatives of the B components to the scheme, which do not cause instability of the difference scheme.

For the numerical solution of the MHD equations, several versions of the scheme were proposed, in which the magnetic flux is conserved with high accuracy. A more accurate visual representation of the arrangement of quantities on a three-dimensional grid of a finite-difference scheme, conservative with respect to the magnetic flux, is done. Images of the location of the values on the grid are shown in Figure 1 and Figure 2. The finite-difference analogue of the magnetic field divergence [divB] for such a scheme should be equal to zero with high accuracy, for our calculations the dimensionless value |[divB]| ~ 10^-4. In connection with the appearance of numerical instabilities near the boundary of the region, additional test calculations were carried out for possible corrections in the choice of the scheme variant with the conservation of the magnetic flux. As a result, taking into account the influence of instabilities near the boundary, the variant in which the iterations converge the fastest remains the most efficient. In this version, the diffusion term in the equation for the magnetic field was taken in the form of a finite-difference analog ΔB ([ΔB]), rather than rot(rotB) ([rot(rotB)]), and the values at the central point of the grid template of the finite-difference scheme were taken on the next iteration.

In the finite-difference scheme, which is conservative with respect to the magnetic flux, the magnetic field is specified not by the field vectors at the grid points, but by the magnetic fluxes through the boundaries of the grid cells, averaged per unit area (Figure 1). Used to specify the magnetic field components, the magnetic fluxes through the cell faces, averaged per unit area, approximate the components of the magnetic field vector perpendicular to the given face, taken at the central point of the given face. The values of plasma density ρ, velocity V, and temperature T were set at the grid points of the finite-difference scheme, i.e., at the vertices of the cubes, which are grid cells.

The dimensionless electric field E= -V×B + ν_m j (here ν_m = Rm^-1 (σ₀/σ) = Rm^-1 T^-3/2 is the dimensionless magnetic viscosity) and the dimensionless electric current density j=rotB is clearly are not included in the equations of magnetic hydrodynamics, but are expressed in terms of the quantities included in these equations. It is convenient to use the values E and j for constructing a difference scheme; they are determined on the edges of the grid cells. The edges of the grid cells, which can be parallel to one of the coordinate axes X, Y or Z, are numbered by the three-dimensional number i, j, k of the end of this face with the minimum value of the corresponding coordinate.

The equation for the magnetic field (1) in the introduced notation has the form $\frac{\partial B}{\partial t}=\mathrm{rot}(-E)$. The flux of the magnetic field B through the surface of a small area S, perpendicular to the vector B is equal to Φ_B= B S, (if there is an arbitrary surface of a large area, then the magnetic flux through it can no longer be expressed through the field vector at any one point, the magnetic flux will be equal to the integral over the surface ${\mathsf{\Phi }}_{\mathrm{B}}={\displaystyle \underset{S}{\iint }(\mathrm{B},ds)}$, where B is the magnetic field vector defined at each point of the surface, ds is the vector perpendicular to the surface, the value of which is equal to the small area ds on the surface under consideration, and (B,ds) is the scalar product of these vectors).

The flux of the vector rot(-E) through the surface of a small area S, perpendicular to the vector rot(-E) is equal to Φ_rot(-E) = rot(-E) S. According to the induction equation $\frac{\partial B}{\partial t}=\mathrm{rot}(-E)$, this flow Φ_rot(-E) represents is a change in time of the magnetic flux Φ_B. Since the divergence of the magnetic field must always be zero divB=0, the magnetic flux through a surface bounding any region in space will be zero. This confirms the induction equation, from which it follows that if for an initial magnetic field its flux through the surface bounding any region in space is zero, it will remain zero for any moment in time, since the flux of the vector field rotor (in this case rot(-E) ), representing the change in time of the magnetic flux, is always zero.

To improve stability, a scheme is constructed in which the flow through the grid cell boundary remains equal to zero with very high accuracy (with the accuracy of iteration convergence if this is an implicit scheme solved by the iteration method, or with the accuracy of representing numbers in a computer if an explicit scheme is used, which less stable). This is a scheme that is conservative with respect to magnetic flux.

The magnetic flux through a face perpendicular to the X axis is equal to the product of the magnetic field component B_x, defined on this face, by the area of this face

Φ_{x, i, j, k} = B_{x, i, j, k} S_x = B_{x, i, j, k} h_y h_z (h_y, h_z - grid steps along the Y and Z axes, Figure 1. bottom row).

Similarly, magnetic fluxes are determined through faces perpendicular to the Y and Z axes. The total flux of the magnetic field emerging through the surface limiting the volume of the computational grid cell with numbers i+1, j+1, k+1 (the cell shown in Figure 1 with the limits of coordinate change (x_i≤x≤x_i+1, y_j≤y≤y_j+1, z_k≤z≤z_k+1) ) can be found knowing the flows through the faces of this cell; for a scheme that is conservative with respect to the magnetic flux, this value should be equal to zero (with high accuracy):

Φ_{x, i+1, j+1, k+1} - Φ_{x,i, j+1, k+1} + Φ_{y,i+1, j+1, k+1} - Φ_{y, i+1, j, k+1} + Φ_{z,i+1, j+1, k+1} - Φ_{z,i+1, j+1, k} = 0

or, substituting the values of magnetic fluxes:

B_{x, i+1, j+1, k+1} S_x - B_{x,i, j+1, k+1} S_x + B_{y,i+1, j+1, k+1} S_y - B_{y, i+1, j, k+1} S_y + B_{z,i+1, j+1, k+1} S_z - B_{z,i+1, j+1, k} S_z =0

Where S_x S_y S_z are the areas of cell faces perpendicular to the vectors of the selected coordinate system e_x e_y e_z : S_x= h_yh_z ; S_y=h_xh_z ; S_z=h_xh_y

The density of the magnetic flux source in a grid cell is the value of the total magnetic flux divided by the volume of the cell. Dividing the total magnetic flux by the volume of the grid cell V=h_xh_yh_z we get

(B_{x, i+1, j+1, k+1} - B_{x,i, j+1, k+1} )/h_x + (B_{y,i+1, j+1, k+1} - B_{y, i+1, j, k+1})/h_y + (B_{z,i+1, j+1, k+1} - B_{z,i+1, j+1, k})/h_z = 0

i.e. finite-difference analog of the magnetic field divergence divB ([divB]) vanishes, as expected

Equality to zero of the total flux of the magnetic field through the cell boundary in a scheme that is conservative with respect to the magnetic flux is achieved by setting the initial magnetic field on the faces of the grid cells with zero total magnetic flux and zero change in the flux with time. The change in the magnetic flux, representing the flux of the vector rot(-E), which goes through the boundary of the grid cell, will be zero if, in the scheme of the components of the vector rot(-E), the finite-difference analogs of the derivatives are expressed in terms of the components of the vectors, the electric field E, specified on the edges of the cells grids, as illustrated in Figure 2.

Despite the use of special methods, simulation in the real scale of time is possible only with the help of parallel computing. They were carried out by parallel computing threads on graphics cards using CUDA technology.

The main problem of MHD simulation above a real active region is the numerical instabilities that arise near the boundary of the computational region. The developed methods for stabilizing the instability made it possible to solve partially the problem for calculations with relatively low viscosities (Re_m= 10⁹, Re= 10⁷), which practically do not suppress the perturbation propagating from the photosphere. The results obtained make it possible to determine ways to further improve the technique for stabilizing numerical instabilities.

3. Results

The configuration of the magnetic field obtained by MHD simulation is so complex that it is often impossible to determine from this configuration the positions of singular lines and the current sheets appearing near them from it. For this purpose, a graphical search system of flare positions has been developed [12] and then modernized. The system is based on the search for current density maxima, which are reached in the middle of the current sheets, and further analysis of the magnetic field configuration near the points of the found maxima.

Figure 3, Figure 4 and Figure 5 present a comparison of the results of MHD simulations with observations of radio emission at a frequency of 17 GHz obtained with the Nobeyama radioheliograph (NoRH) before the M 1.9 flare on May 26, 2003 above the active region of AR 10365 and during this flare (2003- 05-26 04:47:05, 05:32:05 and 05:41:05). The positions of the current density maxima are marked with numbers 1, 2, 3, 4, 5, 11, 13, 14, 15, 16, 62, 64, 97, and 115, the remaining maxima are indicated by green dots (the current density maxima are renumbered in decreasing order of density current). A divergent magnetic field is superimposed on the X-type magnetic field at many current density maxima. However, due to the presence of the X-type configuration, current sheets are formed at such locations even if the divergent field dominates and a deformed divergent magnetic field is obtained as a result of the superposition.

4. Discussion

Comparison of the results of MHD simulations with observations before and during flares showed a general agreement between the positions of the current sheets found by our graphical search system and the regions of intense flare radiation. A significant number of points at which current density maxima are reached are located in the region of strong radiation. Also, a significant part of the points of maxima is located at a relatively small distance from the region of strong radiation (~10 Mm and less), which can be explained by the error of the numerical method and physical processes during the flare. These results confirm the mechanism according to which the energy accumulated in the magnetic field of the current sheet is released during the flare.

There are some problems with the details of such coincidences. Especially if you choose current sheets in which the current directions are almost parallel to the solar surface and they are located at the tops of magnetic loops (or at least near the tops of magnetic loops) at a fairly high altitude of 25,000 km - 45,000 km. These heights are lower than 60,000 - 70,000 km, where it is difficult to expect the accumulation of large magnetic energy due to low magnetic fields, but at the same time, these are not heights of 15,000 - 20,000 km, where, apparently, a current sheet may not appear with a sufficiently small longitudinal field so that explosive instability can occur in it. In the magnetic field of current sheets at the tops of magnetic loops at altitudes of 25,000–45,000 km, apparently, energy for a flare can accumulate best of all, and after a quasi-stationary evolution, such current sheets can go into an unstable state.

The results obtained indicate the need for work to solve the problems that have arisen. It may be necessary to carry out more accurate calculations. This can be helped by a more efficient stabilization of the instability near the boundary of the computational domain, as well as a more accurate setting of the boundary conditions, for which there are also some possibilities.

5. Conclusions

• The methods of MHD simulation in the solar corona above a real active region have been improved. The application made it possible to perform calculations above the active region of AR 10365 during a time interval of about three days, but showed the need for further modernization. For the upwind, absolutely implicit finite-difference scheme, which is conservative with respect to the magnetic flux, developed and implemented in the PERESVET program, a variant of magnetic flux conservation was chosen, for which the stabilization of the numerical instability that occurs near the boundary of the computational domain is best performed. The parallel computing algorithm implemented by computational threads on graphics cards using CUDA technology has been modernized over the past few years, after its first version was developed. The stabilization of the numerical instability that occurs near the boundary of the computational domain has been improved by more accurate selection of artificial viscosity values and by using other methods.
• An analysis of the new results of MHD simulation before solar flares and during solar flares confirmed the formation of singular X-type lines in the magnetic field configuration and the appearance of current sheets in them, in the magnetic field of which energy for solar flares is accumulated.
• The results obtained confirmed the existence of a large number of singular magnetic field lines, in the vicinity of which a divergent magnetic field is superimposed on the X-type field configuration. Even if the divergent magnetic field dominates such that a deformed divergent magnetic field results from the superposition, the presence of an X-type field can produce a sufficiently strong current sheet to explain the appearance of a medium-power flare.
• The results of the analysis of MHD simulation data for the time moments before and during flares confirmed that a large number of current density maxima (more than a hundred), located on a singular magnetic field line in the vicinity of which the current sheet is formed, are in the region of intense flare radio emission at a frequency of 17 GHz, or at a distance from the region of intense radiation not exceeding 10 Mm. This arrangement of current density maxima confirms the solar flare mechanism based on the release of energy stored in the magnetic field of the current sheet. The details of the location of the current density maxima did not always correspond to the observational data, which indicates the need to develop more accurate methods for solving the equations of magnetohydrodynamics and setting the boundary conditions on the solar surface.