Effect of ion-neutral collisions in simulations of emerging active regions
aa r X i v : . [ a s t r o - ph . S R ] D ec Effect of ion-neutral collisions in simulations of emergingactive regions
James E. Leake
College of Science, George Mason University, 4400 University Drive, Fairfax,Virginia [email protected] andMark G. Linton
U.S. Naval Research Lab 4555 Overlook Ave., SW Washington, DC 20375.
ABSTRACT
We present results of 2.5D numerical simulations of the emergence ofsub-surface magnetic flux into the solar atmosphere, with emerging fluxregions ranging from 10 to 10 Mx, representing both ephemeral andactive regions. We include the presence of neutral Hydrogen in the govern-ing equations, improve upon previous models by including the ionizationin the equation of state, and use a more realistic convection zone model.We find that ionization and recombination of plasma during the rise of aconvection zone flux tube reduces the rise speed of the tube’s axis. Thepresence of neutral Hydrogen allows the effective flow of mass across field-lines, by the addition of a Pedersen resistivity to the generalized Ohm’slaw, which dissipates current perpendicular to the magnetic field. Thiscauses an increase of up to 10% in the amount of magnetic in-plane fluxsupplied to the corona and a reduction of up to 89% in the amount ofsub-surface plasma brought up into the corona. However, it also reducesthe amount of free magnetic energy supplied to the corona, and thus doesnot positively affect the likelihood of creating unstable coronal structures.
Subject headings:
CMEs, Flux Emergence, MHD
1. INTRODUCTION
Coronal mass ejections (CMEs) and flares are the most energetic manifestationsof solar activity. It is generally accepted that the energy required for these events 2 –is stored as magnetic energy in the corona. This free magnetic energy is most likelystored in the strongly sheared flux of filament channels (see, e.g., reviews by Forbes(2000); Klimchuk (2001); Linton & Moldwin (2009)). These strongly sheared struc-tures have a significant component of the magnetic field parallel to the neutral line,compared to a potential configuration where the field is perpendicular to the neutralline. In the strong complex regions that are sources of fast CMEs, these filamentchannels form with the active region (Feynman & Martin 1995) and hence the shearfield must emerge with the flux. Given that shear formation is a fundamental driver ofthe eruption in CME models such as the magnetic breakout model (Antiochos et al.1999) and the flux cancellation model (Amari et al. 2000), it is clear that flux emer-gence must be explicitly included in such CME models for them to be physicallyrigorous.Most CME models, such as those of Antiochos et al. (1999), Amari et al. (2000),Chen & Shibata (2000), and Fan & Gibson (2007) are constrained by the need toextend the simulation domain to at least a few solar radii. They therefore do notmodel the lower solar atmosphere. The lower boundary of these simulations has atypical density of 3 × − g/cm and a typical plasma- β < .
2, where β = µ PB , (1) P is the gas pressure, B is the magnetic field strength, and µ is the permeabilityof free space. The photosphere, in contrast, has a density of around 3 × − g/cm (Vernazza et al. 1981) and a β > β convection zone to the low corona whereit is required to drive many CME models.The problem of flux emergence has been studied extensively for a few decades.The difficulty in studying flux emergence is that a simulation must couple the high β ,high density convection zone, where plasma pressure forces dominate, and the corona,which is low β and magnetically dominated. The basic idea of flux emergence is thatsub-surface flux tubes are created by dynamo actions at the base of the convectionzone. Indeed recent studies suggest that deep convective layers control the pattern oflarge-scale solar activity (Arkhypov et al. 2011, 2012). These flux tubes are assumedto have acquired sufficient twist to survive the convective motions as they rise tothe surface. The subsequent expansion into the atmosphere and its interaction withboth a field-free corona and pre-existing coronal fields has been studied extensively inrecent years, and a comprehensive list can be found in the review of Archontis (2008). 3 –In Leake et al. (2010) we showed that emerging two and a half dimensional (2.5D)sub-surface flux tubes did not supply enough free energy to the corona (less than 20%of the magnetic energy was in the shear component of the field) and hence the re-sulting coronal magnetic field did not erupt. We concluded that in order to create aneruption a method is needed which allows the transfer of significant shear field intothe corona by removing the dense material which is inhibiting the rise of the emergingstructures. This can be achieved in two ways. One method is the draining of plasmaalong the axis of the emerging flux tube. Recent three dimensional (3D) simulations offlux emergence have shown that the drainage along strongly arched flux tube axes canaid the rise of magnetic field into the corona (Hood et al. 2009; MacTaggart & Hood2009b). Also, complex 3D motions associated with shear flows during flux emer-gence have helped increase the amount of free energy that emerges into the corona(Manchester et al. 2004; Archontis & T¨or¨ok 2008; MacTaggart & Hood 2009a). Asecond method is the drainage of neutral plasma across field lines. Using a modelfor the support of a prominence in a partially ionized solar atmosphere, Gilbert et al.(2002) predicted the amount of vertical draining of neutral atoms from prominencestructures, and also found evidence of cross-field diffusion of neutral material in solarfilaments to support this model (Gilbert et al. 2007). In this paper we will explore thedrainage of neutral material across fieldlines during flux emergence, using a modifiedMHD model which includes the partial ionization of the solar atmosphere.In the simulations of Leake et al. (2010) and the majority of previous simula-tions, the numerical models use the magnetohydrodynamics (MHD) equations fora fully ionized plasma. For most of the convection zone and corona this is a validapproach. However, in the upper 3 Mm of the convection zone, as well as the pho-tosphere and the chromosphere, the temperature is low enough that the plasma isnot fully ionized. It has been shown that the effects of ion-neutral collisions cre-ate an anisotropic resistivity in the single-fluid equations (Cowling 1957; Braginskii1965). This Pedersen resistivity acts only on currents perpendicular to the mag-netic field, and has been shown to be up to 12 orders of magnitude larger thanthe parallel Spitzer, or Coulomb, resistivity in the chromosphere (Khodachenko et al.2004). The anisotropic dissipation of perpendicular currents by this Pedersen resistiv-ity has been included in the study of the damping of MHD waves in the chromosphere(de Pontieu 1999; Leake et al. 2005), and flux emergence in both 2.5D and 3D simu-lations (Leake & Arber 2006; Arber et al. 2007). In those flux emergence simulations,it was shown that the Pedersen resistivity led to increased rates of flux emergencedue to dissipation by ion-neutral collisions, and increased collisional heating.In the simulations presented here we improve on two of the limitations of the par-tially ionized flux emergence simulations of Leake & Arber (2006) and Arber et al. 4 –(2007). The first is related to the equation of state used. In the simulations ofLeake & Arber (2006) and Arber et al. (2007), the ionization fraction was calculatedas a function of density and temperature, but the specific internal energy density( ǫ ) did not include a contributing term from the ionization fraction: The equationof state in those simulations did not include the change in internal energy densitydue to ionization and recombination. In the simulations presented here, this term isspecifically included in the equation of state. As will be shown in the following sec-tion, this approach therefore also includes a correct description of a partially ionizedconvection zone which is adiabatically stratified.The second improvement is to use flux tubes which initially contain axial mag-netic flux comparable to observations of real active regions on the Sun. The simula-tions of Leake & Arber (2006) and Arber et al. (2007) use flux tubes in the convectionzone with axial fluxes of less than 10 Mx. The axial flux in these tubes is two ordersof magnitude less than that typically measured in active regions on the surface of theSun (10 Mx). More recently, Cheung et al. (2010), Rempel (2011), and Fang et al.(2012) have simulated the emergence of active region size flux tubes ( > Mx).However, those simulations do not include the effects of Pedersen resistivity on theflux emergence process. We model the emergence of large-scale tubes with axialfluxes of 10 Mx, as well as the more typical small scale flux tubes seen in the sim-ulations of Leake & Arber (2006), Arber et al. (2007), and Leake et al. (2010), usingour partially ionized plasma model.
2. NUMERICAL METHOD2.1. Equations
We solve the standard resistive MHD equations, modified to include the effectsof partial ionization. These are solved numerically using the Lagrangian remap code
Lare2d (Arber et al. 2001). The equations are obtained by summing the equationsfor ions ( i ) electrons ( e ) and neutrals ( n ). Hence the total mass density ( ρ ), gaspressure ( P ) and internal specific energy density ( ǫ ) are obtained by summing overthe three species, e.g., ρ = P a ρ a = P a m a n a where m a and n a are the mass andnumber density, respectively, of each species ( a = i, e, n ). The average velocity, v , isdefined by ρ v = P a ρ a v a . The neutral fraction is defined by ξ n = n n / ( n n + n i ) . (2) 5 –The equations are given below in Lagrangian form, using SI units: DρDt = − ρ ∇ . v , (3) D v Dt = − ρ ∇ P + 1 ρ j ∧ B + g + 1 ρ ∇ . S , (4) D B Dt = ( B . ∇ ) v − B ( ∇ . v ) − ∇ ∧ ( η j k ) − ∇ ∧ ( η p j ⊥ ) , and (5) DǫDt = − Pρ ∇ . v + ς ij S ij + ηj k + η p j ⊥ , (6)where j k and j ⊥ are the parallel and perpendicular current vectors, respectively, andare defined as j k = ( B . j ) B | B | , (7) j ⊥ = B ∧ ( j ∧ B ) | B | . (8)Here the current density is defined by j = ∇ ∧ B /µ and µ is the permeability offree space. The total gas density, pressure, and internal specific energy density aredefined at the center of each numerical cell. The magnetic field, B , is defined at cellfaces, and the velocity is defined at cell vertices. This staggered grid preserves ∇ . B during the simulation. The gravitational acceleration is denoted by g and is set to thevalue of gravity at the mean solar surface (274 m / s ). S is the stress tensor which hascomponents S ij = ν ( ς ij − δ ij ∇ . v ), with ς ij = ( ∂v i ∂x j + ∂v j ∂x i ) . The viscosity ν is set to3 × kg m − s − , and δ ij is the Kronecker delta function. The Coulomb resistivity, η , is given by η = m e ( ν ′ ei + ν ′ en ) n e e , (9)where m e and e are the mass and charge of the electron, respectively. The effectivecollisional frequencies for collisions of electrons with ions and neutrals, respectively,are ν ′ ei and ν ′ en . The Pedersen resistivity, η p , is given by η p = η + ξ n | B | α n . (10)The quantity α n is calculated using α n = m e n e ν ′ en + m i n i ν ′ in where ν ′ in is the effec-tive collisional frequency for collisions of ions with neutrals. We refer the reader toLeake et al. (2005) for the method of calculating the effective collisional frequencies.For this partially ionized plasma, the total pressure, P , and the specific internalenergy density, ǫ , can be written as P = ρk B T /µ m , and (11) 6 – ǫ = k B Tµ m ( γ −
1) + (1 − ξ n ) X i m i , (12)respectively, where k B is Boltzmann’s constant, γ is 5/3, and X i = 13 . ξ n , is afunction of temperature, T , itself and so Equation (12) is solved implicitly for T ateach time step. The reduced mass, µ m , is µ m = m i / (2 − ξ n ) . (13)Here m i = m f m p , where m p is the mass of a proton, and m f = 1 .
25 is a pre-factorwhich is designed to include the effects of heavier elements and, as will be shown in thenext section, will help reconcile our initial conditions with more realistic theoreticalmodels of the solar convection zone.To calculate ξ n we use a simple model based on the modified Saha equation(Saha 1921), which takes into account the fact that the chromosphere is not in localthermodynamic equilibrium (Brown 1973). This equation can be solved for the steadystate solution of the ionization equation (Athay & Thomas 1961) to give n i n n = f ( T ) b ( T ) , (14)where (15) f ( T ) = (2 πm e k B T ) γ − h exp( − X i k B T )and (16) b ( T ) = TwT R exp (cid:20) X i k B T ( TT R − (cid:21) . Here T R is the temperature of the photospheric radiation field, w is its dilution factor,and h is the Planck constant. Below the surface, T R = T and w = 1 so that b ( T ) = 1.Above the surface, T R = 6420 K and w = 0 .
5. Given n i /n n , the ratio of the numberdensity of neutrals to number density of ions can be calculated from r = n n n i = 12 − s(cid:18) ρ/m i n i /n n (cid:19)! , (17)and the neutral fraction, ξ n = n n / ( n n + n i ), is ξ n = r r . (18) 7 –For the convection zone and lower solar atmosphere, the fully ionized value ofthe Spitzer/Coulomb resistivity gives a magnetic Reynold’s number R m = U L π/c η (where U and L are typical measures of velocity and length) much larger than unity( > ). Thus typical numerical diffusion in MHD simulations of flux emergenceexceeds the classical diffusion, and hence η = 0 is often assumed. On the otherhand Leake & Arber (2006) showed that the value of the Pedersen resistivity givesan effective magnetic Reynold’s number below unity in the chromosphere for fluxemergence simulations, and that the increased diffusion due to the Pedersen resistivityis significantly larger than both the diffusion due to the Coulomb resistivity andtypical numerical diffusion.The equations are solved in 2.5D: the simulation box is 2D, with x and y beingindependent variables and z being ignorable, but all three components of the vectorvariables are evolved. We require a background stratification which represents the solar convection zoneand atmosphere. We first define a temperature profile which is chosen to resemblethe real Sun. Then with a specific value of density (or pressure) at the surface, thehydrostatic equilibrium equation can be solved to give the density and pressure inthe simulation domain.Above the surface ( y > T ( y ) = T ph + ( T cor − T ph )2 (tanh y − y cor w tr + 1) , (19)where T ph = 5778 K, T cor = 150 T ph , y cor = 3 .
75 Mm and w tr = 0 .
75 Mm. Thisprofile is designed to resemble semi-empirical atmospheric models of the quiet Sundeveloped by Vernazza et al. (1981). Below the surface, an adiabatically stratifiedgas in hydrostatic equilibrium is used as a model for the atmosphere: The vertical( y ) temperature gradient is equal to the adiabatic temperature gradient (Stix 2004).In this way the atmosphere is marginally stable to convection. For a fully ionizedplasma this yields a temperature profile in the convection zone of (cid:18) dTdy (cid:19) a = − µ m gk B γ − γ . (20)Previous simulations of flux emergence which assume the plasma to be fullyionized typically use the neutral limit of the reduced mass of µ m = m i , rather than 8 – µ m = m i /
2, which is strictly the fully ionized limit. As we will show, this choice ismade to ensure realistic coronal densities result from the solution to the hydrostaticequilibrium equation.Removing the assumption that the plasma is fully ionized yields a different resultfor the initial equilibrium. Rising and falling packets of plasma can now ionize andrecombine, and do not expand and contract ideally.We want to derive an equation for (cid:16) dTdy (cid:17) a , the adiabatic temperature gradient fora partially ionized plasma. This can be done by considering the adiabatic change ofa volume V of partially ionized plasma. Using Equations (12) and (13), and defining n = n i + n n , and ζ = n i /n = 1 − ξ n , we can express the energy in volume V as E = ρǫV = m i nV (cid:18) k B T (1 + ζ )( γ − m i + ζ X i m i (cid:19) = N (cid:18) k B T (1 + ζ )( γ −
1) + ζ X i (cid:19) , (21)where N = nV is the total number of ions and neutrals and does not change.Adiabatic evolution assumes no heat exchange: dQ = 0. Using the first law ofthermodynamics ( dQ = dE + P dV ), this implies dE = − P dV. (22)Differentiating Equation (21) gives dE = 1 γ − N k B dT (1 + ζ ) + 1 γ − N k B T dζ + N X i dζ . (23)Writing the total pressure as P = Nk B T (1+ ζ ) V we can differentiate V to obtain − dVV = − (cid:18) − dPP + dTT + dζ ζ (cid:19) . (24)Inserting Equations (23) and (24) into Equation (22), and using P = Nk B T (1+ ζ ) V gives dPP = γγ − dTT + θ dζ ζ , (25)where θ ≡ γγ − X i k B T . (26)This is the condition for the adiabatic change for a partially ionized plasma, witharbitrary ionization state ζ . For practical purposes we would like to eliminate dζ y < ζ : ζ − ζ = A T γγ − P e − X i /k B T , (27)where A contains physical constants only. This can be differentiated to give2 ζ (1 − ζ ) dζ = ζ − ζ (cid:18) − dPP + θ dTT (cid:19) , (28)or dζ ζ = (cid:18) − dPP + θ dTT (cid:19) ζ (1 − ζ )2 , (29)which can be inserted into Equation (25) to give dTT = dPP θ ζ (1 − ζ )2 γγ − + θ ζ (1 − ζ )2 ! . (30)Dividing by dy , and using dPdy = − ρg and P = ρk B T /µ m , gives (cid:18) dTdy (cid:19) a = − µ m gk B θ ζ (1 − ζ )2 γγ − + θ ζ (1 − ζ )2 ! . (31)If ζ = 0 or ζ = 1, then we recover the usual definition for a fully neutral or fullyionized plasma, as expressed in Equation (20).In summary, to construct a model atmosphere which is in hydrostatic equilibriumand whose temperature gradient below the photosphere ( y = 0) matches Equation(31), we must perform the following procedure:1. For y ≥ T ( y ) from Equation (19).2. For the entire domain, initially set ζ = 1. Set T ( y = 0) = 5778 K and ρ ( y =0) = 3 . × − kg / m , taken from the standard solar model presented in Stix(2004).3. For y <
0, numerically solve Equation (31) with boundary condition T ( y = 0),using the current ζ ( y ). This gives T ( y < dPdy = − ρg with boundarycondition ρ ( y = 0), using the current T ( y ) and ζ ( y ), and the equation P = ρk B T /µ m . This gives ρ ( y ) for the whole domain. 10 –5. Use the Saha equation (14) with the current T ( y ) and ρ ( y ) to get ζ ( y ), andhence the reduced mass µ m ( y ) = m i ζ ( y ) for the whole domain.6. Repeat steps 3-5 until the reduced mass µ m ( y ) converges for all y .In order to show how this approach differs from the standard approach used inprevious fully ionized simulations, we present the results of this process for threemodels. A summary of these simulations, along with simulations described in latersections, is shown in Table 1. Simulation 0 uses the fully ionized plasma model withthe ionized limit of µ m = m i /
2. Simulation 1 again uses the fully ionized plasmamodel but with the neutral limit of µ m = m i . This is the case used by previous fullyionized simulations (Archontis 2008). We describe Simulation 2 later, and do notconsider it at this point. Simulation 3 uses the partially ionized plasma model withan ionization which is allowed to change in space and time, based on the modifiedSaha equation. The results are shown in Figure 1. Also shown are curves taken from astandard solar model which includes the transfer of heat by convection and radiationin the partially ionized plasma of the solar convection zone (Stix 2004), and a 1Dsemi-empirical atmospheric model by Vernazza et al. (1981). We can see from Figure1 that the temperature of the partially ionized model (Simulation 3) best matchesthe standard solar model in the convection zone. The temperature of the two fullyionized models (solid and dotted lines) depart from the standard solar model curvein the convection zone, as the ionization cannot change in these two models. In allof the models we use m f = 1 .
25 in order to include the effect of heavier elements,and this ensures that at a height of y = 0 the models have a reduced mass of 1 . m p which is the same as the standard solar model. However, only the partially ionizedmodel is able to capture the variation in reduced mass with height in the convectionzone and corona.While the temperature in the convection zone in the partially ionized model(Simulation 3) matches the standard solar model quite closely, there is some differencebetween the density and pressure profiles in the convection zone for the partiallyionized model and the standard solar model. We expect that this difference is becausethe models here do not include the ionization of Helium, which is taken into accountin the standard solar model. To test whether this difference was caused by the factthat the convection zone is not adiabatically stable throughout its depth, we alsotested a model which had dT /dy = F ( dT /dy ) a where F = F ( y ) was a functionwhich reproduced the non-adiabatic temperature profile of the standard solar modelof Stix (2004). This unstable convection zone profile did not produce any significantimprovements in the density and pressure profiles. 11 –Looking at the density and pressure at and above 2 Mm in Figure 1, there is alarge difference between the fully ionized model with the ionized limit of µ m = m i / − kg / m , whereas Simulation 0 has a density 4 ordersof magnitude higher than this. Simulation 0 was included here to show why previousauthors who have performed fully ionized simulations use a reduced mass of µ m = m i rather than µ m = m i /
2. Using the fully ionized limit of µ m = m i /
3. EMERGENCE OF SMALL SCALE FLUX TUBES3.1. Initial conditions
In order to compare our partially ionized simulations to previous small-scale fluxemergence studies (Magara 2001; Fan 2001; Manchester et al. 2004), we choose aninitial flux tube with a very small axial flux compared to the flux of a real activeregion. In the next section we will address the emergence of flux tubes with activeregion size flux (10 Mx). We insert a cylindrical magnetic flux tube into the modelconvection zone at x = 0 and y = y t = − . B z , and twist field, B θ , for the tube are given by B z = B e − r /a , and (32) B θ = qrB z , (33)respectively, where r = p x + ( y − y t ) is the radial distance from the center of thetube, B is the axial field strength at the center ( r = 0), a = 0 . q = − /a . For a field strength of B = 6000 G, this tube hasan axial flux of 5 × Mx. This tube contains 25% less axial flux than our previous2.5D simulations (Leake et al. 2010). Our simulation domain is -3 to 42 Mm in y and-22.5 to 22.5 Mm in x . The numerical grid contains 640 ×
640 cells, at a uniformresolution of 70 km.To investigate the effects of partial ionization, we perform three separate simula-tions. Simulation 1 is the fully ionized plasma model (with the neutral reduced masslimit of µ m = m i ), as described in the previous section. To relate this work to thework of Leake & Arber (2006) we also include the model used in that paper, and callit Simulation 2. In Simulation 2, the ionization is calculated as a function of ρ and 12 – T using the Saha equation (14), but the second term in Equation (12) is ignored, sothat ionization effects are not included in the energy calculation. In addition, thismeans that the convection zone profile used in Leake & Arber (2006) does not includethe variation of ionization in the equation for the adiabatic temperature gradient: ituses Equation (20) rather than Equation (31). The initial background stratificationfor Simulation 2 is identical to Simulation 1, and so is not included in Figure 1. Sim-ulation 3, as described in the previous section, uses the partially ionized model withthe full equation of state, and includes ionization and recombination in the equationfor the adiabatic temperature gradient. As shown in Figure 1, Simulation 3 has ahigher gas pressure at -1.8 Mm, where the flux tube is located, compared to Simula-tion 1. So for the same field strength, the flux tube will have a higher plasma β inthe partially ionized simulation (3). We therefore vary the magnetic field so that theplasma β has the same value at the center of the tube in all three simulations. Therise speed of a buoyant flux tube scales as 1 / √ β (Longcope et al. 1996), and so bymatching β across simulations we ensure that the rise speeds are initially the same.For Simulations 1 and 2 we use B = 6000 G. For Simulation 3 we use B = 10900 G.The plasma β , calculated as β = µ p ( y tube ) /B ( x, y tube ), for the initial conditions ofthese three models is shown in Figure 2, Panel (a).In all these simulations we add a perturbation P ( r ) to the background plasmapressure P ( y ) such that ( ∇ P ) r = ( j ∧ B ) r , so that the tube is in radial force balance.Assuming that the flux tube is in thermal equilibrium with its surroundings makesthe tube less dense than the surrounding plasma and initiates its buoyant rise to thesurface. With this “isothermal” assumption, the density perturbation, ρ ( r ), whichis added to the background density ρ ( y ), is given by ρ /ρ = P /P ∼ − / β. (34)This perturbation for the three simulations is shown in Figure 2, Panel (b). In the three simulations the flux tubes rise to the surface due to their initialbuoyancy and undergo a degree of horizontal expansion as they meet the convectivelystable photosphere. A contact layer is created when the tube’s field meets the pho-tosphere. As more field builds up at the surface, the layer becomes unstable to themagnetic Rayleigh Taylor instability (Acheson 1979). The criterion for this instabilitycan be written as − H p ∂∂y ln ( B ) > H p k (cid:18) l n (cid:19) − γ βδ (35) 13 –(Archontis et al. 2004) where H p is the local pressure scale height and δ is the super-adiabatic excess (Stix 2004), which is the double logarithmic temperature gradient( ∂lnT /∂lnP ) minus its adiabatic value. In these simulations, δ is zero in the modelconvection zone, but negative in the model photosphere. The wavenumber k is thewavenumber parallel to the magnetic field in the horizontal plane, l is the wavenumberperpendicular to the field in the horizontal plane, and n is the wavenumber in thevertical direction. For the contact layer that is created, the left hand side in Equation(35) is positive, and acts to destabilize the layer. The last term on the right handside is also positive and acts to stabilize the layer.Figure 3 shows the rise of the flux tubes in the convection zone in the threesimulations. The contour lines show fieldlines in the plane, given by constant intervalsof A z , where A z is defined as the vector potential for the 2D vector ( B x , B y ) = ∇∧ A z ˆ e z with boundary condition A z ( ∞ ) = 0. The outer fieldlines shown here are the contoursof A z = 0 . A z ) where A z is always negative. The background color contourshows log( ρ × m / kg). The tubes have risen to the surface with comparable speeds,which is to be expected as they started with the same plasma β . The outer fieldlinesare a different shape in Simulation 3, the partially ionized model simulation, than inSimulations 1 and 2. This is due to an increased horizontal expansion at the top ofthe tube in Simulation 3 relative to Simulations 1 and 2. This vertical gradient inhorizontal expansion is caused by changes in ionization, which are not taken accountin Simulations 1 and 2.Figure 4 shows the height of the center and edge of the tube in all three simu-lations. The center of the flux tube is defined as the location of the local extremaof A z , as A z is initially a negative monotonically increasing function of radius forour twisted flux tubes. The edge of the flux tube is defined as the intersection of A z = 0 . A z ) with the x = 0 line. In Simulation 3 (dot-dashed line), the fluxtube’s center stops 0.5 Mm below the surface. On the other hand, in the fully ionizedsimulation (1) and the partially ionized simulation (2), the tube’s center rises to y = 0(the surface). Figure 4, Panel (b) shows the rise speed in the convection zone in allthree simulations. All three simulations show an identical rise speed initially, as theyhave the same plasma β . The rise speed in Simulation 3 slows down relative to thatof Simulations 1 and 2 due to changes in ionization during the rise.At around 1000-1020 s, all three simulations show the onset of the magneticbuoyancy instability. Figure 5 shows the stabilizing ( − γ βδ ) and the destabilizing( − H p ∂∂y ln ( B )) terms in Equation (35) for the three simulation at two different times.As magnetic field reaches the photosphere, the local plasma- β decreases, which re-duces the stabilizing term − γ βδ , until it becomes less than the destabilizing term 14 – − H p ∂∂y ln ( B ). At this point, the instability starts to develop, and magnetic fieldemerges into the atmosphere. Simulation 2 shows a much smaller gradient in mag-netic field than the other two simulations, indicated by a lower magnitude of thedestabilizing term in Figure 5, Panel (e). In Leake & Arber (2006) we showed thatthe magnetic field was able to slip through the plasma due to the Pedersen resistivity,which dissipated currents perpendicular to the field. In Simulation 2 we see the effectof this as a reduction in the vertical gradient of magnetic field in a layer above a heightof 0.6 Mm above the surface. We do not see this in Simulation 1 as the model doesnot include Pedersen resistivity in that simulation. At this point in time, we do notsee this in Simulation 3, as the rise of the tube in the convection zone in Simulation3 was slowed due to the ionization/recombination effects, but this reduction in thegradient in magnetic field is seen later in the simulation.Figure 6 shows the expansion of the magnetic field into the corona due to themagnetized Rayleigh Taylor instability, and Figure 7 shows the magnetic Reynoldsnumber R m ∼ | v ∧ B || η p j ⊥ | (36)for the two partially ionized simulations (at time 1359 s for Simulation 2, and attime 1724 s for Simulation 3). The fully ionized simulation has a theoretically infiniteReynold’s number (due to η p = 0). At around 0.3 to 0.8 Mm above the surface,the two partially ionized simulations show a value of R m <
1. This means that thefield’s motion is dominated by diffusion rather than advection, and the dissipationincreases the emergence of the field into the atmosphere. This result is in agreementwith the previous simulations in Leake & Arber (2006), even though in Simulation 3we have a different equation for the specific energy density ǫ which allows for changesin ionization.As shown in Figure 4, the upper edge of the tube expands further and fasterinto the corona in Simulations 2 and 3, as the increased dissipation due to ion-neutralcollisions allows the field to diffuse through the plasma without having to lift materialup. In the fully ionized simulation, the field must lift more plasma up with it, whichslows its rise. Figure 8 shows the amount of in-plane flux that has emerged into thecorona for each simulation, normalized to the initial amount of in-plane flux in eachflux tube. This in-plane flux is calculated byΦ( t ) = [ max ( A z ( x, y )) − min ( A z ( x, y ))] y>y (37)where y = 1 . t ) is normalized to | min ( A z ) | at t=0, so that itrepresents the fraction of initial in-plane flux that emerges above the height 1.2 Mm.All three simulations emerge above 60% of the initial in-plane flux above this height. 15 –Due to the Pedersen dissipation, the two partially ionized simulations have peaks inthe in-plane flux above 1.2 Mm higher than Simulation 1, the fully ionized simulation.In the fully ionized model, the expansion of the plasma due to the emergingfield is associated with a significant amount of cooling, and the plasma temperaturedrops to around 100 K within the emerging flux region, which can be seen in Figure9. The transition region is pushed outward from 2 Mm to 14 Mm. The ion-neutralcollision effects in the two partially ionized models (Simulations 2 and 3) are able tocounteract this expansive cooling, and so this amount of cooling is not seen in thetwo partially ionized simulations, although some cooling to 1000 K is still seen andthe transition region is still pushed outwards. There are two effects which contributeto this result. The first is collisional ion-neutral heating. The second effect is theslippage of mass through the field, due to the Pedersen resistivity, which allows thefield to expand without expanding and cooling the plasma. These two effects ensurephotospheric/chromospheric temperatures stay above 1000 K in the partially ionizedsimulations, as originally shown in Leake & Arber (2006). This result gives furtherevidence that ion-neutral collisions are an important source of heating in emergingactive regions. Radiative-MHD numerical investigations by Leenaarts et al. (2011)show that, in 2D simulations, unrealistic cooling occurs which can only be counter-acted by Joule heating. We have shown here that the additional Joule heating dueto ion-neutral collisions is one way to achieve this necessary effect.We now investigate the effect of the increased dissipation due to ion-neutralcollisions on the amount of mass and free energy supplied to the corona during fluxemergence. In Simulation 3, the density is higher at the initial tube location, asshown in Figure 1, and hence the flux tube contains more plasma than in the othertwo simulations (1 and 2). We therefore calculate the change in mass in the corona,above a certain height during the flux emergence, normalized to the amount of massinitially in the flux tube. This diagnostic effectively estimates the fraction of the fluxtube mass which is supplied to the corona. Figure 10 shows this diagnostic above1.2 Mm in the three simulations. This shows that at t=2900 s, the fully ionizedsimulation (1) initially has lifted up approximately 7.7 times more of the tube’s massabove 1.2 Mm than Simulation 3, but only 28% more than Simulation 2.Figure 11 shows the quantity | j ∧ B || j || B | , which is equivalent to sin( θ ), where θ is theangle between the current density and magnetic field, for the three simulations, duringthe post-emergence stage. A value of 0 means that all the currents are aligned withthe field, i.e., the field is force-free. This shows that the angle between the field andthe current density is lower inside the expanding tube in Simulation 3 compared withSimulation 1. The Pedersen dissipation acts on the cross-field currents, thus reducing 16 – | j ∧ B || j || B | . This reduction in the Lorentz forcing component of the current density in partexplains why so much less mass is supplied to the corona in Simulation 3, relativeto Simulation 1. Simulation 2 also shows lower values of | j ∧ B || j || B | inside the emergingactive region, relative to Simulation 1, again, explaining in part the reduction of masssupplied to the corona.It is worth noting here that the Pedersen resistivity in the partially ionized modelonly effectively dissipates perpendicular currents inside the emerging flux structure.The edge of the flux tube is a current sheet, with a large perpendicular current. How-ever, for this current sheet, where the magnetic field decreases to zero, the Pedersenresistivity falls to zero as the current increases, and so there is limited dissipation ofthe current sheet. In fact Arber et al. (2009) showed that for current sheets with noguide field, the Pedersen resistivity acts to sharpen currents, rather than dissipatethem.There are two effects of the ion-neutral collisions working in parallel in the twopartially ionized models during the emergence. Firstly, due to the Pedersen resistivity,the mass can ‘slip’ through the emerging magnetic field, and so less mass is lifted perunit amount of concave up flux in the partially ionized model relative to the fullyionized model. Secondly, cross-field currents are reduced by the Pedersen resistivity,and so the magnetic field cannot lift up as much mass in the partially ionized modelsas in the fully ionized model. For the parameters chosen for these simulation, we findthat the result of these two effects is the reduction of the amount of mass suppliedto the corona in both partially ionized models (2 and 3), relative to the fully ionizedsimulation (1).Figure 10 also shows the amount of energy in the shear component of the field( B z ) normalized to the total amount of energy in the field. E shear = R x =22 . Mmx = − . Mm R y =42 Mmy B z dxdy R x =22 . Mmx = − . Mm R y =42 Mmy | B | dxdy (38)where y = 1 . B z ) is coupled to the mass, and we see more shear field supplied to thecorona than in the other two simulations. At t=2900 s, Simulation 1 has a valueof E shear
10% higher than Simulation 2, and 25% higher than Simulation 3. Wepostulate that Simulation 1 can emerge more shear flux, even though this meansemerging more mass, because it has more Lorentz force.Comparing Simulations 2 and 3 shows that including ionization effects in theequation of state significantly affects the emergence of flux. There is more horizontal 17 –expansion during the emergence process, both in the convection zone and the coronawhen the full equation of state is used. By 2900 s, Simulation 2 has emerged 6 timesthe normalized mass into the corona that Simulation 3 has. Also, Simulation 2 hasemerged 15% more shear (in terms of E shear ) into the corona than Simulation 3.We also note that the conclusions made in Leake & Arber (2006), which used themodel of Simulation 2, are repeated here with Simulation 3, which uses a more self-consistent model for the partially ionized plasma, in that it includes instantaneousionization/recombination in the equation of state. These results include the increasein the amount of in-plane flux injected into the corona (Simulations 2 and 3 emergemore in-plane flux than Simulation 1), and the additional Joule heating which is vitalfor maintaining a chromospheric temperature above 1000 K.In summary, we have performed three simulations of an emerging flux tube witha single value of β . Simulation 1 was a fully ionized model. Simulation 2 was apartially ionized model without the full equation of state, and Simulation 3 was apartially ionized model with the full equation of state, and a correct treatment of theadiabatic convection zone. We have shown that if a fully ionized model is used, theneutral limit for µ m = m i should be used to get the correct density in the corona.These results show that although the presence of neutral Hydrogen can reducethe amount of mass supplied to the corona, which in principle can increase the likeli-hood of eruption, the amount of shear supplied to the corona in the partially ionizedsimulations is actually less than in the fully ionized simulations. Hence the likelihoodof eruption is less in the partially ionized simulations. In all three simulations, theamount of magnetic energy in the shear component of the field is less than 20% ofthe total magnetic energy in the field, and thus the coronal field is non-eruptive, justas in the simulations of Leake et al. (2010). For these flux tubes, the ‘slippage’ offield through the nearly neutral lower atmosphere does not positively increase thelikelihood of eruption.
4. EMERGENCE OF LARGE SCALE FLUX TUBES
The results of the previous section showed that including the ionization in theequation of state of the partially ionized model gave quantitatively different resultsfrom the model without the ionization terms. We performed the comparison in orderto put the results in context with the previous simulations of Leake & Arber (2006).We know that the model used in Simulation 3, which includes the ionization in theequation of state, is the more self-consistent model of the two partially ionized models,and so we drop the model used in Simulation 2, which does not include the ionization 18 –term in the equation of state. We study now the fully ionized model and the partiallyionized model in terms of the emergence of active region size magnetic flux tubes.The axial flux in the tube in the previous section’s simulations was less than 10 Mx, which is much smaller than the 10 Mx size of a typical active region sunspot(Priest 2003; Zirin 1998). We must reconcile our initial conditions with realistic solarvalues to better test the effect of partial ionization on emerging active regions.To create an active-region scale flux tube, we increase the tube radius to 1.5 Mmand place the tube at a depth of -4.95 Mm. The domain is increased to cover a range-56.25 to 56.25 Mm in x and -15 to 97.5 Mm in y . The numerical resolution is kept thesame as before at 70 km. The results of the previous section showed that flux tubeswith initially comparable β in the flux tube start to rise with comparable speeds. Inthis section we run four simulations, which are grouped into two pairs of comparable- β simulations. The first pair consists of a fully ionized model (Simulation 4), and apartially ionized model (Simulation 5). Both simulations start with a flux tube whichhas β = 4. The second pair consists again of one fully ionized model (Simulation 6)and one partially ionized model (Simulation 7), both having a flux tube β of 40. Thesimulations are briefly described in Table 1. The axial flux in the tube in Simulation5 is 1 . × Mx, and so is comparable with observed active region sunspots. Thetwo different values of β of 4 and 40 are chosen to cover a range of active regionformation time. The simulations from the previous section show that flux tubes withan initial plasma β of 2 emerge in less than 30 minutes, and so increasing the plasma β in the tube to 40, we expect emergence on a timescale of approximately 2 hours,which is somewhat closer to the general formation time of active regions on the Sun(12-48 hours). Figure 12 shows the initial β profiles and density perturbations, ρ /ρ ,for these four simulations.Figure 13 shows the in-plane magnetic field as regular contours of A z , and log ofdensity as a color contour, at two different times in Simulations 4 and 5 (the β = 4simulations). The first time (left panels) is when the tube has reached the surface,and built up enough magnetic field to trigger the magnetic buoyancy instability. Thesecond time (right panels) is at a later stage, when the outer 10% of the in-planeflux has reached approximately 10 Mm above the surface. Figure 14 shows the samephases of the emergence process but for the two simulations 6 and 7, where β = 40.The high β in these simulations means that the initial rise-speed is slow comparedto the low β simulations 4 and 5, and it takes longer to build up enough field at thesurface to trigger the instability which drives field into the corona.Figure 15, Panel (a) shows the height of the tube center and tube edge for thelarge scale simulations. As in the small scale simulations, the tube centers are con- 19 –fined to just below the surface, while the edges expand in to the field-free corona.Figure 15, Panel (b) shows the rise speeds in the convection zone. Firstly, we notethat these simulations confirm that the convection zone rise speed scales with 1 / √ β (Longcope et al. 1996); Taking the two fully ionized simulations (4 and 6), and cal-culating the ratio of the maximum rise speeds in the convection zone gives a valueof 3.3, which is approximately equal to the inverse of the ratio of √ β in those twosimulations ( p / β = 4 models, solid and dot-dash line respectively), the rise speed peaks at a value of 2.3 km/s in the partiallyionized model (Simulation 5), and a value of 2.9 km/s in the fully ionized model (Sim-ulation 4). This difference in rise speeds also occurred in the small scale simulations,and is a consequence of the ionization and recombination during the rise of the tubein the partially ionized convection zone. This result is repeated in the β = 40 simu-lations (Simulation 6: double-dot-dash line, and Simulation 7: long-dash line), withthe fully ionized simulation obtaining a higher rise speed than the partially ionizedsimulation. This may explain why the flux tube axis settles lower in the convectionzone in Simulation 7 compared to Simulation 6. This effect may be important for theformation of coronal structures in 3D, as it is thought that the evolution of the fluxtube axis is important for the formation of sheared structures in 3D.The amount of normalized in-plane flux (Φ) supplied to the corona above 1.2 Mmis shown in Figure 16. Comparing Φ( t ) for the low and high β simulations, we seethat more of the original in-plane flux emerges in the low β simulations (4 and 5) thanin the high β simulations (6 and 7). Still there is almost 50% of the original in-planeflux supplied to the corona in the high β case, which is comparable to the amountof flux which emerges above 1.2 Mm in the small-scale simulations of the previoussection. Comparing the fully ionized simulations to the partially ionized simulations(4 vs. 5 and 6 vs. 7 ) it is clear that the flux emerges earlier in the partially ionizedsimulations (5 and 7), compared to the fully ionized simulations (4 and 6).The resulting temperature profiles are shown in Figure 17. As in the small scalesimulations, the partially ionized simulations do not see the drastic cooling below1000 K in the chromosphere that is seen in the fully ionized simulations. This isa result of both the collisional heating, and the dissipation by Pedersen resistivityreducing the expansion of the plasma as the field expands, a result which we found inthe smaller-scale simulations of the previous section. We have therefore shown thatthis result applies to the larger active region and slower emergence of these large scalesimulations (4 through 7). 20 –Figure 18 shows the normalized change in mass above 1.2 Mm in the four largescale simulations. Comparing the two β = 4 simulations (4 and 5), we find that thefully ionized simulation supplies approximately 10 times more of the flux tube massto the corona than the partially ionized simulation, a result which is repeated for thetwo β = 40 simulations (6 and 7). This, again is a direct consequence of the Pedersendissipation in the partially ionized models.It is worth comparing the amount of flux tube mass supplied to the corona inthese large scale simulations to that in the small scale simulations. The two largescale β = 4 simulations supply approximately the same amount of flux tube massabove 1.2 Mm (10 − for the fully ionized model and 10 − for the partially ionized)as the small scale simulations, which have β ≈
2. The β = 40 large scale simulationssupply more than an order of magnitude less flux tube mass to the corona than the β = 4 simulations (5 × − for the fully ionized and 4 × − for the partially ionizedmodel). So in general the amount of flux tube mass supplied above 1.2 Mm duringflux emergence is inversely proportional to the initial β in the flux tubes at t=0, andhence is dependent on the strength of the magnetic field in the tubes, as one wouldexpect.Also shown in Figure 18 is E shear above 1.2 Mm. For the β = 4 simulations,the fully ionized simulation supplies 15% more free energy at t= 2900 s than thepartially ionized simulation. There is no significant difference between the two β = 40simulations. In these large scale simulations, less than 25% of the magnetic energyin the corona is in the shear component, which broadly agrees with the results of thefully ionized simulations of Leake et al. (2010). Although the free energy does exceed20% of the magnetic energy in the β = 4 simulations, it quickly falls below this value.The β = 40 simulations stay below 20% for the duration of the simulations. Theresulting coronal configurations in these 2.5D simulations are non-eruptive, despitethe decrease in mass supplied to the corona due to the presence of neutral materialin the lower solar atmosphere.The successful emergence of flux tubes with active region values of magnetic flux(10 Mx) has been achieved, and been shown to be qualitatively similar to the emer-gence of smaller scale flux tubes, which have been studied extensively. Furthermore,the rise speed of these large scale flux tubes is affected by the partially ionized regionof the solar atmosphere, as is the amount of mass and shear supplied to the corona,even though the initial tube radius (1.5 Mm) is comparable to the size of the regionwhere the Pedersen resistivity is important (from the surface to about 1 Mm abovethe surface). 21 –
5. CONCLUSIONS
We have investigated the effect of neutral Hydrogen in the lower solar atmosphereon the emergence of sub-surface magnetic flux into the solar atmosphere, using 2.5DMHD models modified to include a variable ionization state, and the anisotropyin Ohm’s law caused by ion-neutral collisions. Previous simulations using a fullyionized MHD model (Leake et al. 2010) have shown that the amount of free energyis insufficient to drive a CME when flux tubes emerge in this 2.5D cartesian setup.The presence of too much mass and the lack of a strong enough shear field yieldedno further expansion of the flux rope structure, which was ultimately constrained byoverlying field. Typically the maximum amount of energy in the shear field was 10-15% of the magnetic energy. The main aim of this paper was to investigate the effectof the partially ionized atmosphere on the likelihood of creating eruptive magneticfield from emerging convection zone flux tube.The single-fluid MHD equations were modified to include the presence of neu-tral Hydrogen. This led to a modified induction equation, where perpendicular cur-rents were dissipated by ion-neutral collisions and parallel currents were dissipatedby electron-ion and electron-neutral collisions. It also led to a source term relatedto ionization state in the equation for the specific internal energy density. Pre-vious simulations (Leake et al. 2005; Leake & Arber 2006; Arber et al. 2007) withthese effects used a simple equation of state which did not take into account ioniza-tion/recombination, and thus used a model for an adiabatic convection zone whichdid not include partial ionization. In this paper we included the change in ionizationin our equation of state, and therefore a more realistic convection zone profile.In the partially ionized simulations, the dissipation due to the Pedersen resistivityand the associated collisional heating was concentrated in the lower atmosphere, belowthe transition region. In the fully ionized simulations, the rapid expansion of theemerging field was associated with a rapid expansion and cooling of the plasma in theemerging region, which resulted in temperatures below 1000 K, which we consider tobe too low for the real Sun. In the partially ionized simulations, we did not see suchcooling below 1000 K. There are two effects which contribute to this. The first is thefact that as the field emerges into the corona in the partially ionized simulations, thePedersen resistivity allows the field to ‘slip’ through the plasma, and this reduces theamount of expansion and cooling seen in the emerging region. The second is thatany cooling will be countered by the collision dissipation. This result gives furtherevidence that ion-neutral collisions are an important source of heating in emergingactive regions, as was suggest in the radiative-MHD simulations of Leenaarts et al.(2011). 22 –The result of the increased dissipation in the lower atmosphere in the partiallyionized simulations led to a number of effects. Firstly, for the small scale simulationsand the β = 40 large scale simulations, the amount of in-plane flux supplied to thecorona increased compared to the fully ionized models. The Pedersen dissipation inthe partially ionized simulations led to a decrease in the Lorentz forcing componentof the current density inside emerging active regions, and a reduction in the masslifted into the corona by the emerging flux. In the fully ionized simulations, there wasmore of the flux tube mass supplied to the corona. Coupled to this was an increasein shear field supplied to the corona in the fully ionized simulation. As a result, thelikelihood of eruption was not increased by the effects of neutral Hydrogen in thelower atmosphere, but instead reduced.We performed studies with both small scale (radius of 0.3 Mm and fluxes of5 − × Mx), and large scale (radius of 1.5 Mm and fluxes of 10 − Mx)flux tubes. The larger flux tubes resulted in active regions of size 20 Mm, based onthe separation of the two intersections of the 10% flux contour with the surface y = 0.For the smaller flux tubes the active region size was typically 10 Mm. Interestingly,increasing the initial flux tube radius by a factor of 5 only resulted in an increase inthe resultant active region size by a factor of 2.Results of our small scale flux emergence simulations showed that by includingthe ionization in the equation of state, the rise speed was effectively reduced byionization and recombination in the upper convection zone, leaving the axis of the tubelower down in the convection zone after the emergence. Also in the simulations wherethe ionization was included in the equation of state, the amount of mass and shearflux supplied to the corona was decreased relative to the partially ionized simulationwhich did not include ionization in the equation of state.For the larger scale simulations we used two different vales for β in the flux tube.The timescale of emergence depended on the β value. For β = 4 the flux tube edgestook 50 minutes to reach a height of 35 Mm, whereas for β = 40 the emergence took100 minutes to reach the same height. This delay was almost entirely accounted for bythe slower rise speed in the convection zone, which scaled with p /β . In all the lowand high β simulations, more than 55% of the initial in-plane flux emerges into thecorona. Even though the flux tube’s radius in the large scale simulation was 10 timeslarger than the photospheric scale height of 150 km, the emergence was qualitativelysimilar to the emergence of smaller flux tubes.For these 2.5D simulations, the presence of neutral Hydrogen in the lower atmo-sphere, and the associated ‘slippage’ of emerging magnetic field (or equivalently theassociated motion of plasma across fieldlines) did not increase the amount of shear 23 –flux supplied to the corona, and therefore we conclude that the effects of neutral Hy-drogen do not increase the likelihood of eruption. We also conclude that 3D plasmamotions along the axis of the flux tubes are more likely to destabilize coronal mag-netic field in emerging active regions. However, the stability of flux ropes formedduring flux emergence, and the subsequent likelihood of eruption, will be dependenton the currents associated with the emerging fields, and therefore should be stronglydependent on the nature of the current dissipation mechanism. We propose in futurework to investigate both 3D effects and the effects of partial ionization. Note thatthese conclusions address how the presence of neutral Hydrogen affects the supply offree energy, or sheared field, into the corona by the emergence of magnetic flux into afield-free corona. Hence the only possible source of free energy in these simulations isthe newly emerging flux. In the case where a pre-existing sheared structure is alreadyformed in the corona, the emergence of new flux can play the role of destabilizing themagnetic field and causing an eruption.Acknowledgements: This work has been supported by the NASA Living With a Star,Solar & Heliospheric Physics programs, the ONR 6.1 Program, and the NRL-Hinodeanalysis program. The simulations were performed under a grant of computer timefrom the DoD HPC program. REFERENCES
Acheson, D. J. 1979, Sol. Phys., 62, 23Amari, T., Luciani, J. F., Mikic, Z., & Linker, J. 2000, ApJ, 529, L49Antiochos, S. K., DeVore, C. R., & Klimchuk, J. A. 1999, ApJ, 510, 485Arber, T. D., Botha, G. J. J., & Brady, C. S. 2009, ApJ, 705, 1183Arber, T. D., Haynes, M., & Leake, J. E. 2007, ApJ, 666, 541Arber, T. D., Longbottom, A. W., Gerrard, C. L., & Milne, A. M. 2001, Journal ofComputational Physics, 171, 151Archontis, V. 2008, Journal of Geophysical Research (Space Physics), 113, 3Archontis, V., Moreno-Insertis, F., Galsgaard, K., Hood, A., & O’Shea, E. 2004,A&A, 426, 1047Archontis, V. & T¨or¨ok, T. 2008, A&A, 492, L35 24 –Arkhypov, O. V., Antonov, O. V., & Khodachenko, M. L. 2011, Sol. Phys., 270, 1Arkhypov, O. V., Antonov, O. V., & Khodachenko, M. L. 2012, Sol. Phys., , 244Athay, R. G. & Thomas, R. N. 1961, Physics of the solar chromosphere, (InterscienceMonographs and Texts in Physics and Astronomy, New York: IntersciencePublication, 1961)Braginskii, S. I. 1965, Reviews of Plasma Physics, 1, 205Brown, J. C. 1973, Sol. Phys., 29, 421Chen, P. F. & Shibata, K. 2000, ApJ, 545, 524Cheung, M. C. M., Rempel, M., Title, A. M., & Sch¨ussler, M. 2010, ApJ, 720, 233Cowling, T. 1957, Magnetohydrodynamics, (New York: Cambridge University Press)de Pontieu, B. 1999, A&A, 347, 696Fan, Y. 2001, ApJ, 554, L111Fan, Y. & Gibson, S. E. 2007, ApJ, 668, 1232Fang, F., Manchester, IV, W., Abbett, W. P., & van der Holst, B. 2012, ApJ, 745,37Feynman, J. & Martin, S. F. 1995, J. Geophys. Res., 100, 3355Forbes, T. G. 2000, J. Geophys. Res., 105, 23153Gary, G. A. 2001, Sol. Phys., 203, 71Gary, G. A. & Alexander, D. 1999, Sol. Phys., 186, 123Gilbert, H., Kilper, G., & Alexander, D. 2007, ApJ, 671, 978Gilbert, H. R., Hansteen, V. H., & Holzer, T. E. 2002, ApJ, 577, 464Hood, A. W., Archontis, V., Galsgaard, K., & Moreno-Insertis, F. 2009, A&A, 503,999Khodachenko, M. L., Arber, T. D., Rucker, H. O., & Hanslmeier, A. 2004, A&A,422, 1073Klimchuk, J. A. 2001, in Space Weather, ed. . H. S. P. Song, G. Siscoe, volumeGeophysical Monograph 125, (Washington: AGU), 142 25 –Leake, J. E. & Arber, T. D. 2006, A&A, 450, 805Leake, J. E., Arber, T. D., & Khodachenko, M. L. 2005, A&A, 442, 1091Leake, J. E., Linton, M. G., , & Antiochos, S. K. 2010, ApJ, 722, 550Leake, J. E., Linton, M. G., & Antiochos, S. K. 2010, ApJ, 722, 550Leenaarts, J., Carlsson, M., Hansteen, V., & Gudiksen, B. V. 2011, A&A, 530, A124Linton, M. G. & Moldwin, M. B. 2009, JGR Space Physics, 114Longcope, D. W., Fisher, G. H., & Arendt, S. 1996, ApJ, 464, 999MacTaggart, D. & Hood, A. W. 2009a, A&A, 508, 445MacTaggart, D. & Hood, A. W. 2009b, A&A, 507, 995Magara, T. 2001, ApJ, 549, 608Manchester, IV, W., Gombosi, T., DeZeeuw, D., & Fan, Y. 2004, ApJ, 610, 588Priest, E. R. 2003, Solar magnetohydrodynamics, pages 38–42, (Kluwer AcademicPublishers)Rempel, M. 2011, ApJ, 740, 15Saha, M. N. 1921, Royal Society of London Proceedings Series A, 99, 135Stix, M. 2004, The sun : an introduction, (Springer-Verlag)Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635Zirin, H. 1998, The Astrophysics of the Sun, (Cambridge University Press, 1998)
This preprint was prepared with the AAS L A TEX macros v5.2.
26 –Name Model Reduced Eq. of state CZ profile Tube radius B ( kG ) β tube Fluxmass µ m Eq.(12) (cid:16) ∂T∂y (cid:17) a × Mx0 FIP m f m p / m f m p - Eq.(20) 0.3 Mm 6 2 5.42 PIP Eq.(13) No 2nd term Eq.(20) 0.3 Mm 6 2 5.43 PIP Eq. 13) 2nd term Eq.(31) 0.3 Mm 10.9 2 9.84 FIP m f m p - Eq.(20) 1.5 Mm 15 4 340.5 PIP Eq.(13) 2nd term Eq.(31) 1.5 Mm 61.2 4 1380.6 FIP m f m p - Eq.(20) 1.5 Mm 0.49 40 11.7 PIP Eq.(13) 2nd term Eq.(31) 1.5 Mm 19.2 40 430.Table 1: Simulations used in this paper. Simulation 0 is used only to highlight thechoice of the reduced mass µ m and is not run beyond t = 0 s. Simulation 2 hasthe same background stratification as the model used in Leake & Arber (2006). FIPstands for ‘fully ionized plasma’, and PIP for ‘partially ionized plasma’ . 27 – −2 0 2 4 6Height (Mm)10 T e m pe r a t u r e ( K ) −2 0 2 4 6Height (Mm)10 −12 −10 −8 −6 −4 −2 D en s i t y ( k g / m ^ ) −2 0 2 4 6Height (Mm)0.60.81.01.21.4 R edu c ed m a ss / p r o t on m a ss −2 0 2 4 6Height (Mm)10 −2 P r e ss u r e ( P a ) Sim 0Sim 1Sim 3Solar ModelVALC Model (a) (b)(c) (d)
Fig. 1.— Background 1D stratification for three models (Simulations 0, 1, 3). Dottedline: Simulation 0, the fully ionized simulation with µ m = m i /
2. Solid line: Simula-tion 1, the fully ionized simulation with µ m = m i . Dot-dashed line: Simulation 3, thepartially ionized model. The red line is from a standard model of the solar convectionzone (Stix 2004). The red crosses are the 1D semi-empirical atmospheric model ofVernazza et al. (1981). Panel (a): Temperature. Panel (b): Gas density. Panel (c):Reduced mass normalized by the proton mass. Panel (d): Gas pressure. 28 – −2 −1 0 1 2Radius (Mm)110100 p / B a t y = y t ube −2 −1 0 1 2Radius (Mm)−0.10−0.08−0.06−0.04−0.020.00 D en s i t y pe r t u r ba t i on ρ / ρ a t y = y t ube (a) (b) Fig. 2.— Panel (a): Initial p /B at y = y tube for the initial magnetic flux tubeembedded in the convection zone for the three small scale simulations, Simulation 1(the fully ionized simulation), Simulation 2 (the partially ionized simulation withoutionization effects in the equation of state), and Simulation 3 (the complete partiallyionized model). Panel (b): Initial density perturbation ρ /ρ , as in Equation (34), inthe tube for the same three simulations. 29 – −2−1012−2−1012 −2 −1 0 1 2−2−1012 −2 −1 0 1 2 −12 −10 −8 −6 −4 −2 a) : Sim1 t= 874 s b) : Sim1 t= 1093 sc) : Sim2 t= 874 s d) : Sim2 t= 1093 se) : Sim3 t= 874 s f) : Sim3 t= 1093 s Y(Mm)X(Mm) log(Density x m /kg) Fig. 3.— Evolution of three small scale (ephemeral region size) flux tubes. The whitelines show the in-plane field (constant contours of A z ) and the color contours showthe log of density. The A z contours are at seven values regularly spaced between(and inclusive of) the extrema of A z and 0.1 of this extrema (the extrema contouris a single point located at the center of the tube). Top row: Simulation 1, thefully ionized simulation. Middle row: Simulation 2, the partially ionized simulationwithout ionization effects in the equation of state. Bottom row: Simulation 3, thecomplete partially ionized model. 30 – H e i gh t ( Mm ) −1.8 −1.6 −1.4 −1.2 −1.0 −0.8 −0.6Height (Mm)0.00.51.01.5 R i s e s peed ( k m / s ) (a) Flux tube edges Flux tube centers (b)
Sim1Sim2Sim3
Fig. 4.— Panel (a): Height of the center and edge of the flux tubes as a function oftime. The center is defined as the location of the minimum in the vector potential A z . The edge is defined as the intersection of A z = 0 . A z ) with the x = 0 line.Panel (b): Rise speed in the convection zone of the centers of the flux tubes for thesame three simulations. The curves for Simulations 1 and 2 lie on top of each otherin panel (b). 31 – DeStabilizing TermStabilizing Term (d) Sim 1: 1020 s (b) Sim 2: 947 s(e) Sim 2: 1020 s (c) Sim 3: 874 s(f) Sim 3: 995 s
Fig. 5.— Line plots along x = 0 for the stabilizing ( − γ βδ ) and destabilizing( − H p ∂∂y ln ( B )) terms in the magnetic buoyancy instability inequality (35). Left: Sim-ulation 1, the fully ionized simulation. Center: Simulation 2, the partially ionizedsimulation without ionization effects in the equation of state. Right: Simulation 3,the complete partially ionized model. 32 – −2024681012−2024681012 −6 −4 −2 0 2 4 6−2024681012 −6 −4 −2 0 2 4 6 −12 −10 −8 −6 −4 −2 a) : Sim1 t= 1457 s b) : Sim1 t= 1967 sc) : Sim2 t= 1457 s d) : Sim2 t= 1967 se) : Sim3 t= 1457 s f) : Sim3 t= 1967 s Y(Mm)X(Mm) log(Density x m /kg) Fig. 6.— Same as Figure 3 but at later times of 1457 s and 1967 s, showing laterstage emergence of the flux tube in Simulations 1, 2 and 3. 33 – M agne t i c R e y no l d s N u m be r Sim2Sim3
Fig. 7.— Magnetic Reynolds number (using the Pedersen resistivity) as a function ofheight at x = 0 for the two partially ionized simulations. The plot is taken at time1359 s for Simulation 2, and at time 1724 s for Simulation 3. 34 –
500 1000 1500 2000 2500 3000Time (s)0.00.10.20.30.40.50.6 I n − p l ane f l u x abo v e . Mm Sim 1Sim 2Sim 3
Fig. 8.— In-plane flux above 1.2 Mm as a function of time, normalized to the in-planeflux in the initial flux tube, for Simulations 1, 2 and 3. 35 – T e m pe r a t u r e ( K ) Sim 1: t= 1991 sSim 2: t= 1359 sSim 3: t= 1675 st=0 s
Fig. 9.— Temperature as a function of height along the line x = 0. The timesin each simulation are chosen so that the transition regions in each simulation areapproximately co-spatial in height. 36 – × 10 −6 −5 −5 −5 N o r m a li z ed t o t a l c hange i n m a ss E _ s hea r (a) (b) Sim1Sim2Sim3
Fig. 10.— Panel (a): Change in mass above 1.2 Mm, normalized to the total mass inthe initial flux tube, as function of time. Panel (b): E shear , from Equation (38) forthe same three simulations. (a) (b) (c) Fig. 11.— The Lorentz force normalized by current density and magnetic field mag-nitudes | j ∧ B || j || B | , for Simulation 1 at 2377 s (panel a), Simulation 2 at 1560 s (panel b),and Simulation 3 at 1920 s (panel c). The times are shown so that in each figure the0 . A z ) contour has reached approximately the same height. 37 – −2 −1 Radius (Mm)110100 p / B a t y = y t ube −0.10−0.08−0.06−0.04−0.020.00 D en s i t y pe r t u r ba t i on ρ / ρ a t y = y t ube (a) (b) Sim 4Sim 5Sim 6Sim 7
0 1 2 −2 −1Radius (Mm) 0 1 2
Fig. 12.— Panel (a): Initial p /B at y = y tube for the large scale flux tube simulations.Solid line is Simulation 4, the β = 4 fully ionized model. Dot-dash line is Simulation5, the β = 4 partially ionized model. Double-dot dash line is Simulation 6, the fullyionized model with β ≈
40. Long-dashed line is Simulation 7, the partially ionizedmodel with β ≈
40. Panel (b): Initial density perturbation in the tube for the samethree simulations. 38 – −10−505101520 −15 −10 −5 0 5 10 15−10−505101520 −15 −10 −5 0 5 10 15 −12 −10 −8 −6 −4 −2 0 (a) : Sim4 t= 1578 s (b) : Sim4 t= 2307 s(c) : Sim5 t= 1093 s (d) : Sim5 t= 1821 s
Y(Mm)X(Mm) log(Density x m /kg) Fig. 13.— Emergence of large scale (active region size) flux tubes, showing in-planefield (constant contours of A z ) and color contours of the log of density. The A z contours are at seven values regularly spaced between (and inclusive of) the extremaof A z and 0.1 of this extrema (the extrema contour is a single point located at thecenter of the tube). Top row: Simulation 4, the fully ionized simulation with β = 4.Bottom row: Simulation 5, the partially ionized simulation with β = 4. 39 – −10−505101520 −15 −10 −5 0 5 10 15−10−505101520 −15 −10 −5 0 5 10 15 −12 −10 −8 −6 −4 −2 0 (a) : Sim6 t= 4615 s (b) : Sim6 t= 5829 s(c) : Sim7 t= 3594 s (d) : Sim7 t= 4687 s Y(Mm)X(Mm) log(Density x m /kg) Fig. 14.— Emergence of large scale (active region size) flux tubes, showing in-planefield (constant contours of A z ) and color contours of the log of density. Top row:Simulation 6, the fully ionized simulation with β = 40. Bottom row: Simulation 7,the partially ionized simulation β = 40. 40 – Time (s)−10010203040 H e i gh t ( Mm ) −5 −4 −3 −2 −1Height (Mm)0.00.51.01.52.02.5 R i s e s peed ( k m / s ) (a) Flux tubeedgesFlux tube centers (b)
Sim 4Sim 5 Sim 6Sim 7
Fig. 15.— Panel (a): Height of the center and edges of the flux tube as function oftime. The curves for the tube centers in Simulations 4 and 5 lie on top of each other.Panel (b): Axial rise speeds in the convection zone for the same four simulations. 41 –
Time (s)0.00.20.40.60.8 I n − p l ane f l u x abo v e . Mm Sim 4Sim 5Sim 6Sim 7
Fig. 16.— The amount of in-plane flux above 1.2 Mm as a function of time, normalizedto the integral of A z below y = 0 at t = 0 s, for Simulations 4, 5, 6 and 7. 42 – Height (Mm)10 T e m pe r a t u r e ( K ) Sim 4: t= 2185 sSim 5: t= 1772 sSim 6: t= 5439 sSim 7: t= 3933 st=0 s
Fig. 17.— Temperature as a function of height at x = 0 for all four large scalesimulations. The times in each simulation are chosen so that the transition regionsin each simulation are approximately co-spatial in the y-direction. 43 – -9 -8 -7 -6 -5 N o r m a li z ed t o t a l c hange i n m a ss E _ s hea r Sim 4Sim 5Sim 6Sim 7(a) (b)
Fig. 18.— Panel (a): Change in mass above 1.2 Mm, normalized to the total massin the initial flux tube, as a function of time. Solid line: Simulation 4, the fullyionized simulation with β = 4. Dot–dashed line: Simulation 5, the partially ionizedsimulation with the same β as Simulation 4. Double-dot-dashed line: Simulation 6,the fully ionized model with β = 40. Long dashed line: Simulation 7, the partiallyionized simulation with the same β as Simulation 6. Panel (b): E shearshear