EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 638 (June 2020) A&A, 638 (2020) A79 Full HTML
Free Access
Issue
A&A
Volume 638, June 2020
Article Number A79
Number of page(s) 10
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202037809
Published online 15 June 2020

© ESO 2020

1. Introduction

When modeling astrophysical systems, the simplest magnetohydrodynamic (MHD) approximation is frequently used where the plasma is considered a single fluid with total coupling between its constituent microscopic species. This assumption is able to satisfactorily describe the physics of many phenomena in different astrophysical contexts; however, the approximation may no longer be valid when the plasma is partially ionized and ions and neutrals drift with respect to each other. This is the case for the interstellar medium (e.g., Spitzer 1978; Zweibel 2002), molecular clouds (e.g., Zweibel & Josafatsson 1983; Padoan et al. 2000; Basu & Ciolek 2004; Crutcher 2012), protoplanetary disks (e.g., Wardle 1999; Salmeron & Wardle 2008; Gressel et al. 2015; Tomida et al. 2015), star formation (e.g., Mestel & Spitzer 1956; Shu et al. 1987; Kudoh & Basu 2008), and the solar chromosphere (e.g., Goodman 2004; Zweibel et al. 2011; Khomenko & Collados 2012; Martínez-Sykora et al. 2015; Zweibel 2015; Shelyag et al. 2016), among others.

It is possible to relax the MHD approximation to deal with partially ionized gases, considering the relative speed and associated friction between neutrals, ions, and electrons, and still treating the plasma as a single fluid: the Generalized Ohm’s Law (see the seminal books by Braginskii 1965; Mitchner & Kruger 1973; Cowling 1976). This way, the departure of the MHD approximation can be handled by just extending the induction and energy equations by adding the ambipolar diffusion term, which concerns the decoupling of neutral and charged components, and the Hall effect, which takes the drift velocities between ions and electrons into account. This extension has been applied in different codes by Mac Low et al. (1995), Leake et al. (2005), O’Sullivan & Downes (2007), Cheung & Cameron (2012), Martínez-Sykora et al. (2012), Masson et al. (2012), Tomida et al. (2015), González-Morales et al. (2018), Grassi et al. (2019), among others, and has been shown to be important to better understand the role of the ambipolar diffusion and Hall terms in astrophysics. However, the inclusion of partial ionization effects into advanced numerical codes confronts the modeler with different difficulties. On the one hand, the importance of the new effects sensitively depends on the microscopic constitution of the plasma, namely, on the abundances, the chemistry, the ionization degree and the collisions between different species. For instance, in the solar atmosphere, Martínez-Sykora et al. (2012) showed that the approximation chosen to determine the values of collision cross sections and frequencies is crucial for ion-neutral interaction effects: there are significant discrepancies in the ambipolar diffusion coefficient depending on the assumption considered that lead to different results for the thermal properties, primarily in the chromosphere. In protostellar disc formation, Zhao et al. (2016) found that reducing the number of very small grains enhances ambipolar diffusion. In molecular clouds, Grassi et al. (2019) showed that cosmic rays can impact on the ionization level of the molecular gases, thus modifying the importance of the ambipolar diffusion. Those are a few examples of how the inclusion of proper physics is essential to obtain a realistic outcome when addressing partially ionized plasma. On the other hand, the computations including partial ionization effects (albeit addressed from a single-fluid approach, thus avoiding the complexity of multifluid equations; see, e.g., Leake et al. 2012; Alvarez Laguna et al. 2016) are very slow when they are solved through explicit methods, due to the strong constraints with respect to the timestep. According to the Courant-Friedrichs-Lewy (CFL) criterion (Courant et al. 1928), the maximum timestep, ΔtCFL, for the numerical solution of parabolic (e.g., diffusion) problems using explicit schemes decreases as the square of the spatial resolution Δx (i.e., ΔtCFL ∝ Δx2/D, where D is the diffusion coefficient). Wherever high spatial resolution is required, this quadratic dependence can strongly limit the calculation speed in comparison with non-diffusive MHD computations, whose Δt is linearly dependent on Δx. As a consequence, high-resolution experiments of diffusion problems are virtually impossible to perform explicitly. Different strategies have been carried out to alleviate this problem. For instance, Nakamura & Li (2008), Li et al. (2011), Masson et al. (2012) use different thresholds to decrease the ambipolar term to avoid strongly restrictive timesteps when needed. Other authors, such as Mac Low et al. (1995), Mellon & Li (2009), Leake & Arber (2006) adopt a sub-cycling method in which the induction equation is evolved separately from the rest of MHD equations when the timestep corresponding to the ambipolar diffusion is smaller than the dynamical timestep. An extension of this method is used by Martínez-Sykora et al. (2012); Martínez-Sykora et al. (2017a,b) to also consider sub-cycling the evolution of the energy equation because of the dissipation due to ambipolar diffusion. Another technique is super time stepping (STS; Alexiades et al. 1996), which allows the restrictive CFL criterion to be relaxed to speed up the explicit calculation of parabolic problems. This technique was shown to efficiently accelerate heat conduction calculations (see also Meyer et al. 2012; Iijima & Yokoyama 2015), and since then it has been extensively used in ambipolar diffusion contexts (Choi et al. 2009; Commerçon et al. 2011; Tomida et al. 2015; Gressel et al. 2015; González-Morales et al. 2018). The drawback is that the STS method has two free input parameters, so it is necessary to carefully choose their values not only to optimize the performance, but also to avoid the destabilization of the scheme which may lead to meaningless results (for more details about numerical approaches in partially ionized systems see the recent review by Ballester et al. 2018).

The purpose of this paper is to introduce a numerical tool that allows us to confront the numerical challenges due to ambipolar diffusion in the solar atmosphere. To this end, we have developed a new module in the Bifrost code (Gudiksen et al. 2011), taking care, among other things, of the number of elements included in the calculations and their ionization state. Due to the numerical stiffness imposed by the ambipolar diffusion, its numerical implementation must be efficient to be able to calculate complex problems in which high resolution is mandatory.

The layout of this work is as follows. Section 2 details the relevant equations of the Generalized Ohm’s Law to establish the context for the subsequent parts of the paper. Section 3 describes the implementation of the collision cross sections and frequencies necessary to compute the ambipolar diffusion coefficient. Section 4 addresses the computation of the ionization state for the ambipolar diffusion term when assuming local thermodynamic equilibrium (LTE), or nonequilibrium (NEQ) ionization and recombination of hydrogen and helium. In Sect. 5 we explain the STS method, together with its implementation in the Bifrost code. Section 6 contains the recipes of the hyperdiffusion terms to guarantee stability for the code. Section 7 presents the validation test. Finally, Sect. 8 summarizes the main conclusions of the present work.

2. Generalized Ohm’s law

The generalized Ohm’s law is a relation between the electric field and the electric current that makes it possible to overcome the difficulties of dealing with multifluid plasmas, such as extremely high magnetic field mediated wave speeds, a large number of equations, and stiff systems (e.g., Ballester et al. 2018, and references therein). This is doable by using a one-fluid approximation that requires a high level of (but not infinite) coupling between neutrals and the charged species. In a reference frame locally moving with a plasma element, it can be shown that this relation is given by

E ′ = η J ′ − η amb ( J ′ × B ′ ) × B ′ | B ′ | 2 + η Hall ( J ′ × B ′ ) | B ′ | , $$ \begin{aligned} \boldsymbol{E^{\prime} } = \eta \boldsymbol{J^{\prime} } - \eta _{\mathrm{amb}}\frac{(\boldsymbol{J^{\prime} } \times \boldsymbol{B^{\prime} }) \times \boldsymbol{B^{\prime} }}{|\mathbf {B^{\prime} }|^2} + \eta _{\mathrm{Hall}}\frac{(\boldsymbol{J^{\prime} } \times \boldsymbol{B^{\prime} })}{|\mathbf {B^{\prime} }|}, \end{aligned} $$(1)

where E′ is the electric field, B′ the magnetic field, and J′ the current density all measured in that reference frame. The coefficient η is the standard ohmic diffusion given by

η = m e ( ν en + ν ei ) n e q e 2 ; $$ \begin{aligned} \eta = \frac{m_{\rm e} (\nu _{\rm en} + \nu _{\rm ei} )}{n_{\rm e} q_{\rm e}^2}; \end{aligned} $$(2)

the ambipolar diffusion coefficient, ηamb, is defined as

η amb = ( ρ N / ρ ) 2 | B ′ | 2 Σ n Σ i ρ n ν ni ∗ = η amb ∗ | B ′ | 2 ; $$ \begin{aligned} \eta _{\mathrm{amb}}= \frac{(\rho _N/\rho )^2|\mathbf {B^{\prime} }|^2}{\Sigma _{\rm n} \Sigma _i \rho _{\rm n} \nu ^{*}_{\rm ni}} = \eta ^{*}_{\mathrm{amb}}|\mathbf {B^{\prime} }|^2 ;\end{aligned} $$(3)

and the Hall coefficient ηHall is given by

η Hall = | B ′ | q e n e = η Hall ∗ | B ′ | . $$ \begin{aligned} \eta _{\mathrm{Hall}}= \frac{|\mathbf {B^{\prime} }|}{q_{\rm e} n_{\rm e}} = \eta ^{*}_{\mathrm{Hall}}|\mathbf {B^{\prime} }| .\end{aligned} $$(4)

Here qe is the electron charge; ne the density number of electrons; me the electron mass; ρn is the density of the neutral element n; ρN the total neutral mass density (ρN = Σnρn); ρ the total mass density; νen and νei the collision frequency of electrons with neutrals and ions, respectively; and ν ni * $ \nu^{\ast}_{\mathrm{ni}} $ is the reduced neutral-ion collision frequency, which we explain in Sect. 3. It is important to realize that to obtain Eq. (1), the relative velocity between charged species is considered negligible, and similarly between neutral species (see Zaqarashvili et al. 2011; Khomenko et al. 2014; Shelyag et al. 2016, among others). This is an assumption that has been widely made in the literature to simplify the multifluid equations to a single fluid equation, and it implies that the collision times between charged species and the neutral-neutral species must both be much smaller than the mixed charged-neutral collision times, and that all of the foregoing must be much smaller than the macroscopic timescales. On the other hand, we note that the focus of this paper is on the ambipolar diffusion term and its new implementation in the Bifrost code. The implementation of the Hall term has been presented in the papers by Martínez-Sykora et al. (2012); Martínez-Sykora et al. (2017b).

2.1. Induction equation

We now discuss the laboratory reference frame where the plasma element moves with a velocity u, and the electric field, the electric current and the magnetic field are given by E = E′ − u × B, J = J′, and B = B′, respectively. Using Faraday’s induction equation and Eq. (1), we obtain a generalized induction equation:

∂ B ∂ t = ∇ × [ u × B − η Hall ∗ ( J × B ) + η amb ∗ ( J × B ) × B − η J ] $$ \begin{aligned} \frac{\partial \boldsymbol{B}}{\partial t} = \nabla \times \Big [ \boldsymbol{u} \times \boldsymbol{B} - \eta ^{*}_{\mathrm{Hall}}(\boldsymbol{J} \times \boldsymbol{B}) + \eta ^{*}_{\mathrm{amb}}(\boldsymbol{J} \times \boldsymbol{B}) \times \boldsymbol{B} - \eta \boldsymbol{J}\Big ] \end{aligned} $$(5)

This induction equation contains two additional terms in comparison to the one from the classic resistive MHD: one term proportional to J × B associated with the Hall effect and another term proportional to (J × B)×B associated with the ambipolar diffusion. These terms cannot lead to changes in the magnetic topology, and consequently neither produce magnetic reconnection. This can be seen if we define

u Hall = − η Hall ∗ J $$ \begin{aligned} \boldsymbol{u}_{\rm Hall} = - \eta ^{*}_{\mathrm{Hall}}\boldsymbol{J} \end{aligned} $$(6)

and

u amb = η amb ∗ ( J × B ) $$ \begin{aligned} \boldsymbol{u}_{\rm amb} = \eta ^{*}_{\mathrm{amb}}(\boldsymbol{J} \times \boldsymbol{B}) \end{aligned} $$(7)

for the Hall and ambipolar terms, respectively. With these notations,

∂ B ∂ t = ∇ × [ ( u + u Hall + u amb ) × B − η J ] , $$ \begin{aligned} \frac{\partial \boldsymbol{B}}{\partial t} = \nabla \times \Big [ \big ( \boldsymbol{u} + \boldsymbol{u}_{\mathrm{Hall}} + \boldsymbol{u}_{\rm amb} \big ) \times \boldsymbol{B} - \eta \boldsymbol{J} \Big ] ,\end{aligned} $$(8)

which means that the magnetic field is no longer frozen into the plasma flow, but it is frozen into a pseudo-flow with speed u + uHall + uamb. Even though these terms cannot lead to changes in the magnetic topology, they can significantly change the behavior of reconnection, for example by a rapid thinning of the current sheet or an interplay with the plasmoid instability (see, e.g., Huang et al. 2011; Ni et al. 2015); when η = 0, the topology is preserved and there is no magnetic reconnection.

2.2. Energy equation

From basic electrodynamics, the power exerted by the electromagnetic field on the plasma is given by J′ ⋅ E′. Using the generalized Ohm’s law (Eq. (1)), in addition to the classic ohmic dissipation ηJ2, a new term appears that leads to an irreversible entropy increase, i.e.,

Q amb = η amb J ⊥ 2 , $$ \begin{aligned} Q_{\mathrm{amb}} = \eta _{\mathrm{amb}}J_{\perp }^2, \end{aligned} $$(9)

where J⊥ is the current component perpendicular to the magnetic field. This dissipation term is associated with the collisions between neutrals and ions, hereafter referred to as ambipolar diffusion heating. As a consequence, when taking ion-neutral interaction effects into account, a new entropy source has to be added to the energy equation as

∂ e ∂ t = − ∇ · ( e u ) − P ∇ · u + η J 2 + η amb J ⊥ 2 + Q rad + Q Spitz , $$ \begin{aligned} \frac{\partial e}{\partial t} = - \nabla \cdot (e \boldsymbol{u}) - P \nabla \cdot \boldsymbol{u} + \eta J^2 + \eta _{\mathrm{amb}}J_{\perp }^2 + Q_{{\mathrm{rad} }} + Q_{{\mathrm{Spitz} }} ,\end{aligned} $$(10)

where e is the internal energy per unit volume, P the gas pressure, Qrad represents all the entropy sources due to radiation, and QSpitz is the entropy source due to the thermal Spitzer conductivity. The Hall term does not cause any energy dissipation since it is perpendicular to J.

3. Collision frequency and cross section

Collision frequency and cross section are important quantities when determining the rate at which electrons and the different ions and neutrals interact with each other. For this reason, we have to take the state-of-the-art models and measurements of the mentioned parameters into consideration.

3.1. Collision frequency

The reduced neutral-ion collision frequency, ν ni * $ \nu^{\ast}_{\mathrm{ni}} $, is given by

ν ni ∗ = m ni m n n i σ ni ( 8 K B T π m ni ) 1 / 2 , $$ \begin{aligned} \nu ^{*}_{\rm ni} = \frac{m_{\rm ni}}{m_{\rm n}} n_{{i}} \sigma _{\rm ni} \left( \frac{8 K_{\rm B}T}{\pi m_{\rm ni}}\right)^{1/2} ,\end{aligned} $$(11)

where ni is the ion number density, KB the Boltzmann constant, T the temperature, and mni = mnmi/(mn + mi) the reduced mass of the neutral and ion species. In the solar atmosphere the most frequent collisions are those of the most abundant elements, H and He, with ions of other elements (Khomenko et al. 2014). Therefore, we consider the following ion-neutral interactions: neutral hydrogen (H) and helium (He) atoms colliding with singly ionized ions of the 16 most important elements (either high abundance or low ionization potential) in the Sun and electrons, as well as neutral hydrogen molecules (H2) colliding with protons (p) and electrons (e). Thus, the denominator of the ambipolar diffusion coefficient (Eq. (3)) in our case is

Σ n Σ i ρ n ν ni ∗ = ρ H ( Σ i ν H , i ∗ + ν H , e ∗ ) + ρ He ( Σ i ν He , i ∗ + ν He , e ∗ ) + ρ H 2 ( ν H 2 , p ∗ + ν H 2 , e ∗ ) $$ \begin{aligned} \begin{split} \Sigma _n \Sigma _i \rho _n \nu ^{*}_{\rm ni} = \rho _{\rm H} (\Sigma _i\nu ^{*}_{\mathrm{H} ,i} + \nu ^{*}_{\mathrm{H} ,e}) + \rho _{\rm He} (\Sigma _i\nu ^{*}_{\mathrm{He} ,i} + \nu ^{*}_{\mathrm{He} ,e}) \\ + \rho _{\rm H_2} ( \nu ^{*}_{\mathrm{H_2} ,\mathrm{p} } + \nu ^{*}_{\mathrm{H_2} ,\mathrm{e} } ) \end{split} \end{aligned} $$(12)

3.2. Cross section

The cross section for elastic scattering between a given neutral and a charged particle, σni, is implemented through tables calculated on the basis of values given by different authors:

Neutral hydrogen atoms with protons (H–p). From 10−4 to 102 eV in the center of mass of the collision (ECM), the cross section values are based on quantum-mechanical indistinguishability calculations, additionally including charge transfer (see Krstić & Schultz 1999a; Glassgold et al. 2005; Vranjes & Krstic 2013)1. In Fig. 1 we plot this cross section, as a solid red curve, in the range ECM = [∼ 0.1, 100] eV, which corresponds to typical solar temperatures T = [103, ∼106] K (1 eV = 11 604 K). For comparison purposes, we have also plotted, as a dashed red line, the constant cross section for collisions between neutral hydrogen atoms and protons used by other authors such as Khodachenko et al. (2004), Leake & Arber (2006), Soler et al. (2009), Khomenko & Collados (2012), among others. For temperatures below 104 K, the constant definition underestimates, by approximately an order of magnitude, the temperature-dependent cross section.

thumbnail Fig. 1.

Elastic scattering cross section, σni, implemented in the Bifrost code to calculate ηamb. The different lines correspond to H–p collisions including temperature dependency and charge exchange (solid red line), H–e (pink), He–p (yellow), He–He+ (blue), He–e (green), H2–p (black), and H2–e (gray). For comparison purposes, the constant σni curve has been added, corresponding to H–p collisions frequently used in the literature (dashed red line).

Neutral hydrogen atoms with electrons (H–e). The values are extracted from Vranjes & Krstic (2013) and references therein. In Fig. 1, this cross section is shown as a pink line from 0.1 to 10 eV.

Neutral helium atoms with protons (He–p). The elastic scattering cross sections are taken from Vranjes & Krstic (2013) and shown in Fig. 1 as a yellow curve.

Neutral helium atoms with singly ionized helium (He–He+). The values are obtained from the second chapter of the book by Franz (2009). In Fig. 1, this cross section is shown in blue from ∼4 to 100 eV.

Neutral helium atoms with electrons (He–e). These values are extracted from Vranjes & Krstic (2013) and references therein. In Fig. 1, this cross section is plotted in green from 0.1 to 10 eV.

H2 molecules with protons (H2–p). This cross section is based on the fully quantum-mechanical calculations by Krstic & Schultz (1999b) and is shown in Fig. 1 as a black line.

H2 molecules with electrons (H2–e). Values of this elastic scattering cross section are extracted from a data base compiled by Yoon et al. (2008) and shown in Fig. 1 as a gray line.

The collision cross sections for hydrogen (or helium) and heavier elements are not well known, so we adopt the same assumption as made by Vranjes et al. (2008): we take the cross section between hydrogen (or helium) and protons multiplied by mi/mH (or mi/mHe), where mi is the mass of the heavier element.

4. Ionization state

The number of elements considered and their ionization states are also crucial to properly determine the ambipolar diffusion coefficient. For this reason, the module presented in this paper has been developed to compute ηamb in a consistent way with the different existing available equation of state (EOS) modules in the Bifrost code. We have set up two methods to calculate ηamb according to the EOS used: obtaining the ambipolar diffusion coefficient through precomputed tabulated values, which is only possible when using an EOS that assumes LTE (Sect. 4.1), or getting the ambipolar diffusion coefficient on the fly, which is mandatory for NEQ ionization and recombination calculations (Sect. 4.2).

4.1. LTE

For LTE numerical experiments, Bifrost has an EOS based on tables generated by the Uppsala Opacity Package (Gustafsson et al. 1975), which computes instantaneous molecular dissociation equilibria and the LTE ionization balance of the 16 most important elements: elements with a high abundance (greater than 7 on a logarithmic scale) complemented with elements with an abundance down to 5 but with low ionization potential, and therefore important as electron donors at low temperatures (see Gudiksen et al. 2011, for details). Table 1 shows the list of those 16 elements together with the abundances, A, and atomic mass. Abundances were kept at the original values of Gustafsson et al. (1975) to ensure compatibility with simulations by Stein & Nordlund (2000).

Table 1.

Abundances (A) (Gustafsson et al. 1975) and atomic mass (Z) for the 16 elements (El.) used in the EOS of Bifrost.

In order to tabulate the ambipolar diffusion coefficient, we create a table for the η amb ∗ $ {\eta^{\ast}_{{\rm amb}}} $ coefficient, defined through Eq. (3), which does not depend on the magnetic field, and can therefore be given as a function of the density, ρ, and internal energy per unit volume, ϵ. When executing numerical experiments, interpolations are carried out to get the ambipolar diffusion coefficient at the input energy and density. Figure 2 shows η amb ∗ $ {\eta^{\ast}_{{\rm amb}}} $ as a function of ρ and ϵ from the EOS table with temperature isocontours superimposed. The zero-point for the excitation/ionization/dissociation energy part of the internal energy is arbitrary. We chose this energy to be 5 eV for a completely neutral gas with no molecules. This ensures a positive internal energy even when all hydrogen is bound in molecular form.

thumbnail Fig. 2.

Tabulated η amb ∗ $ {\eta^{\ast}_{{\rm amb}}} $ as a function of density, ρ, and internal energy per unit volume, ϵ. The values for ϵ are shifted by 5 eV to obtain a positive internal energy when having hydrogen molecules (see main text). Temperature isocontours are superimposed as dashed lines for (from left to right) 3 × 103 K, 6 × 103 K, 104 K, 2 × 104 K, and 105 K.

We find that the ambipolar diffusion coefficient has a strong dependence on the number of elements considered in the EOS. In order to show this important fact, in Fig. 3 we compare η amb ∗ $ {\eta^{\ast}_{{\rm amb}}} $, as a function of ρ and T, obtained using modified versions of the EOS of Bifrost. The figure shows the complete EOS of Bifrost (panel A); an EOS only considering atoms of H, He, Na, Si, Mg, K, and molecules of H2 (panel B); an EOS just taking H, H2, and He into account (panel C); and an EOS only including atomic H (panel D). Comparing panels A and D, we see that considering a pure hydrogen plasma leads to an overestimation of the role of the ambipolar diffusion by several orders of magnitude in the coolest temperatures. The reason is that as the temperature decreases and hydrogen becomes neutral, the total number of ions drops more quickly than when heavier elements are present, resulting in larger values of the ambipolar diffusion coefficient. Assuming a pure hydrogen plasma also underestimates the ambipolar diffusion above log(T)∼3.8, in this case because as the temperature increases, hydrogen becomes ionized leading to a lack of neutrals. Looking at panels A and C now, it is clear that introducing helium helps obtain the correct values of η amb ∗ $ {\eta^{\ast}_{{\rm amb}}} $ above log(T)∼3.8. Introducing molecules of H2 slightly modifies the behavior in the range of temperature below log(T)∼3.4 (see the changes in the isocontours of the color scale); however, the values of η amb ∗ $ {\eta^{\ast}_{{\rm amb}}} $ are still overestimated in this regime due to the lack of the main electron donors. Panel B shows that we can start to get the right order of magnitude of the ambipolar diffusion coefficient once we introduce some of the main electron donors in the chromosphere, such as sodium, silicon, magnesium, and potassium.

thumbnail Fig. 3.

Tabulated η amb ∗ $ {\eta^{\ast}_{{\rm amb}}} $ as a function of the density, ρ, and temperature, T, for different EOS. Panel A: Complete EOS of Bifrost. Panel B: EOS only considering atoms of H, He, Na, Si, Mg, K, and molecules of H2. Panel C: EOS taking only H, H2, and He into account. Panel D: EOS including only atomic H.

4.2. NEQ

In the previous section we show that the inclusion of elements heavier than hydrogen hugely impacts on the ambipolar diffusion coefficient calculation. This is illustrated through different LTE EOS; nonetheless, in the solar atmosphere, there are important departures from ionization equilibrium (see the seminal papers by Klein et al. 1976, 1978; Kneer 1980; Carlsson & Stein 1992, 2002), mainly in the chromosphere and transition region (e.g., Bradshaw & Mason 2003; Bradshaw & Cargill 2006; Olluri et al. 2013, 2015; Martínez-Sykora et al. 2016; Nóbrega-Siverio et al. 2018; Kerr et al. 2019; Rutten 2019). This means that it is not only important to include all the relevant elements, but also to obtain the ion and neutral number density for them through a nonequilibrium ionization calculation. Moreover, it is crucial to compute the ambipolar diffusion coefficient due to its strong dependency on the number density of ions and neutrals. In the following, the cases of nonequilibrium ionization of hydrogen and helium are commented separately.

4.2.1. NEQ of hydrogen

The Bifrost code has a module available that calculates the ionization and recombination of atomic hydrogen and the formation and dissociation of molecular hydrogen, H2, under NEQ conditions (see Leenaarts et al. 2007, 2011, for details about this module). The ambipolar diffusion module presented in this paper can take the atomic and molecular hydrogen number densities from this NEQ module to compute the collision frequency (Eq. (11)) and the ambipolar diffusion coefficient (Eq. (3)), while the ionization state of the rest of the elements is assumed in LTE. Using the NEQ of hydrogen module together with the new ambipolar diffusion module presented here, Nóbrega-Siverio et al. (2020) recently showed that the LTE assumption can lead to an important underestimation of the ionization fraction in magnetic flux emerging regions of up to 2–3 orders of magnitude. This means that flux emergence experiments carried out under the LTE assumption significantly overestimate the role of the ambipolar diffusion.

4.2.2. NEQ of hydrogen and helium

The Bifrost code also has the capability of computing the NEQ ionization and recombination of helium by means of a module developed by Golding et al. (2016). The NEQ ionization and recombination of helium is important in the upper chromosphere, where the helium ion fractions considerably depart from their equilibrium values. This module always works together with the NEQ module of hydrogen, so when running ambipolar diffusion experiments, we take the atomic and molecular hydrogen, and the atomic helium number densities from this module to calculate νni and ηamb, while the rest of elements are assumed in LTE. The combination of these modules largely impacts the ambipolar diffusion distribution and changes the ambipolar heating in comparison to computations under the LTE assumption (Martínez-Sykora et al. 2020a), which can provide a completely new interpretation of the observations in the chromosphere (Martínez-Sykora et al. 2020b).

5. Super time stepping (STS) method

Super time stepping (STS; Alexiades et al. 1996) is a technique that can be used to accelerate explicit parabolic calculations. It is based on the stability properties of the Chebyshev polynomials, which allow us to relax the CFL condition and take a larger timestep, ΔtSTS. The method uses two input parameters, n, a positive integer, and ν ∈ (0, 1), which is a damping factor. The timestep allowed by the STS method is then given by

Δ t STS = n 2 ν [ ( 1 + ν ) 2 n − ( 1 − ν ) 2 n ( 1 + ν ) 2 n + ( 1 − ν ) 2 n ] Δ t AD , CFL , $$ \begin{aligned} \Delta t_{\mathrm{STS}} = \frac{n}{2 \sqrt{\nu }} \left[ \frac{ (1 + \sqrt{\nu })^{2n} - (1 - \sqrt{\nu })^{2n}}{(1 + \sqrt{\nu })^{2n} + (1 - \sqrt{\nu })^{2n}} \right] \Delta t_{\mathrm{AD,CFL}} , \end{aligned} $$(13)

and is divided into n sub-timesteps τi as

Δ t STS = ∑ i = 1 n τ i , $$ \begin{aligned} \Delta t_{\mathrm{STS}} = \sum _{i=1}^{n} \tau _i ,\end{aligned} $$(14)

with

τ i = Δ t AD , CFL [ ( ν − 1 ) cos ( π ( 2 i − 1 ) 2 n ) + ( ν + 1 ) ] − 1 , $$ \begin{aligned} \tau _i = \Delta t_{\mathrm{AD,CFL}} \left[ (\nu - 1) \cos {\left( \frac{\pi (2i - 1)}{2n} \right)} + (\nu + 1) \right]^{-1}\;, \end{aligned} $$(15)

where ΔtAD,CFL is the timestep of the parabolic problem (in our case, ambipolar diffusion) given by the classic CFL condition. It is important to note that the intermediate values computed along the n sub-timesteps have no approximating properties: it is only after the whole ΔtSTS has been reached that the results approximate the solution and, consequently, have a physical meaning.

The above expressions can be easily implemented to improve the calculation efficiency of ambipolar diffusion experiments; nevertheless, the drawback is that the STS is a first-order scheme in time. In addition, the method has two free parameters, n and ν, and it is necessary to carefully choose their values to optimize the performance while keeping numerical stability. The maximum ratio ΔtSTS/ΔtCFL that can be reached for any given n corresponds to the following limit:

lim ν → 0 Δ t STS Δ t AD , CFL = lim ν → 0 n 2 ν [ ( 1 + ν ) 2 n − ( 1 − ν ) 2 n ( 1 + ν ) 2 n + ( 1 − ν ) 2 n ] = n 2 $$ \begin{aligned} \lim _{\nu \rightarrow 0} \frac{\Delta t_{\mathrm{STS}} }{\Delta t_{\mathrm{AD,CFL}} } = \lim _{\nu \rightarrow 0} \frac{n}{2 \sqrt{\nu }} \left[ \frac{ (1 + \sqrt{\nu })^{2n} - (1 - \sqrt{\nu })^{2n}}{(1 + \sqrt{\nu })^{2n} + (1 - \sqrt{\nu })^{2n}} \right] = n^2 \end{aligned} $$(16)

In this limit, the CFL criterion would need n2 steps to reach one ΔtSTS, while the STS method requires just n, as explained above (Eq. (15)). Assuming a similar computing load for each step, whether in the STS or in the CFL-limited calculation, this implies that the STS method would be n times faster than the simple CFL-limited one. However, it is necessary to impose a lower threshold for ν because ν = 0 is a stability limit, and choosing low values of ν can make the STS method very sensitive to round-off errors (Alexiades et al. 1996). In the next subsection, we explain our choices made for ν and n.

5.1. Implementation of the STS method in the Bifrost code

In this section we show the STS implementation in the Bifrost code and the choice of the n and ν parameters. To this end, we use an operator splitting scheme, which is a standard technique to advance in time the solution at each timestep separately for groups of terms that contain different physics, hence with different numerical requirements, while keeping codes modular. In our case, the induction and energy equations (Eqs. (5) and (10), respectively) are separated and solved over a timestep Δt determined by the minimum of the two following values: (1) the standard timestep given by the MHD and radiation terms considered in the Bifrost code (see Gudiksen et al. 2011, for further details); and (2) the timestep imposed by the ambipolar and Hall term that appear in the generalized Ohm’s law from Eq. (1). In particular, we take advantage of the operator splitting to solve separately the ambipolar diffusion term using, when necessary, the STS method. In order to determine whether the STS method can be efficiently applied, let us assume a system in which

Δ t MHD , CFL = C Δ t AD , CFL , C ∈ R + . $$ \begin{aligned} \Delta t_{\mathrm{MHD,CFL}} = C\; \Delta t_{\mathrm{AD,CFL}} , \quad C \in \mathbb{R} ^+. \end{aligned} $$(17)

In this equation ΔtMHD,CFL is the timestep given by the CFL criterion for the standard MHD part, and ΔtAD,CFL is the timestep for the ambipolar diffusion term. For any acceleration method like the STS to be of interest, it is obviously necessary that ΔtMHD,CFL >  ΔtAD,CFL (i.e., C >  1).

5.2. Choice of the free parameters n and ν

In the literature, even though STS has been broadly implemented in various codes and physical processes (see Sect. 1), there is no thorough analysis of the choice of the n and ν parameters, so the use of the STS method is not precisely specified. We have investigated several properties of these STS free parameters. In the following, we detail various aspects that should be considered:

1. ΔtSTS should not be larger than the timestep imposed by the general MHD and radiation terms (ΔtMHD,CFL) since the time advance of the system is limited by the minimum of all applicable timestep conditions. Consequently,

Δ t AD , CFL ≤ Δ t STS ≤ Δ t MHD , CFL . $$ \begin{aligned} \Delta t_{\mathrm{AD,CFL}} \le \Delta t_{\mathrm{STS}} \le \Delta t_{\mathrm{MHD,CFL}} . \end{aligned} $$(18)

This condition can be rewritten using Eqs. (13) and (17) as

1 ≤ n 2 ν [ ( 1 + ν ) 2 n − ( 1 − ν ) 2 n ( 1 + ν ) 2 n + ( 1 − ν ) 2 n ] ≤ C $$ \begin{aligned} 1 \le \frac{n}{2 \sqrt{\nu }} \left[ \frac{ (1 + \sqrt{\nu })^{2n} - (1 - \sqrt{\nu })^{2n}}{(1 + \sqrt{\nu })^{2n} + (1 - \sqrt{\nu })^{2n}} \right] \le C \end{aligned} $$(19)

with ν ∈ (0, 1) and C >  1. The above is the minimum requirement to apply the STS method. For n = 1 the function at the center of the inequality is below 1 for all ν >  0; therefore, since n is an integer, Eq. (19) implies n ≥ 2.

2. In addition, we tighten the lower bound of Eqs. (18) and (19) as follows: for the STS method to be computationally advantageous, its n substeps must permit a greater advance in time than n times the timestep ΔtAD,CFL. Thus, in Eq. (18) we set nΔtAD,CFL <  ΔtSTS, so the lower bound of Eq. (19) becomes

n < n 2 ν [ ( 1 + ν ) 2 n − ( 1 − ν ) 2 n ( 1 + ν ) 2 n + ( 1 − ν ) 2 n ] ≡ f ( n , ν ) $$ \begin{aligned} n < \frac{n}{2 \sqrt{\nu }} \left[ \frac{ (1 + \sqrt{\nu })^{2n} - (1 - \sqrt{\nu })^{2n}}{(1 + \sqrt{\nu })^{2n} + (1 - \sqrt{\nu })^{2n}} \right] \equiv f(n,\nu ) \end{aligned} $$(20)

For simplicity, we have introduced the symbol f(n, ν) for the term on the right-hand side of the expression. It can be shown that the constraint (20) imposes a maximum value for ν. As proof, in Fig. 4 we have plotted f(n, ν) against n for different values of ν (dashed lines). In the image the solid red curve is f(n, ν) = n2, which corresponds to the limit of stability (ν = 0). The solid black curve is simply the straight f(n, ν) = n that allows us to identify which combinations of n and ν verify the above inequality. As seen in the figure, f(n, ν) = n roughly coincides with the curve ν = 0.25. Through a limit analysis of f(n, ν)/n,

lim n → ∞ 1 2 ν [ ( 1 + ν ) 2 n − ( 1 − ν ) 2 n ( 1 + ν ) 2 n + ( 1 − ν ) 2 n ] = 1 2 ν , $$ \begin{aligned} \lim _{n \rightarrow \infty } \frac{1}{2 \sqrt{\nu }} \left[ \frac{ (1 + \sqrt{\nu })^{2n} - (1 - \sqrt{\nu })^{2n}}{(1 + \sqrt{\nu })^{2n} + (1 - \sqrt{\nu })^{2n}} \right] = \frac{1}{2\sqrt{\nu }}, \end{aligned} $$(21)

thumbnail Fig. 4.

Plot showing f(n, ν) vs. n for different values of ν (dashed curves). Solid lines represent f(n, ν) = n2 (red) and f(n, ν) = n (black).

we find that the constraint f(n, ν)/n = 1 gives us an asymptotic limit of ν = 0.25, the maximum value of ν that results from the optimization constraint (20). In the code, to avoid working with asymptotic limits, we use as the upper limit of ν the value obtained for the lowest possible n ( i . e . , n = 2 ) : ν < 5 − 2 $ (i.e., n=2): \nu < \sqrt{5} - 2 $.

3. There is also a strong constraint concerning the use of large n. In principle it is mathematically possible to have any value of n; nonetheless, numerically, large values of n can be problematic. This is true because, during the n sub-cycling shown in Eq. (15), the inner values obtained within the loop do not have any physical meaning, as already mentioned, and can reach very large values that compromise the accuracy of the calculations due to the limited precision. There are three ways to alleviate this problem: (a) increasing from single to double precision, which allows calculations with larger values of n, but increases the memory requirements; (b) introducing an elaborate normalization, which helps to avoid Infinity or !NAN values, but increases the number of operations and complicates its implementation; (c) imposing a strict maximum of n according to the experience. The last is the solution we adopted; it requires some testing to determine the valid n that will keep the simulations stable. We find that n starts to be problematic above approximately 12. In order to ensure greatest stability, we decided to impose as a maximum n = 10. This, together with the previous considerations, gives us the two following ranges we must fulfill: 0 < ν < 5 − 2 $ 0 < \nu < \sqrt{5} - 2 $ and 2 ≤ n ≤ 10. In Fig. 5, the valid parameter-space derived from these ranges is shown within the three solid red lines of the image.

thumbnail Fig. 5.

Plot showing f(n, ν) versus n for different values of ν (dashed curves). Solid lines represent the boundaries of the parameter-space considering the constraints described in the text related to optimization and stability, namely f(n, ν) = n2, f(n = 10, ν) and, f ( n , ν = 5 − 2 ) $ f(n,\nu=\sqrt{5}-2) $.

4. The final step to determine the optimum combination of n and ν can now be taken. The ideal situation would be to reach ΔtMHD, CFL in the minimum number of stable n steps without compromising the accuracy of the solution. To this end, we impose

Δ t STS = Δ t MHD , CFL , $$ \begin{aligned} \Delta t_{\mathrm{STS}} = \Delta t_{\mathrm{MHD,CFL}} ,\end{aligned} $$(22)

which implies the following relationship between n and ν:

f ( n , ν ) = n 2 ν [ ( 1 + ν ) 2 n − ( 1 − ν ) 2 n ( 1 + ν ) 2 n + ( 1 − ν ) 2 n ] = C $$ \begin{aligned} f(n,\nu ) = \frac{n}{2 \sqrt{\nu }} \left[ \frac{ (1 + \sqrt{\nu })^{2n} - (1 - \sqrt{\nu })^{2n}}{(1 + \sqrt{\nu })^{2n} + (1 - \sqrt{\nu })^{2n}} \right] = C \end{aligned} $$(23)

We could obtain the minimum number of substeps n in the limit ν → 0, so, from Eq. (16) and considering that n has to be an integer, n = int ( C ) $ n=\mathrm{int}(\sqrt{C}) $. However, since ν = 0 corresponds to the limit of stability, we need to consider n > int ( C ) $ n > \mathrm{int}(\sqrt{C}) $. The first option is to use n = C + 1 $ n = \sqrt{C} + 1 $; nevertheless, the associated ν could still be small to affect our calculations by round-off errors as pointed out by Alexiades et al. (1996). For this reason, in our STS calculations we impose

n = int ( C ) + 2 with C > 4 $$ \begin{aligned} n = \mathrm{int} (\sqrt{C}) + 2 \quad \text{ with}\quad C > 4 \end{aligned} $$(24)

For the cases when 2 ≤ C ≤ 4, instead of STS we apply the sub-cycling method described by Martínez-Sykora et al. (2017b), that is, only the induction and energy equation are evolved separately from the rest of MHD equations; when 1 ≤ C ≤ 2 we follow the CFL criterion.

Having determined n, to apply the STS method, we can see that for each n there is a limited range of values of ν that satisfies Eq. (23) for C >  4. This is illustrated in Fig. 6 through a 2D map of C/n showing the range of ν that satisfy that condition for each n as a colored patch. The corresponding values of C/n are given in the color bar. In the image, we overplot a red horizontal line to delimit ν = 5 − 2 $ \nu = \sqrt{5} - 2 $, and a black one that delineates the minimum value of ν for the maximum n = 10 (ν ≈ 0.00185). We see that with this implementation, we can speed up the calculations by a factor 8. Using the criteria just indicated for the choice of the ν and n parameters, we can get the best performance, optimizing the speed for the ambipolar diffusion calculations in the Bifrost code.

thumbnail Fig. 6.

Two-dimensional map of C/n showing the combinations of n and ν that verify Eq. (23) for C >  4. In the image the horizontal dashed lines indicate the values ν = 5 − 2 $ \nu = \sqrt{5} - 2 $ (red) and the value of ν for the maximum n = 10 (ν ≈ 0.00185, black).

6. Hyperdiffusion

In order to maintain the stability of the code, we have implemented hyperdiffusion terms related to the ambipolar diffusion. The diffusive operator for the induction equation follows the idea originally described by Galsgaard & Nordlund (1995) and later implemented in the Bifrost code (Gudiksen et al. 2011) and in its preexisting non-STS generalized Ohm’s law module (Martínez-Sykora et al. 2017b). In the following, for brevity, we only describe the hyperdiffusion operator for the induction equation in the x direction. Equivalent explanations can be given for the other two directions.

∂ B x ∂ t = … + ∂ ∂ z [ Δ x 2 ( ν x ( 1 ) + ν x ( 2 ) + ν x ( 3 ) ) Q ( ∂ B z ∂ x ) J y ] + ∂ ∂ z [ Δ z 2 ( ν z ( 1 ) + ν z ( 2 ) + ν z ( 3 ) ) Q ( ∂ B x ∂ z ) J y ] − ∂ ∂ y [ Δ x 2 ( ν x ( 1 ) + ν x ( 2 ) + ν x ( 3 ) ) Q ( ∂ B y ∂ x ) J z ] − ∂ ∂ y [ Δ y 2 ( ν y ( 1 ) + ν y ( 2 ) + ν y ( 3 ) ) Q ( ∂ B x ∂ y ) J z ] , $$ \begin{aligned} \begin{split} \frac{\partial B_x}{\partial t} = \ldots&+ \frac{\partial }{\partial z}\left[ \frac{\Delta x}{2} \left(\nu ^{(1)}_x + \nu ^{(2)}_x + \nu ^{(3)}_x \right)Q\left(\frac{\partial B_z}{\partial x} \right)J_{ y} \right] \\&+ \frac{\partial }{\partial z}\left[ \frac{\Delta z}{2} \left(\nu ^{(1)}_z + \nu ^{(2)}_z + \nu ^{(3)}_z \right)Q\left(\frac{\partial B_x}{\partial z} \right)J_{ y} \right] \\&- \frac{\partial }{\partial { y}}\left[ \frac{\Delta x}{2} \left(\nu ^{(1)}_x + \nu ^{(2)}_x + \nu ^{(3)}_x \right)Q\left(\frac{\partial B_{ y}}{\partial x} \right)J_z \right] \\&- \frac{\partial }{\partial { y}}\left[ \frac{\Delta y}{2} \left(\nu ^{(1)}_{ y} + \nu ^{(2)}_{ y} + \nu ^{(3)}_{ y} \right)Q\left(\frac{\partial B_x}{\partial { y}} \right)J_z \right] \end{split},\end{aligned} $$(25)

where Δx, Δy, and Δz are the grid spacing in the x, y, and z direction, respectively; Q is a positive definite quenching factor (see Galsgaard & Nordlund 1995; Gudiksen et al. 2011; Martínez-Sykora et al. 2017b); and ν j (1) $ \nu^{(1)}_j $, ν j (2) $ \nu^{(2)}_j $, and ν j (3) $ \nu^{(3)}_j $ are the three hyperdiffusion terms in the x, y, or z direction that we describe in the following.

– ν j (1) $ \nu^{(1)}_j $ is a hyperdiffusion term proportional to the ambipolar diffusion coefficient that damps waves with the shortest wavelengths, i.e., the largest wavenumber: the Nyquist wavenumber. The term is

ν j ( 1 ) = ν 1 η amb π Δ j , $$ \begin{aligned} \nu ^{(1)}_j = \nu _1 \eta _{\mathrm{amb}}\frac{\pi }{\Delta j},\end{aligned} $$(26)

where ν1 is a dimensionless coefficient of order 10−2.

– ν j (2) $ \nu^{(2)}_j $ is a hyperdiffusion term related to the advective transport of magnetic field attached to the charged species, or more precisely to the electrons. It is defined as

ν j ( 2 ) = ν 2 | u amb , j | , $$ \begin{aligned} \nu ^{(2)}_j = \nu _2 | u_{{\mathrm{amb} }, j} | ,\end{aligned} $$(27)

where ν2 is a dimensionless coefficient of order 10−2. The Bifrost code already has a hyperdiffusion term because of the advection of the magnetic field with the velocity u, and for that reason here we only need to consider the extra part due to ambipolar diffusion.

– ν j (3) $ \nu^{(3)}_j $ is a hyperdiffusion term related to magnetic shocks. The magnetic field is only influenced by the perpendicular velocity field, so it is possible to adopt a hyperdiffusion term that depends on the convergence of the velocity field. As in the previous term, Bifrost has already implemented a term for u⊥, so here we only have to implement one for the ambipolar diffusion velocity. Since this velocity field is already perpendicular, the ν j (3) $ \nu^{(3)}_j $ term is

ν j ( 3 ) = ν 3 Δ j | ∇ · u amb | j , $$ \begin{aligned} \nu ^{(3)}_j = \nu _3 \Delta j | \nabla \cdot {{\boldsymbol{u}}_{{\mathrm{amb} }}} |_j,\end{aligned} $$(28)

where ∇ is a first-order finite difference operator and ν3 a dimensionless coefficient of order 10−2.

7. Validation test

After programming the STS following the parameter choice described in Sect. 5 and explaining the hyperdiffusion terms in Sect. 6, we need to verify that the STS is correctly implemented in the Bifrost code. To that end, we performed a test that has the advantage of allowing us to verify the method in 2.5D MHD scenarios for different initial conditions. This is possible thanks to the existence of a self-similar analytical solution (see Pattle 1959). In the following, we describe this analytical solution and the goodness of our STS implementation.

7.1. Self-similar solution for ambipolar diffusion problems

The test we intend to carry out is based on the induction equation when only taking the ambipolar term into account, in other words the pure ambipolar diffusion case. If we also assume that the ambipolar coefficient η amb ∗ $ {\eta^{\ast}_{{\rm amb}}} $ is homogeneous, the induction equation (Eq. (5)) becomes

∂ B ∂ t = η amb ∗ ∇ × [ ( J × B ) × B ] . $$ \begin{aligned} \frac{\partial \boldsymbol{B}}{\partial t} = \eta ^{*}_{\mathrm{amb}}\nabla \times \left[ (\boldsymbol{J \times B}) \times \boldsymbol{B} \right]. \end{aligned} $$(29)

In a 2D domain with polar coordinates, this equation becomes a nonlinear diffusion equation when considering an axisymmetric function for B. Pattle’s self-similar solution is of the form

B y ( r , t ) ⟶ ( Φ 4 π η amb ∗ t ) 1 / 3 [ 1 − r 2 R 2 ( t ) ] 1 / 2 $$ \begin{aligned} B_{ y} (r, t) \longrightarrow \left(\frac{\Phi }{4 \pi \eta ^{*}_{\mathrm{amb}}t} \right)^{1/3} \left[ 1 - \frac{r^2}{R^2 (t)} \right]^{1/2} \end{aligned} $$(30)

for r <  R(t), and 0 otherwise, where

R 2 ( t ) ⟶ 3 2 [ 4 η amb ∗ t ( Φ π ) 2 ] 1 / 3 . $$ \begin{aligned} R^2 (t) \longrightarrow \frac{3}{2} \left[ 4\eta ^{*}_{\mathrm{amb}}t \left(\frac{\Phi }{\pi } \right)^2 \right]^{1/3} .\end{aligned} $$(31)

Here Φ is the magnetic flux, y is the symmetry axis, and r2 = (x − x0)2 + (z − z0)2. From these expressions it is clear that By ∝ t−1/3 and R ∝ t1/6. The problem has the additional feature that any single-lobed initial condition with total flux Φ is expected to tend, asymptotically in time, to the shape given in Eqs. (30) and (31) (Zel’dovich & Raizer 1967). This is a robust multidimensional test to check whether the STS method was properly implemented in the code.

7.2. Numerical test

In order to test the STS method with that self-similar solution, we create a 2.5D snapshot of 1024 × 1024 points. The physical domain spans from −5.0 ≤ x ≤ 7.5 Mm to −5.0 ≤ z ≤ 7.5 Mm, with a numerical resolution of Δx = Δz = 12 km. In that domain, we define as initial condition the following axisymmetric function for the magnetic field:

B y 0 ( x , z ) = B 0 e − [ ( x − x 0 ) 2 + ( z − z 0 ) 2 ] / w 2 , $$ \begin{aligned} B_{{ y}0} (x,z) = B_0 e^{ - \left[ (x-x_0)^2 + (z-z_0)^2 \right]/w^2} ,\end{aligned} $$(32)

where B0 = 0.8 G, x0 = 0 Mm, z0 = 0 Mm, and w = 0.5 Mm. The rest of parameters are representative of the chromosphere: the initial internal energy per unit volume is e0 = 4 × 10−3 erg cm−3; the initial density ρ0 = 10−9 g cm−3; and the ambipolar diffusion coefficient is set to a constant, namely η amb ∗ = 10 16 $ {\eta^{\ast}_{{\rm amb}}}= 10^{16} $ cm2 s−1 G−2. This initial condition does indeed converge to the self-similar solution with the same integrated magnetic flux after an initial transient phase. The results for the self-similar stage are shown in Fig. 7. In the image the top panel illustrates the evolution of the maximum of the magnetic field with time through black asterisks. The power-law fit to that curve is

B y = 0.100 t − 0.332 , $$ \begin{aligned} B_{ y} = 0.100 t^{-0.332}, \end{aligned} $$(33)

thumbnail Fig. 7.

Results of the STS validation test for the self-similar solution. Top panel: evolution of the maximum of the magnetic field with time (black asterisks) and power-law fit (red line). Bottom panel: evolution of Xb with time (black asterisks) and power-law fit (red line). Above each panel is shown the resulting fit formula.

and it is shown as a red curve in the figure. The exponent in (33) agrees quite well with the predicted value, −1/3, shown in Eq. (30). On the other hand, the bottom panel contains the evolution of Xb (one of the Cartesian components of R) with time (black asterisks) and its power-law fit (red). The power-law fit yields

X b = − 1.732 t 0.166 , $$ \begin{aligned} X_{b} = -1.732 t^{0.166}, \end{aligned} $$(34)

which again agrees with the expected value 1/6 from Eq. (31). This tests indicates that our STS implementation is correct. In addition, to ascertain the speed-up factor that we gain by means of the STS, we ran the same test with the same number of CPUs under the CFL criterion and also using the sub-cycling method described by Martínez-Sykora et al. (2017b). The results show that for this test, STS can lead to a speed-up factor of 8 and 4 with respect to the CFL and sub-cycling method, respectively.

8. Conclusions

In this paper, we described the implementation of a new module in the Bifrost code to efficiently address ambipolar diffusion problems in the solar atmosphere with the state-of-the-art microphysics, pursuing a fourfold purpose:

– We reviewed the existing literature to implement the most accurate collision cross sections and frequencies, which are crucial for proper calculation of the ambipolar diffusion coefficient. For instance, we considered the temperature-dependent cross section of the collisions between hydrogen and protons (Vranjes & Krstic 2013), which can be roughly one order of magnitude larger than the constant value usually taken in the literature (e.g., Khodachenko et al. 2004; Leake & Arber 2006; Soler et al. 2009) for temperatures below 104 K. In addition, we included the cross sections of collisions with hydrogen molecules that may be important in the coolest regions of the Sun, such as sunspot umbrae, where estimations from observations reveal that there is a substantial H2 molecule formation present (Jaeggli et al. 2012).

– Through an analysis of the elements included in the EOS, we showed that in the solar atmosphere it is necessary to include heavier elements than hydrogen, primarily helium and the most important electron donors, to properly estimate the role of the ambipolar diffusion. Considering a pure hydrogen plasma can lead to an overestimation of the ambipolar diffusion coefficient by several orders of magnitude.

– Since the correct determination of the ionization fraction is also relevant for the ambipolar diffusion calculation, we also made this new module compatible with previous Bifrost modules that compute the NEQ ionization and recombination of hydrogen and helium as well as the H2 formation under NEQ conditions (Leenaarts et al. 2011; Golding et al. 2016). This way, we opened a new avenue to address partial ionization effects together with departures of the ionization state from the LTE in the solar atmosphere. In fact, NEQ ionization and recombination effects of hydrogen are shown to be key for the understanding of the ambipolar diffusion in magnetic flux emergence processes in the Sun (Nóbrega-Siverio et al. 2020). In addition, Martínez-Sykora et al. (2020a), combining the module presented in this paper together with the two above-mentioned NEQ modules of hydrogen and helium, have shown the important variations in the thermodynamics in the chromosphere and transition region with respect to calculations carried out assuming LTE.

– We also implemented the STS method in the Bifrost code, carrying out a thorough analysis of the choice of the two free parameters of the STS method to obtain the best performance. This choice allows us to speed up single-precision computations including the ambipolar diffusion term up to a factor 8 in comparison with the same calculation run under the CFL criterion. In addition, we described the hyperdiffusion terms that maintain the stability of the code when ambipolar diffusion is taken into account.

The numerical techniques presented in this paper may give hints to address ambipolar diffusion problems in other astrophysical systems.

1

For higher values of the energy ECM, i.e., 102 to 106 eV, see Schultz et al. (2008).

Acknowledgments

This research is supported by the Research Council of Norway through its Centres of Excellence scheme, project number 262622, and through grants of computing time from the Programme for Supercomputing. It is also supported by the Spanish Ministry of Science, Innovation and Universities through projects AYA2014-55078-P and PGC2018-095832-B-I00, as well as through the Synergy Grant number 810218 (ERC-2018-SyG) of the European Research Council; by NASA through grants NNX17AD33G, 80NSSC18K1285 and contract NNG09FA40C (IRIS) and by the NSF grant AST1714955. The authors thankfully acknowledge the computer resources provided at the Pleiades cluster through the computing projects s1061, and s2053 from the High End Computing (HEC) division of NASA. In addition, this study has been discussed within the activities of team 399 “Studying magnetic-field-regulated heating in the solar chromosphere” at the International Space Science Institute (ISSI) in Switzerland.

References

  1. Alexiades, V., Amiez, G., & Gremaud, P.-A. 1996, Commun. Numer. Methods Eng., 12, 31 [CrossRef] [Google Scholar]
  2. Alvarez Laguna, A., Lani, A., Deconinck, H., Mansour, N. N., & Poedts, S. 2016, J. Comput. Phys., 318, 252 [NASA ADS] [CrossRef] [Google Scholar]
  3. Ballester, J. L., Alexeev, I., Collados, M., et al. 2018, Space Sci. Rev., 214, 58 [NASA ADS] [CrossRef] [Google Scholar]
  4. Basu, S., & Ciolek, G. E. 2004, ApJ, 607, L39 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bradshaw, S. J., & Cargill, P. J. 2006, A&A, 458, 987 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Bradshaw, S. J., & Mason, H. E. 2003, A&A, 401, 699 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Braginskii, S. I. 1965, Rev. Plasma Phys., 1, 205 [Google Scholar]
  8. Carlsson, M., & Stein, R. F. 1992, ApJ, 397, L59 [NASA ADS] [CrossRef] [Google Scholar]
  9. Carlsson, M., & Stein, R. F. 2002, ApJ, 572, 626 [NASA ADS] [CrossRef] [Google Scholar]
  10. Cheung, M. C. M., & Cameron, R. H. 2012, ApJ, 750, 6 [NASA ADS] [CrossRef] [Google Scholar]
  11. Choi, E., Kim, J., & Wiita, P. J. 2009, ApJS, 181, 413 [NASA ADS] [CrossRef] [Google Scholar]
  12. Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Courant, R., Friedrichs, K., & Lewy, H. 1928, Math. Ann., 100, 32 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  14. Cowling, T. G. 1976, Magnetohydrodynamics (Crane Russak & Co) [Google Scholar]
  15. Crutcher, R. M. 2012, ARA&A, 50, 29 [NASA ADS] [CrossRef] [Google Scholar]
  16. Franz, G. 2009, Low Pressure Plasmas and Microstructuring (Technology (Springer)) [CrossRef] [Google Scholar]
  17. Galsgaard, K., & Nordlund, Å. 1995, A 3D MHD code (for Parallel Computers) [Google Scholar]
  18. Glassgold, A. E., Krstić, P. S., & Schultz, D. R. 2005, ApJ, 621, 808 [NASA ADS] [CrossRef] [Google Scholar]
  19. Golding, T. P., Leenaarts, J., & Carlsson, M. 2016, ApJ, 817, 125 [NASA ADS] [CrossRef] [Google Scholar]
  20. González-Morales, P. A., Khomenko, E., Downes, T. P., & de Vicente, A. 2018, A&A, 615, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Goodman, M. L. 2004, A&A, 416, 1159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Grassi, T., Padovani, M., Ramsey, J. P., et al. 2019, MNRAS, 484, 161 [CrossRef] [Google Scholar]
  23. Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84 [NASA ADS] [CrossRef] [Google Scholar]
  24. Gudiksen, B. V., Carlsson, M., Hansteen, V. H., et al. 2011, A&A, 531, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Gustafsson, B., Bell, R. A., Eriksson, K., & Nordlund, A. 1975, A&A, 42, 407 [NASA ADS] [Google Scholar]
  26. Huang, Y.-M., Bhattacharjee, A., & Sullivan, B. P. 2011, Phys. Plasmas, 18, 072109 [NASA ADS] [CrossRef] [Google Scholar]
  27. Iijima, H., & Yokoyama, T. 2015, ApJ, 812, L30 [NASA ADS] [CrossRef] [Google Scholar]
  28. Jaeggli, S. A., Lin, H., & Uitenbroek, H. 2012, ApJ, 745, 133 [NASA ADS] [CrossRef] [Google Scholar]
  29. Kerr, G. S., Carlsson, M., Allred, J. C., Young, P. R., & Daw, A. N. 2019, ApJ, 871, 23 [CrossRef] [Google Scholar]
  30. Khodachenko, M. L., Arber, T. D., Rucker, H. O., & Hanslmeier, A. 2004, A&A, 422, 1073 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Khomenko, E., & Collados, M. 2012, ApJ, 747, 87 [NASA ADS] [CrossRef] [Google Scholar]
  32. Khomenko, E., Díaz, A., de Vicente, A., Collados, M., & Luna, M. 2014, A&A, 565, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Klein, R. I., Stein, R. F., & Kalkofen, W. 1976, ApJ, 205, 499 [NASA ADS] [CrossRef] [Google Scholar]
  34. Klein, R. I., Stein, R. F., & Kalkofen, W. 1978, ApJ, 220, 1024 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Kneer, F. 1980, A&A, 87, 229 [NASA ADS] [Google Scholar]
  36. Krstić, P. S., & Schultz, D. R. 1999a, Phys. Rev. B, 60, 2118 [NASA ADS] [CrossRef] [Google Scholar]
  37. Krstic, P. S., & Schultz, D. R. 1999b, J. Phys. B At. Mol. Phys., 32, 2415 [NASA ADS] [CrossRef] [Google Scholar]
  38. Kudoh, T., & Basu, S. 2008, ApJ, 679, L97 [NASA ADS] [CrossRef] [Google Scholar]
  39. Leake, J. E., & Arber, T. D. 2006, A&A, 450, 805 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Leake, J. E., Arber, T. D., & Khodachenko, M. L. 2005, A&A, 442, 1091 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Leake, J. E., Lukin, V. S., Linton, M. G., & Meier, E. T. 2012, ApJ, 760, 109 [NASA ADS] [CrossRef] [Google Scholar]
  42. Leenaarts, J., Carlsson, M., Hansteen, V., & Rutten, R. J. 2007, A&A, 473, 625 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Leenaarts, J., Carlsson, M., Hansteen, V., & Gudiksen, B. V. 2011, A&A, 530, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180 [NASA ADS] [CrossRef] [Google Scholar]
  45. Mac Low, M.-M., Norman, M. L., Konigl, A., & Wardle, M. 1995, ApJ, 442, 726 [NASA ADS] [CrossRef] [Google Scholar]
  46. Martínez-Sykora, J., De Pontieu, B., & Hansteen, V. 2012, ApJ, 753, 161 [NASA ADS] [CrossRef] [Google Scholar]
  47. Martínez-Sykora, J., De Pontieu, B., Hansteen, V., & Carlsson, M. 2015, R. Soc. London Philos. Trans. Ser. A, 373, 40268 [Google Scholar]
  48. Martínez-Sykora, J., De Pontieu, B., Carlsson, M., & Hansteen, V. 2016, ApJ, 831, L1 [NASA ADS] [CrossRef] [Google Scholar]
  49. Martínez-Sykora, J., De Pontieu, B., Hansteen, V. H., et al. 2017a, Science, 356, 1269 [NASA ADS] [CrossRef] [Google Scholar]
  50. Martínez-Sykora, J., De Pontieu, B., Carlsson, M., et al. 2017b, ApJ, 847, 36 [NASA ADS] [CrossRef] [Google Scholar]
  51. Martínez-Sykora, J., Leenaarts, J., De Pontieu, B., et al. 2020a, ApJ, 889, 95 [CrossRef] [Google Scholar]
  52. Martínez-Sykora, J., De Pontieu, B., de la Cruz Rodriguez, J., & Chintzoglou, G. 2020b, ApJ, 891, L8 [CrossRef] [Google Scholar]
  53. Masson, J., Teyssier, R., Mulet-Marquis, C., Hennebelle, P., & Chabrier, G. 2012, ApJS, 201, 24 [NASA ADS] [CrossRef] [Google Scholar]
  54. Mellon, R. R., & Li, Z.-Y. 2009, ApJ, 698, 922 [NASA ADS] [CrossRef] [Google Scholar]
  55. Mestel, L., & Spitzer, L. J. 1956, MNRAS, 116, 503 [NASA ADS] [CrossRef] [Google Scholar]
  56. Meyer, C. D., Balsara, D. S., & Aslam, T. D. 2012, MNRAS, 422, 2102 [NASA ADS] [CrossRef] [Google Scholar]
  57. Mitchner, M., & Kruger, C. H. 1973, Partially ionized gases (John Wiley & Sons Inc.) [Google Scholar]
  58. Nakamura, F., & Li, Z.-Y. 2008, ApJ, 687, 354 [NASA ADS] [CrossRef] [Google Scholar]
  59. Ni, L., Kliem, B., Lin, J., & Wu, N. 2015, ApJ, 799, 79 [NASA ADS] [CrossRef] [Google Scholar]
  60. Nóbrega-Siverio, D., Moreno-Insertis, F., & Martínez-Sykora, J. 2018, ApJ, 858, 8 [NASA ADS] [CrossRef] [Google Scholar]
  61. Nóbrega-Siverio, D., Moreno-Insertis, F., Martínez-Sykora, J., Carlsson, M., & Szydlarski, M. 2020, A&A, 633, A66 [CrossRef] [EDP Sciences] [Google Scholar]
  62. Olluri, K., Gudiksen, B. V., & Hansteen, V. H. 2013, AJ, 145, 72 [NASA ADS] [CrossRef] [Google Scholar]
  63. Olluri, K., Gudiksen, B. V., Hansteen, V. H., & De Pontieu, B. 2015, ApJ, 802, 5 [NASA ADS] [CrossRef] [Google Scholar]
  64. O’Sullivan, S., & Downes, T. P. 2007, MNRAS, 376, 1648 [NASA ADS] [CrossRef] [Google Scholar]
  65. Padoan, P., Zweibel, E., & Nordlund, Å. 2000, ApJ, 540, 332 [NASA ADS] [CrossRef] [Google Scholar]
  66. Pattle, R. E. 1959, Q. J. Mech. Appl. Math., 12, 407 [CrossRef] [MathSciNet] [Google Scholar]
  67. Rutten, R. J. 2019, Sol. Phys., 294, 165 [NASA ADS] [CrossRef] [Google Scholar]
  68. Salmeron, R., & Wardle, M. 2008, MNRAS, 388, 1223 [NASA ADS] [Google Scholar]
  69. Schultz, D. R., Krstic, P. S., Lee, T. G., & Raymond, J. C. 2008, ApJ, 678, 950 [NASA ADS] [CrossRef] [Google Scholar]
  70. Shelyag, S., Khomenko, E., de Vicente, A., & Przybylski, D. 2016, ApJ, 819, L11 [NASA ADS] [CrossRef] [Google Scholar]
  71. Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23 [Google Scholar]
  72. Soler, R., Oliver, R., & Ballester, J. L. 2009, ApJ, 699, 1553 [NASA ADS] [CrossRef] [Google Scholar]
  73. Spitzer, L. 1978, Physical Processes in the Interstellar Medium (John Wiley & Sons) [Google Scholar]
  74. Stein, R. F., & Nordlund, Å. 2000, Sol. Phys., 192, 91 [NASA ADS] [CrossRef] [Google Scholar]
  75. Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117 [NASA ADS] [CrossRef] [Google Scholar]
  76. Vranjes, J., & Krstic, P. S. 2013, A&A, 554, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Vranjes, J., Poedts, S., Pandey, B. P., & de Pontieu, B. 2008, A&A, 478, 553 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Wardle, M. 1999, MNRAS, 307, 849 [NASA ADS] [CrossRef] [Google Scholar]
  79. Yoon, J.-S., Song, M.-Y., Han, J.-M., et al. 2008, J. Phys. Chem. Ref. Data, 37, 913 [CrossRef] [Google Scholar]
  80. Zaqarashvili, T. V., Khodachenko, M. L., & Rucker, H. O. 2011, A&A, 534, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Zel’dovich, Y. B., & Raizer, Y. P. 1967, Physics of Shock Waves and High-temperature Hydrodynamic Phenomena (New York: Academic Press) [Google Scholar]
  82. Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050 [NASA ADS] [CrossRef] [Google Scholar]
  83. Zweibel, E. G. 2002, ApJ, 567, 962 [NASA ADS] [CrossRef] [Google Scholar]
  84. Zweibel, E. G. 2015, in Magnetic Fields in Diffuse Media, eds. A. Lazarian, E. M. de Gouveia Dal Pino, & C. Melioli, Astrophys. Space Sci. Libr., 407, 285 [NASA ADS] [CrossRef] [Google Scholar]
  85. Zweibel, E. G., & Josafatsson, K. 1983, ApJ, 270, 511 [CrossRef] [Google Scholar]
  86. Zweibel, E. G., Lawrence, E., Yoo, J., et al. 2011, Phys. Plasmas, 18, 111211 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Abundances (A) (Gustafsson et al. 1975) and atomic mass (Z) for the 16 elements (El.) used in the EOS of Bifrost.

In the text

All Figures

thumbnail Fig. 1.

Elastic scattering cross section, σni, implemented in the Bifrost code to calculate ηamb. The different lines correspond to H–p collisions including temperature dependency and charge exchange (solid red line), H–e (pink), He–p (yellow), He–He+ (blue), He–e (green), H2–p (black), and H2–e (gray). For comparison purposes, the constant σni curve has been added, corresponding to H–p collisions frequently used in the literature (dashed red line).
In the text
thumbnail Fig. 2.

Tabulated η amb ∗ $ {\eta^{\ast}_{{\rm amb}}} $ as a function of density, ρ, and internal energy per unit volume, ϵ. The values for ϵ are shifted by 5 eV to obtain a positive internal energy when having hydrogen molecules (see main text). Temperature isocontours are superimposed as dashed lines for (from left to right) 3 × 103 K, 6 × 103 K, 104 K, 2 × 104 K, and 105 K.
In the text
thumbnail Fig. 3.

Tabulated η amb ∗ $ {\eta^{\ast}_{{\rm amb}}} $ as a function of the density, ρ, and temperature, T, for different EOS. Panel A: Complete EOS of Bifrost. Panel B: EOS only considering atoms of H, He, Na, Si, Mg, K, and molecules of H2. Panel C: EOS taking only H, H2, and He into account. Panel D: EOS including only atomic H.
In the text
thumbnail Fig. 4.

Plot showing f(n, ν) vs. n for different values of ν (dashed curves). Solid lines represent f(n, ν) = n2 (red) and f(n, ν) = n (black).
In the text
thumbnail Fig. 5.

Plot showing f(n, ν) versus n for different values of ν (dashed curves). Solid lines represent the boundaries of the parameter-space considering the constraints described in the text related to optimization and stability, namely f(n, ν) = n2, f(n = 10, ν) and, f ( n , ν = 5 − 2 ) $ f(n,\nu=\sqrt{5}-2) $.
In the text
thumbnail Fig. 6.

Two-dimensional map of C/n showing the combinations of n and ν that verify Eq. (23) for C >  4. In the image the horizontal dashed lines indicate the values ν = 5 − 2 $ \nu = \sqrt{5} - 2 $ (red) and the value of ν for the maximum n = 10 (ν ≈ 0.00185, black).
In the text
thumbnail Fig. 7.

Results of the STS validation test for the self-similar solution. Top panel: evolution of the maximum of the magnetic field with time (black asterisks) and power-law fit (red line). Bottom panel: evolution of Xb with time (black asterisks) and power-law fit (red line). Above each panel is shown the resulting fit formula.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.