• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Numerical Studies of Magnetic Reconnection and Heating Mechanisms for the Ellerman Bomb

    2023-05-26 08:31:22MingyuLiuLeiNiGuanChongChengUdoZieglerandJunLin

    Mingyu Liu, Lei Ni,3, Guan-Chong Cheng, Udo Ziegler, and Jun Lin,3

    1 Yunnan Observatories, Chinese Academy of Sciences, Kunming 650216, China; leini@ynao.ac.cn

    2 University of Chinese Academy of Sciences, Beijing 100049, China

    3 Center for Astronomical Mega-Science, Chinese Academy of Sciences, Beijing 100101, China

    4 Astrophysikalisches Institut Potsdam D-14482 Potsdam, Germany

    Abstract An Ellerman Bomb (EB) is a kind of small scale reconnection event, which is ubiquitously formed in the upper photosphere or the lower chromosphere.The low temperature(<10,000 K)and high density(~1019–1022)plasma there makes the magnetic reconnection process strongly influenced by partially ionized effects and radiative cooling.This work studies the high β magnetic reconnection near the solar temperature minimum region based on high-resolution 2.5D magnetohydrodynamics simulations.The time-dependent ionization degree of hydrogen and helium are included to realize more realistic diffusivities,viscosity and radiative cooling in simulations.Numerical results show that the reconnection rate is smaller than 0.01 and decreases with time during the early quasi-steady stage,then sharply increases to a value above 0.05 in the later stage as the plasmoid instability takes place.Both the large value of ηen (magnetic diffusion caused by the electron-neutral collision) and the plasmoid instability contribute to the fast magnetic reconnection in the EB-like event.The interactions and coalescence of plasmoids strongly enhance the local compression heating effect,which becomes the dominant mechanism for heating in EBs after plasmoid instability appears.However, the Joule heating contributed by ηen can play a major role to heat plasmas when the magnetic reconnection in EBs is during the quasi-steady stage with smaller temperature increases.The results also show that the radiative cooling effect suppresses the temperature increase to a reasonable range, and increases the reconnection rate and generation of thermal energy.

    Key words: magnetic reconnection – magnetohydrodynamics (MHD) – Sun: chromosphere – Sun: activity

    1.Introduction

    The interaction among the neutrals and ionized plasmas that affects magnetic reconnection is one of the major scientific challenges in understanding explosive phenomena in magnetized plasmas in the universe (e.g., Ni et al.2020; Ji et al.2022).The relatively cold and dense environment in the low solar atmosphere provides a natural laboratory that helps us to understand magnetic reconnection in partially ionized plasmas(Wang 1995).As the observational techniques develop, a growing number of small scale reconnection events have been observed,such as Ellerman Bombs(EBs)(e.g.,Ellerman 1917;Ding et al.1998; Georgoulis et al.2002), ultraviolet bursts(e.g., Peter et al.2014; Tian et al.2016; Chen et al.2019),spicules and chromospheric jets (e.g., dePontieu et al.2007;Samanta et al.2019; Shen 2021; Schmieder et al.2022),nanoflares and campfires(e.g.,Parnell&Jupp 2000;Song et al.2020; Berghmans et al.2021; Chen et al.2021; Joshi et al.2021), and some other explosive events (e.g., Huang et al.2019; Rast et al.2021; Xue et al.2021).

    Ellerman (1917) discovered an EB for the first time by analyzing Hα spectral line profiles in an active spot-group,where a sudden brightening in a very narrow band extending a few angstroms on either side of the spectral line was observed.EBs usually show elongated structures at Hα-0.8 ? (e.g.,Georgoulis et al.2002;Pariat et al.2007),and can possibly be composed of fine structures (Hashimoto et al.2010).Plenty of observations affirm that they also have strong emissions in Ca II lines (e.g., Fang et al.2006; Pariat et al.2007; Vissers et al.2013, 2015).Recent high-quality observations in Hβ indicate that EBs are ubiquitous in the quiet Sun and appear everywhere(Joshi & Rouppe van der Voort 2022).The size of an EB is about 1″ and the typical lifetime is a few minutes (e.g.,Georgoulis et al.2002;Pariat et al.2007).Statistical studies of their energy spectra show that the released energy in an EB is between 1025and 1028erg (e.g., Georgoulis et al.2002; Fang et al.2006).The temperature increase in EBs is about 600–3000 K according to the results based on semi-empirical models and numerous spectral analyses (e.g., Georgoulis et al.2002; Fang et al.2006; Hong et al.2017).

    Magnetic reconnection is believed to be the main mechanism that accounts for an EB.When EBs happen,the strength of the measured magnetic fields at the solar surface of the corresponding reconnection region is about a few hundred Gauss (e.g.,Pariat et al.2004,2007;Li et al.2015).However,the reconnection diffusion region in EBs is still poorly understood because of the limited resolutions of the existing solar telescopes.The multi-wavelength spectral studies and radiative hydrodynamic simulations demonstrate that EBs always occur in the upper photosphere or the lower chromosphere(e.g.,Hong et al.2017,2021)where the plasma density is around 1019–1022cm-3and it is weakly ionized according to the C7 solar atmosphere model (Avrett & Loeser 2008).Therefore, the strong radiative cooling and partially ionized effects have to be considered when we study magnetic reconnection in this region.

    Previous single-fluid two-dimensional (2D) magnetohydrodynamics(MHD)simulations have studied magnetic reconnection between anti-parallel magnetic fields in the low solar atmosphere (Chen et al.2001; Xu et al.2011), where the temperature increase and duplicated Hα and Ca II line profiles can be qualitatively compared with observations.Magnetic reconnections between the U-shaped magnetic fields triggered by Parker instability have also been shown to occur in both 2D and three-dimensional (3D) MHD simulations (Isobe et al.2007; Archontis & Hood 2009), which are consistent with the scenario that EBs are usually observed in the regions with serpentine magnetic field lines (Pariat et al.2004).Recent 3D radiation MHD simulations studied magnetic reconnection between the emerging magnetic field and the background magnetic field below the temperature minimum region (TMR)(Danilovic et al.2017; Danilovic 2017), and the synthesized images in the Hα wing seen from different angles are similar to the flame-like structures observed, and the synthetic magnetograms obtained from Fe I 6301 ? are also consistent with observations.However, we should point out that all these previous MHD simulations about modeling EBs in the stratified solar atmosphere usually focus on comparisons with observations, and the heating mechanisms and small scale physics in the reconnection region were not examined.

    The theoretical and numerical studies about magnetic reconnection in partially ionized plasmas indicate that the ambipolar diffusion caused by the decoupling of ions and neutrals(Brandenburg&Zweibel 1994),the ion recombination effect (Leake et al.2012, 2013) and radiative cooling(Uzdensky et al.2010) can accelerate reconnection when the guide field vanishes.Both single-fluid and multi-fluid MHD simulations show that turbulent reconnection mediated by plasmoids also appears in the low solar atmosphere (Leake et al.2012,2013;Murphy&Lukin 2015;Ni et al.2015,2016;Ni&Lukin 2018;Singh et al.2019;Murtas et al.2021,2022),and plasmoid instability can lead to fast magnetic reconnection(Ni et al.2015;Ni&Lukin 2018),especially in the case of low β (Ni &Lukin 2018).Which of these mechanisms govern fast magnetic reconnection in an EB is still an open question.Previous simulations also indicate that Joule heating is very important for converting magnetic energy into heat (Ni et al.2016)in such a dense plasma.The heating caused by ambipolar diffusion can be ignored in a low β reconnection process around the solar TMR(Ni et al.2021,2022),but it might play a role in a high β reconnection event such as an EB (Chae &Litvinenko 2021).When plasmoid instability appears in the reconnection current sheet, many small scale slow-mode and fast-mode shocks appear inside the plasmoid (Ni et al.2015, 2016), which can also heat plasmas in the reconnection region.The strong horizontal flows at the solar surface are usually observed during the EB formation process.Recent 3D MHD simulations demonstrate that these horizontal flows caused by external forces can drive the reconnection process and cause compression heating in EBs (Cheng et al.2021).Nevertheless, the heating mechanism in EBs still needs more precise simulations to verify.

    In this work, we study reconnection and the heating mechanism in EBs based on single-fluid MHD simulations.The MHD models have been improved by including both hydrogen and helium, and the time-dependent ionization degree of these two elements, which makes the magnetic diffusion caused by electron-neutral collision, the ambipolar diffusion and the radiative cooling effect more realistic than those in previous simulations (Ni et al.2015, 2016).The numerical methods are described in Section 2.We present the numerical results in Section 3.Discussions about fast reconnection and heating mechanisms in EB are provided in Section 4.We give our conclusions and outlooks in Section 5.

    2.Numerical Methods

    2.1.MHD Equations

    In this work, we use the optimized single-fluid MHD code NIRVANA (Ziegler 2008, 2011) to perform the 2.5D MHD simulations.All the five components including hydrogen atoms, helium atoms, electrons, and ions contributed by ionized hydrogen and helium are considered to be well coupled as a single fluid.The solved MHD equations are the same as those in the recent work by Ni et al.(2022):

    Table 1 Collision Cross Sections, the Unit is m2

    The variables ρ,v,B,p,T,e,YiHand YiHeare the mass density,fluid velocity, magnetic field, thermal pressure, temperature,total energy density, and ionization fractions of hydrogen and helium, respectively.EADis the ambipolar diffusion electric field, τSis the stress tensor, and Qradand H are the radiative cooling and heating terms respectively, which will be more detailedly described in the next subsections.The magnetic diffusivity η is contributed by both electron-ion collision (ηei)and electron-neutral collision (ηen).The constant values of mi=1.66057×10-27kg, kB=1.3806×10-23J K-1,μ0=4π×10-7V s A-1m-1and γ=5/3 are the mass of a proton, Boltzmann constant, magnetic permeability coefficient and ratio of specific heats, respectively.The International System of Units (SI) are applied in the simulations.

    We study magnetic reconnection around the solar TMR,where the ion-neutral collision frequency is about 108s-1, and the ion inertial length and ion-neutral collision mean free path are about 1 m and 100 m (e.g., Leake et al.2013; Ni &Lukin 2018), respectively.Therefore, it is reasonable to use such a single-fluid model when we study an EB-like event with a length scale of 200 km.However, we should also point out that the single-fluid model will eventually break down when the reconnection length scale decreases down to the ion-neutral collision mean free path.In this work,the effects of neutrals are embodied by the magnetic diffusion caused by electron-neutral collision and ambipolar diffusion.The number density of the total helium is assumed to be 10% of the total hydrogen.We only consider the primary ionization of helium, and it is reasonable in our low temperature simulations.Since the current sheet in our simulations is assumed to be parallel to the solar surface,we ignore gravity and the initial plasma density is taken as a constant.

    2.2.Diffusions and Viscosity

    The high plasma density in EBs causes strong collisions between different kinds of particles.In this subsection,we will introduce the magnetic diffusion, ambipolar diffusion and viscosity applied in this work.The magnetic diffusion is contributed by both electron-neutral collision and electron-ion collision, while the ambipolar diffusion mainly relates to the ion-neutral collisions, and the viscosity is mainly caused by ion–ion collisions and neutral–neutral collisions.The different cross sections used to calculate the diffusions and viscosity are listed in Table 1; their values are given according to the previous work(Barata&Conde 2010;Vranjes&Krstic 2013),here σe-nH, σe-nHe, σiH-nH, σiH-nHe, σiHe-nH, σiHe-nHeand σnH-nHare the collision cross sections of electron-neutral hydrogen collision, electron-neutral helium collision, ion hydrogen-neutral hydrogen collision, ion hydrogen-neutral helium collision, ion helium-neutral hydrogen collision, ion helium-neutral helium collision and neutral hydrogen-neutral hydrogen collision, respectively.

    The magnetic diffusivity η contributed by two different kinds of collisions is expressed as

    where ne=ρ(YiH+0.1YiHe)/(1.4mi) is the number density of electrons, expressed in cgs units (cm-3), and T is expressed in eV.

    The ambipolar diffusion electric field EADin the energy equation,Equation(3),and induction equation,Equation(4),is given by

    where ηADis the ambipolar diffusion coefficient.The formula of ηADis as below (Ni et al.2020, 2022)

    where the unit of ηADis m3s kg-1.Since the plasmas are composed of hydrogen and helium, we can obtain:

    where ρiH=ρYiH/1.4 is the ionized hydrogen density and ρiHe=0.4ρYiHe/1.4 is the ionized helium density.According to the previous inference (Ni et al.2022), the electron collision part ρeνencan be ignored.

    The stress tensor in the momentum equation, Equation (2),and the energy equation, Equation (3), is given as

    where ξ is the dynamic viscosity coefficient and its unit is kg m-1s-1.The expression for the dynamic viscosity can be simply written as below (Ni et al.2022)

    The detailed derivation process of these diffusions and viscosity can be found in the recent work about magnetic reconnection in ultraviolet bursts (Ni et al.2022).The ionization degrees of hydrogen and helium are based on the RADYN test atmosphere results by solving the radiative transfer equations (Carlsson & Leenaarts 2012), which vary with plasma temperature as plotted in Figure 1(a).The ionization degree of helium is a simple exponential function of temperature, YiHe=1-100.325571-0.0000596T.According to this expression, the value of YiHewill become negative if the temperature is lower than 5413 K.Therefore, we assume YiHe=0.00010084814 when T <5413 K.The curve of YiHe(as displayed in Figure 1) is a good fit down to 50% ionization at T=104K.

    2.3.Radiative Cooling Models

    The radiative cooling effect is very strong in an EB, which cannot be ignored when we study magnetic reconnection in this region.As we know, the photosphere can be considered to be in local thermodynamic equilibrium (LTE) and the chromosphere is in a non-LTE state,which is much more complicated when we try to solve the radiative transfer equations in this region (Rutten 2003).Gan & Fang (1990) derived a simple radiative cooling expression according to the observational experiences, which can roughly measure the radiative cooling effect in the low solar atmosphere.The chromospheric radiative energy balance is dominated by a small number of strong lines from neutral hydrogen, singly ionized calcium and singly ionized magnesium (Vernazza et al.1981).Carlsson &Leenaarts (2012) provided a simple radiative cooling model based on the three lines, which is considered as the most accurate simple radiative cooling for the chromosphere,especially for the upper chromosphere.The reconnection region we study is around the solar TMR.We cannot tell which model is more accurate for this region.Therefore, we simulate different cases with the two different radiative cooling models and compare the results.

    As clearly described in the recent work Ni et al.(2022), the first radiative cooling model (Carlsson & Leenaarts 2012) is written as below

    Figure 1.(a)Temperature-dependent ionization fractions of hydrogen(YiH)and helium(YiHe).(b)Field lines and current density Jz(background color)in the whole simulation domain for Case III at t=0 s.

    2.4.Initial Setups

    The simulation length scale L0=200 km, which is comparable with the length scale of an EB.The 2D simulation domain is 0 ≤x ≤L0and -0.5L0≤y ≤0.5L0, and open boundary conditions are applied in both directions.According to the solar atmosphere model, we assume the initial temperature T0=4400 K and the initial total plasma density ρ0=1.66057×10-6kg m-3, which are the typical values around TMR.The horizontal force-free Harris sheet is applied as the initial equilibrium magnetic fields:

    The small perturbations for both magnetic field and velocity at the beginning are given as below:

    where bpertand vpertare small fractions of b0and vA0respectively, and randomnis the random noise function in the simulation code.

    In order to study the effects of radiative cooling, different diffusions, initial perturbations and resolutions on magnetic reconnection around the solar TMR, we have tested many different cases.In this paper,we will show the results of seven cases, and the differences among these cases are listed in Table 2.The adaptive mesh refinement(AMR)level means the highest level of mesh refinement applied in the simulations.We can see that the only difference between Cases I, II and III is the radiative cooling model.We will compare the results from these three cases to study the effect of radiative cooling on magnetic reconnection.The only difference between Cases III,IV and V is the different diffusion terms included in these cases.The resolution in Case VI is much lower than that in the other cases, which makes the numerical diffusion dominateover the physical diffusions in this case.Compared with the other cases, a much smaller initial perturbation of magnetic field is applied in Case VII.

    Table 2 Differences Among Seven Simulation Cases

    3.Numerical Results

    3.1.Results with Different Radiative Cooling Models

    In order to investigate the effects of radiative cooling on the magnetic reconnection process in EBs, we will detailedly describe and compare the results in Cases I, II and III in this subsection.The radiative cooling term is turned off in Case I,while Qrad2and Qrad1are separately included in the energy equation in Case II and Case III.

    The evolutions of the current sheet in different cases are shown in Figure 2.The 2D distributions of the temperature and plasma density at three different times are presented in this figure, in which the two different background colors in half of the upper and lower panels represent the temperature and density, respectively.The time evolutions of the maximum temperature Tmax, the maximum velocity in the x-direction vxmaxand the reconnection rate γ in the three cases are presented in Figure 3.The reconnection rate is calculated as:γ=vy-aver/vA-aver,where vy-averis the average velocity in the y-direction and vA-averis the average Alfvén velocity inside a small region around the main reconnection X-point; the red boxes in Figure 2 are these small regions for calculating the reconnection rate.

    The initial Harris sheet always becomes thinner and thinner,and then breaks down to form the plasmoid instability at around t=37 s in all the three cases.These plasmoids collide and coalesce to form bigger and bigger ones.As discussed in the previous papers Ni et al.(2016,2022),many slow and fast mode shocks are formed inside the plasmoids during the collision process, which can also be identified in Figure 2.As demonstrated in Figures 2 and 3, the evolutions of all the variables are very similar in the three cases before plasmoid instability appears.

    After that, the radiative cooling strongly affects the temperature and density distributions.The temperature sharply increases with time from 4400 K to above 9400 K in Case I without the radiative cooling effect.Such a temperature increase is too large to match the measured temperature increase in a typical EB(~600–3000 K)(e.g.,Georgoulis et al.2002; Fang et al.2006; Hong et al.2017).The temperature increase becomes slower and slower in Case II and Case III after plasmoid instability appears, and the maximum temperature stays at about 7000 K in Case II and 6800 K in Case III,eventually.Therefore, the temperature increases in Case II and Case III are in good agreement with that in a typical EB.The slightly higher temperature in Case II indicates that the Carlsson & Leenaarts (2012) radiative cooling model is stronger than the Gan & Fang (1990) model.Figure 2 also shows that the strong radiative cooling results in the more nonuniform density distributions inside the plasmoids, and the maximum density in Case III is much higher than that in Case I at t=70.32 s.In Figure 3(b) and (c), we can also see that the radiative cooling efficiently accelerates the outflow velocity and increases the reconnection rate to a higher value during the later stage.The maximum outflow velocity in Case III stays around 9 km s-1,which is consistent with the inverted outflow velocity of a chromospheric reconnection event based on the observational results (Díaz Baso et al.2021).

    As we know,the thermal energy equation in the MHD model is stated as below:

    Figure 2.The 2D distributions of temperature and plasma density at three typical times in six different cases.The two different background colors in half of the upper and lower panels correspond to the temperature and density, respectively.The black solid lines represent the magnetic field lines.

    where ethis the thermal energy density and the “tr”in the second term on the right-hand side means evaluating the trace of a matrix.Figure 4 shows the evolutions of the average power densities contributed by different terms on the right-hand side of Equation (25) in different cases.Variables Qcomp, Qvis, Qei,Qen, QAmpand Qradare the average power densities of the compression heating, heating caused by the viscosity ξ, Joule heating contributed by magnetic diffusions ηeiand ηen,heating contributed by the ambipolar diffusion ηADand lost thermal energy by the radiative cooling effect, respectively.They are the average values calculated in the area of the main reconnection region within -0.026L0≤y ≤0.026L0and 0 ≤x ≤L0.We should point out that the compression heating Qcomphas been calculated by using the average values of both v·?pand-p?·vinthemainreconnectionregion,because·vdFdt(Birn et al.2012), where ∫dF represents the integration in the x and y directions.

    Figure 3.The time evolutions of the maximum temperature Tmax(a),the maximum outflow velocity in the x-direction vxmax(b),and the reconnection rate γ in Cases I,II and III.

    Figure 4.The time evolutions of the average power density contributed by different heating terms and the radiative cooling in six different cases.

    Figure 5.The time evolutions of the different diffusion coefficients (ηei, ηen and ηAmp) at the main reconnection X-point in Case III and Case VII.

    Comparing the results in Figure 4(a), (b) and (c), one can find that the evolutions of different heating terms are very similar in Cases I,II and III.Qenis larger than the other heating terms during the early reconnection stage.However, Qensharply decreases to a much lower value after plasmoid instability appears,but the general trends of all the other terms increase with time.Eventually,Qcompdominates over the other terms, Qendecreases to a value with the same order of magnitude as Qeiand QAmp, and Qvisis always much smaller than the other terms and can be ignored in all the cases.When we look into the details,we can find the differences among the three cases.The value of Qcompin Case III reaches above 10 erg cm-3s-1, while this value is smaller than 10 erg cm-3s-1in Case I.The lower temperature in Case III also makes the values of Qen, QAmpand Qeiin Case III be higher than those in Case I.Therefore, the radiative cooling effect efficiently promotes the generation of more thermal energy.As demonstrated in Figure 4(b) and (c), the radiative cooling effect is much weaker than total heating effect before plasmoid instability appears in both Case II and Case III.After that,the radiative cooling term Qradsharply increases,and then it gradually becomes stable after t=60 s in both cases.The radiative cooling effect has been increased by four orders of magnitude during the whole reconnection process of 100 s.The final average power density of the radiative cooling term Qradis also ~10 erg cm-3s-1in both cases,and the radiative losses at the reconnection site are then calculated as Lrad=Qrad×L0=200 kW m-2, which is very close to the calculated radiative loss (160 kW m-2) of the chromospheric reconnection event based on the observational results (Díaz Baso et al.2021).

    3.2.Results with Different Diffusions

    The collisions between different kinds of particles in EBs are also strong because of the high plasma density there.Such strong collisions might efficiently affect the magnetic reconnection process through the different diffusion terms in the energy equation, Equation (3), and the magnetic induction equation, Equation (4).Figure 5 displays the time-dependent diffusivities at the main reconnection X-point inside the current sheet.In order to unify the units of three diffusion coefficients,we evaluate the ambipolar diffusion via ηAmp=B2ηAD/μ0.As shown in Figure 5, the magnetic diffusivity contributed by electron-neutral collision ηenand the ambipolar diffusion coefficient ηAmpare above 105m2s-1at the beginning, and they are almost two orders of magnitude higher than the magnetic diffusivity contributed by electron-ion collision ηei.However, ηenand ηAmpat the main reconnection X-point sharply decrease during the reconnection process because of the increasing temperature and ionization degree;they decrease to a value that is smaller than ηeiduring the later stage.Therefore, it is necessary to detailedly check how these time dependent diffusivities affect the whole reconnection process.The results in Cases III, IV, V and VI are compared to check the effects of different diffusions.We only include the ηeiterm in Case IV,but we include both the ηeiand ηenterms in Case V.The terms related to the ambipolar diffusion ηADare turned off in both Case IV and Case V.Though all the diffusion terms are included in Case VI, the limited resolution causes the numerical diffusion in Case VI to be much larger than the three physical diffusions.

    The 2D distributions of the temperature and plasma density at three different times in Cases III,IV,V and VI are presented in Figure 2(c),(d),(e)and(f),respectively.The time evolutions of the maximum temperature Tmax,maximum velocity in the xdirection vxmaxand reconnection rate γ in the four cases are presented in Figure 6.Comparing Figure 2(c) and (e), we can find that the missing ambipolar diffusion effect slightly changes the evolutions of the plasmoid instability;the plasmoids appear a little bit later and they are slightly smaller in Case V than those in Case III at the same point in time.As shown in Figure 6, the time evolutions of Tmax, vxmaxand γ are very similar in Case III and Case V.These results indicate that the ambipolar diffusion is not important in such a reconnection process.Comparing Figure 2(d) and (e), one can clearly see that the missing ηenterms cause strong influences on the reconnection process.Before t=60 s, the temperature inside the current sheet in Case IV is much lower than that in Case III,and the plasmoids appear about 5 s later in Case IV.The missing ηenterms in Case IV also result in the smaller outflow velocity vxmaxand smaller reconnection rate γ as displayed in Figure 6(b) and (c).These results indicate that the magnetic diffusion caused by the electron-neutral collision plays an important role in the magnetic reconnection process of an EB.

    Figure 6.The time evolutions of the maximum temperature Tmax(a),the maximum outflow velocity in the x-direction vxmax(b),and the reconnection rate γ in Cases III, IV, V and VI.

    Figures 2(f)and 6 show that the large numerical diffusion in Case VI makes the evolution of the whole reconnection process be very different.The plasmoids appear much later and much fewer plasmoids appear in this case.The temperature increase and reconnection outflow velocity are also strongly underestimated in Case VI.Therefore, the realistic physical diffusions and high resolutions in numerical simulations are very necessary for investigating the reconnection and heating mechanisms in EBs.

    4.Discussion

    4.1.Mechanisms Leading to Fast Magnetic Reconnection

    The high plasma density in EBs makes the collisions between different kinds of particles and the radiative cooling effect be strong in the magnetic reconnection region.The strong collisions will lead to the large diffusion terms in the energy equation, Equation (3), and the magnetic induction equation, Equation (4).The strong radiative cooling will also affect the evolution of the energy equation, Equation (3).

    Previous theoretical results (Brandenburg & Zweibel 1994)demonstrated that the ambipolar diffusion caused by the decoupling of ions and neutrals can result in a thinner current sheet and much higher reconnection rate when there is no guide field.Theoretical work in Uzdensky et al.(2010)indicated that the strong radiative cooling can also accelerate the reconnection when the guide field vanishes.The 2D numerical results in Ni et al.(2015)then demonstrated that the ambipolar diffusion and radiative cooling indeed make the thinning of the current sheet faster and the reconnection rate higher before plasmoid instability appears, in the case with zero guide field.However,the reconnection rate always sharply increased to a value~0.02–0.03 after plasmoid instability appears (Ni et al.2015),no matter whether the ambipolar diffusion and radiative cooling effects are included in the simulations or not.

    Figures 3(c) and 6(c) show the time evolutions of the reconnection rate in different cases.We can see that the reconnection rate is smaller than 0.01 before plasmoid instability appears.After t=25 s, the reconnection rate decreases to an even lower value.The reason is that the reconnection rate γ scales with the Lundquist number S as γ ~S-1/2in the earlier quasi-steady stage with the Sweet-Parker like current sheet,where S=LvA/η.The decreasing η as shown in Figure 5 causes the decreasing reconnection rate during this stage.After plasmoid instability appears, the reconnection rate is always sharply increased in all the cases,which is consistent with the previous simulations (Ni et al.2015;Ni &Lukin 2018).However, there are some differences between this work and the previous one (Ni et al.2015).Comparing the results in Case I and Case III, we can see that the strong radiative cooling effect still slightly increases the reconnection rate after plasmoid instability appears.The results in Figure 6(c)demonstrate that ηenalso plays an important role that results in the fast magnetic reconnection; the reconnection rates in Case III and Case V are about two times larger than that in Case IV,even after plasmoid instability appears.Since ηenis about 25 times larger than ηeibefore t=30 s in Case III, the reconnection rate in Case III is expected to be about 5 times higher than that in Case IV without ηenduring this stage.However, the reconnection rate in Case III is only slightly higher than that in Case IV as shown in Figure 6(c).The reason is that the AMR level is smaller than 4 before t=30 s, which causes the low resolution during this stage.The low resolution then results in the high numerical diffusion, and the reconnection rate in Case IV is overestimated during this stage.The similar results in Case III and Case V demonstrate that the ambipolar diffusion effect can be ignored in the whole reconnection process.

    In order to confirm these results and conclusions, we have calculated the reconnection rate in larger boxes in the upstream region around the main reconnection X-point.The derived reconnection rate in the larger box is relatively lower but the whole evolution trend is the same as those shown in Figures 3(c) and 6(c).The reconnection rate at the main-X point is also calculated based on the values of ηJz(Leake et al.2012, 2013; Ni & Lukin 2018), which exhibits larger fluctuations than the results presented in Figures 3(c) and 6(c) after plasmoid instability appears.However, the main conclusions are still not changed.

    Our numerical results indicate that both the large magnetic diffusion caused by electron-neutral collisions and the plasmoid instability can lead to fast magnetic reconnection in an EB-like event.The strong radiative cooling in EBs can also accelerate magnetic reconnection even when the strong guide field is included.However, we should point out that the research about the ion recombination effect on magnetic reconnection (Leake et al.2012, 2013) in the EB-like event is outside the scope of this work.

    4.2.Heating Mechanisms

    Though magnetic reconnection is believed to be the main mechanism accounting for EBs, the energy converting mechanism inside the reconnection diffusion region is still not well understood.The high resolution MHD simulations with time dependent ionization degree, more realistic diffusions and radiative cooling model allow us to deeply analyze the heating mechanisms in EBs.

    Figure 4 plots the time evolutions of the average power density contributed by different heating terms in the thermal energy equation, Equation (25).The results in Case III show that the Joule heating contributed by electron-neutral collision(Qen) is larger than the other heating terms during the early reconnection stage.After plasmoid instability appears, Qensharply decreases to a much lower value,but the general trends of all the other terms increase with time.The compression heating Qcompexceeds the Joule heating Qenafter 25 s.Eventually, Qendecreases to ~0.4 erg cm-3s, Qeiand QAmpincrease to a value with the same order of magnitude as Qen,and Qcompincreases to a value above ~10 erg cm-3s.The numerical results in Case VII show similar behaviors as those in Case III.However,the five times smaller initial perturbation in Case VII causes the slower evolution of the current sheet and the plasmoids appear much later.In Case VII, Qenis much larger than Qcompand the Joule heating caused by ηenis the most important heating term before plasmoid instability appears.Then, the compression heating Qcompincreases and exceeds the Joule heating Qenafter plasmoid instability appears(~90 s).We can expect that Qcompwill also increase to a value with Q=10 erg cm-3s in Case VII during the later stage.When such a turbulent reconnection process happens in an EB with the typical size of V=700×700×700 km3and the lifetime of t=5 minutes, the total thermal energy supplied by magnetic reconnection is E=Q×V×t ?1027erg,which is in the energy range of a typical EB (~1025–1028erg) measured from observations (e.g., Georgoulis et al.2002; Fang et al.2006).

    The numerical results indicate that the sharp increase of the compression heating is triggered by the appearance of the plasmoids.The compression heating always dominates to supply the thermal energy in EBs as long as the turbulent reconnection mediated by plasmoids appears.As affirmed in Figure 2,the interactions and coalescence of the plasmoids can strongly enhance the local compression inside the plasmoids,many slow-mode and fast-mode shocks appear (Ni et al.2015, 2016, 2022), and the kinetic energy generated in the reconnection process is then converted into thermal energy in this process (Ni et al.2022).Figure 7(a) presents the divergence of the velocity in a plasmoid shown in Figure 2(c); the large compressions appear in the regions with obvious blue and red colors.The vertical white dash-dotted line in Figure 7(a)is the same line as drawn in Figure 2(c).The two fast-mode shock fronts are located at around the two heads of the plasmoid with blue colors.The vertical white dash-dotted line crosses a pair of slow-mode shock fronts with red colors.Figure 7(b) plots the distributions ofρ ρmax,T Tmaxandp pmaxalong the vertical white dash-dotted line.The results in Figure 7(b)demonstrate that the variables suddenly increase to higher values around the two slow-mode shock fronts.Figure 7(c) displays the time evolutions of Qcompand Qkin,where Qkinis the average power density of the kinetic energy calculated inside the reconnection region within-0.026L0≤x ≤0.026L0and 0 ≤y ≤L0.We can see that the sharp increase of Qcompalways corresponds to the sharp decrease of Qkin.Since the small shock structures inside the plasmoids strongly vary with time and location and the other compression heating processes always combine with shocks at the same time, it is difficult to evaluate how much thermal energy is converted from kinetic energy only via shocks.

    Figure 7.(a)The divergence of the velocity in a plasmoid shown in Figure 2(c);the vertical white dash-dotted line is the same line as drawn in Figure 2(c).(b)The distributions of ρ ρm ax , T Tmax and p pmax along the vertical white dash-dotted line in (a), where ρm ax, Tmax and pmax are the maximum values of mass density,temperature and pressure along this line respectively.(c) The time evolutions of Qcomp and Qkin, where Qcomp and Qkin are the average power densities of the compression heating and kinetic energy, respectively.

    Such a result is different from the previous work in Ni et al.(2016),which showed that Joule heating dominates to heat the plasmas in the magnetic reconnection process in EBs with a similar plasma β as in this work.However,we should point out that the fixed ionization degree caused magnetic diffusion contributed by electron-neutral collision and ambipolar diffusion to be overestimated,and the radiative cooling process was underestimated in the reconnection region in Ni et al.(2016).Therefore, the Joule heating and temperature increase were overestimated in Ni et al.(2016).The current model in this work has been improved significantly by including the time dependent ionization degree of hydrogen and helium, and a more realist radiative cooling model.In addition, the coalescence of the plasmoids lasts for a longer time to form much bigger plasmoids in this work, which lead to a larger compression heating effect.Recently, Ni et al.(2022) studied the low β magnetic reconnection process in ultraviolet bursts based on the same MHD model as in this work, and the main difference between the current model in this work and the one in Ni et al.(2022)is the magnetic field strength,and hence the plasma β.Their results demonstrated that the thermal energy accounting for ultraviolet bursts around the solar TMR is also mainly supplied by the local compression heating triggered by plasmoid instability.The much lower plasma β in Ni et al.(2022) leads to the much stronger compression heating effect.Therefore, we can conclude that such a mechanism plays an important role for different kinds of local heating events related to magnetic reconnection in the low solar atmosphere.

    However, we should point out that the plasmoid instability might not always appear in the EB-like events,or in the whole stage of the EB events.Then, the Joule heating contributed by ηen(the magnetic diffusion caused by electron-neutral collision) will play the most important role for supplying thermal energy in a lower temperature EB with a quasi-steady state reconnection process.

    5.Conclusions and Outlooks

    In this work, we have studied the reconnection and heating mechanisms in EBs based on a single-fluid 2.5D MHD framework.Compared with previous works (Ni et al.2015, 2016, 2021), the time dependent ionization degrees of hydrogen and helium are included to result in more realistic diffusivities and viscosity, and two improved radiative cooling models have been applied in the simulations.The reconnection current sheet is assumed to be located around the solar TMR with the initial temperature of 4400 K and density of 1.66057×10-6kg m-3.We have tested seven different cases to study the effects of radiative cooling, different diffusions,initial perturbations and resolutions on magnetic reconnection.The main conclusions are summarized as below:

    1.In the weakly ionized plasmas around the solar TMR,the reconnection rate is smaller than 0.01 and decreases with time during the quasi-steady state Sweet-Parker like current sheet stage, then sharply increases to a value higher than 0.05 during the later unstable reconnection stage with plasmoid instability.Both the large value of ηenand plasmoid instability contribute to the fast magnetic reconnection in an EB-like event.

    2.The compression heating always dominates to supply thermal energy in EBs as long as turbulent reconnection mediated by plasmoids appears.The local compression heating in the reconnection region was efficiently enhanced when the plasmoids collide and coalescence with each other, and part of the generated kinetic energy is converted to thermal energy during this process.However,the Joule heating contributed by ηencan play a major role to heat plasmas when magnetic reconnection in EBs is during the quasi-steady stage with smaller temperature increases.

    3.The average power density caused by the compression heating can reach above 10 erg cm-3s,which can supply a total energy of 1027erg for an EB with size of 1″ and lifetime of 5 min.The strong radiative cooling in EBs results in radiative loss at the reconnection site being 200 kW m-2, which is close to the measured radiative loss (160 kW m-2) of the observed chromospheric reconnection event (Díaz Baso et al.2021).

    4.The strong radiative cooling effect in the reconnection region constrains the temperature increase to a reasonable value (~2400 K), which is consistent with the results based on semi-empirical models and numerous spectral analyses from observations (e.g., Georgoulis et al.2002;Fang et al.2006;Hong et al.2017).The radiative cooling effect also promotes the generation of more thermal energy,accelerates the outflow velocity and increases the reconnection rate.

    Previous two-fluid MHD simulations showed that the ion recombination effect can result in fast magnetic reconnection in weakly ionized plasmas(Leake et al.2012,2013).Future multifluid MHD simulations with a more realistic radiative cooling model is still necessary to study the magnetic reconnection process in EBs.The recent 3D radiation magnetohydrodynamics (RMHD) simulations have successively simulated the formation of EBs in the photosphere and synthesized Hα images and spectral line profiles (Danilovic et al.2017;Danilovic 2017).However, these 3D RMHD simulations always showed a single X-point structure inside the reconnection region.The high resolution observational results indicate that EBs might include many fine structures (Hashimoto et al.2010; Díaz Baso et al.2021), and the recent 3D MHD simulations also demonstrated that a small scale flux rope(corresponding to the plasmoid in 2D) can be generated in the reconnection process of an EB (Cheng et al.2021).Therefore,future higher resolution 3D RMHD simulations with more realistic diffusions are very important to reveal the magnetic reconnection and heating mechanism, and the fine structures inside an EB.The higher resolution and more precise observational results from the Daniel K.Inouye Solar Telescope(DKIST) (Rast et al.2021), Chinese H-Alpha Solar Explorer(CHASE) (Li et al.2022) and other telescopes should be combined with simulations to reveal the hidden smaller scale physics in EBs.

    Acknowledgments

    This research is supported by the National Key R&D Program of China No.2022YFF0503800;the National Natural Science Foundation of China (NSFC, Grant Nos.11973083 and 11933009); the Strategic Priority Research Program of CAS with grants XDA17040507; the outstanding member of the Youth Innovation Promotion Association CAS (No.Y2021024); the Applied Basic Research of Yunnan Province in China Grant 2018FB009;the Yunling Talent Project for the Youth; the project of the Group for Innovation of Yunnan Province grant 2018HC023;the Yunling Scholar Project of the Yunnan Province and the Yunnan Province Scientist Workshop of Solar Physics;Yunnan Key Laboratory of Solar Physics and Space Exploration (No.202205AG070009).The numerical calculations and data analysis were done at the Hefei Advanced Computing Center and the Computational Solar Physics Laboratory of Yunnan Observatories.

    精品欧美国产一区二区三| 久久香蕉国产精品| 一a级毛片在线观看| www日本黄色视频网| 最好的美女福利视频网| 欧美+亚洲+日韩+国产| 日韩欧美精品v在线| 亚洲美女黄片视频| 日韩欧美 国产精品| 特级一级黄色大片| 亚洲 欧美一区二区三区| 女人被狂操c到高潮| 亚洲一卡2卡3卡4卡5卡精品中文| 中国美女看黄片| 国产精品 国内视频| 91av网站免费观看| 久久精品国产清高在天天线| 午夜a级毛片| 久久久国产欧美日韩av| 全区人妻精品视频| 精品一区二区三区四区五区乱码| 国产日本99.免费观看| 国产精品亚洲av一区麻豆| 亚洲国产中文字幕在线视频| 国产aⅴ精品一区二区三区波| 久久精品亚洲精品国产色婷小说| 国产精品亚洲一级av第二区| 每晚都被弄得嗷嗷叫到高潮| 国产精品影院久久| 亚洲欧美日韩东京热| 在线观看免费日韩欧美大片| 一个人免费在线观看电影 | 欧美不卡视频在线免费观看 | 18禁裸乳无遮挡免费网站照片| 久久天躁狠狠躁夜夜2o2o| 一本久久中文字幕| 欧美色视频一区免费| 久久久久性生活片| 女人被狂操c到高潮| www.熟女人妻精品国产| 欧美大码av| 成人欧美大片| 精品电影一区二区在线| 91成年电影在线观看| 在线观看免费午夜福利视频| 国产亚洲精品av在线| 国产av在哪里看| 91成年电影在线观看| 国产成人一区二区三区免费视频网站| 国产激情偷乱视频一区二区| 午夜福利高清视频| 中文字幕熟女人妻在线| 老司机福利观看| 特级一级黄色大片| 黄色女人牲交| 欧美色视频一区免费| 国产探花在线观看一区二区| 日韩免费av在线播放| 国产黄a三级三级三级人| 一级作爱视频免费观看| 天堂av国产一区二区熟女人妻 | 一边摸一边抽搐一进一小说| 久久国产精品人妻蜜桃| 日韩免费av在线播放| 久久久久久久午夜电影| 97超级碰碰碰精品色视频在线观看| 9191精品国产免费久久| 巨乳人妻的诱惑在线观看| 两个人视频免费观看高清| 久久久久免费精品人妻一区二区| 精华霜和精华液先用哪个| 国产高清有码在线观看视频 | 2021天堂中文幕一二区在线观| 叶爱在线成人免费视频播放| 成年女人毛片免费观看观看9| 老司机深夜福利视频在线观看| 午夜影院日韩av| 免费看美女性在线毛片视频| 午夜成年电影在线免费观看| 亚洲一区二区三区色噜噜| 天堂动漫精品| 99久久精品国产亚洲精品| 欧美乱码精品一区二区三区| 亚洲片人在线观看| 高清毛片免费观看视频网站| 日日夜夜操网爽| 97碰自拍视频| 亚洲aⅴ乱码一区二区在线播放 | 欧美日韩中文字幕国产精品一区二区三区| 国产97色在线日韩免费| 亚洲精品一卡2卡三卡4卡5卡| 久久天堂一区二区三区四区| 午夜老司机福利片| 亚洲一区中文字幕在线| 国产蜜桃级精品一区二区三区| 又大又爽又粗| 50天的宝宝边吃奶边哭怎么回事| 午夜影院日韩av| 亚洲一码二码三码区别大吗| 日韩大尺度精品在线看网址| 男女午夜视频在线观看| 午夜福利在线观看吧| 国内久久婷婷六月综合欲色啪| 精品少妇一区二区三区视频日本电影| 亚洲精品美女久久久久99蜜臀| 久久久国产成人免费| 一区二区三区国产精品乱码| 五月伊人婷婷丁香| 制服诱惑二区| 一本综合久久免费| 少妇被粗大的猛进出69影院| 制服诱惑二区| 国产亚洲精品一区二区www| 国产精品久久久av美女十八| 午夜福利在线观看吧| 亚洲在线自拍视频| 少妇人妻一区二区三区视频| 国产成人影院久久av| 亚洲在线自拍视频| 久久国产乱子伦精品免费另类| 国产成人精品久久二区二区91| 黄色丝袜av网址大全| 国产一区二区三区视频了| 色综合站精品国产| 国产精品 欧美亚洲| 亚洲电影在线观看av| 19禁男女啪啪无遮挡网站| 国内精品久久久久久久电影| 麻豆成人av在线观看| 91麻豆精品激情在线观看国产| 好男人在线观看高清免费视频| 精品高清国产在线一区| 色在线成人网| 色播亚洲综合网| 日韩欧美国产在线观看| 最好的美女福利视频网| 男人的好看免费观看在线视频 | 两个人视频免费观看高清| 国产片内射在线| 国产亚洲精品久久久久久毛片| АⅤ资源中文在线天堂| 成熟少妇高潮喷水视频| 久久久久久人人人人人| 日日干狠狠操夜夜爽| 欧美成人免费av一区二区三区| 天堂av国产一区二区熟女人妻 | 男人舔女人下体高潮全视频| 丰满人妻一区二区三区视频av | 国内久久婷婷六月综合欲色啪| 国产伦一二天堂av在线观看| ponron亚洲| 国产精品av视频在线免费观看| 叶爱在线成人免费视频播放| 香蕉国产在线看| 精品久久久久久成人av| 久久久久免费精品人妻一区二区| 老汉色∧v一级毛片| 国产精品久久久人人做人人爽| 亚洲七黄色美女视频| 两人在一起打扑克的视频| 男女那种视频在线观看| 又黄又爽又免费观看的视频| 亚洲欧美日韩无卡精品| 国产精华一区二区三区| 搞女人的毛片| 国产黄a三级三级三级人| 欧美丝袜亚洲另类 | 国产免费av片在线观看野外av| 国产欧美日韩一区二区精品| 欧美日韩瑟瑟在线播放| 90打野战视频偷拍视频| 国产激情久久老熟女| 国产一区二区三区视频了| 国产精品精品国产色婷婷| √禁漫天堂资源中文www| 亚洲一区二区三区色噜噜| 久久中文字幕一级| 啦啦啦韩国在线观看视频| 国产野战对白在线观看| 久久久国产成人免费| 老司机深夜福利视频在线观看| 国产私拍福利视频在线观看| 欧美中文日本在线观看视频| 日本黄色视频三级网站网址| 精品免费久久久久久久清纯| 国产激情偷乱视频一区二区| 18禁国产床啪视频网站| 少妇的丰满在线观看| 成人特级黄色片久久久久久久| a级毛片a级免费在线| 特大巨黑吊av在线直播| 欧美又色又爽又黄视频| 1024视频免费在线观看| 超碰成人久久| 日本精品一区二区三区蜜桃| 成人国语在线视频| 日韩免费av在线播放| 色播亚洲综合网| 日日爽夜夜爽网站| 久久中文看片网| 在线看三级毛片| 黑人巨大精品欧美一区二区mp4| 淫秽高清视频在线观看| 欧美久久黑人一区二区| 欧美激情久久久久久爽电影| 九色成人免费人妻av| 99久久综合精品五月天人人| 神马国产精品三级电影在线观看 | 三级国产精品欧美在线观看 | 日韩欧美在线二视频| 国内少妇人妻偷人精品xxx网站 | 高清在线国产一区| 亚洲自拍偷在线| 青草久久国产| av福利片在线观看| 成在线人永久免费视频| 日日夜夜操网爽| 久久精品夜夜夜夜夜久久蜜豆 | 国产99白浆流出| 91av网站免费观看| 成人精品一区二区免费| 亚洲国产欧洲综合997久久,| 欧美成人免费av一区二区三区| 老司机靠b影院| 最近最新中文字幕大全免费视频| 成人一区二区视频在线观看| 香蕉av资源在线| 两人在一起打扑克的视频| 国产精品亚洲一级av第二区| 丝袜美腿诱惑在线| 人妻丰满熟妇av一区二区三区| av国产免费在线观看| 亚洲一卡2卡3卡4卡5卡精品中文| 熟女电影av网| 亚洲五月天丁香| 亚洲av成人一区二区三| 美女高潮喷水抽搐中文字幕| 精品久久久久久久久久久久久| 女人高潮潮喷娇喘18禁视频| av福利片在线| 不卡一级毛片| 欧美成人一区二区免费高清观看 | 国产av在哪里看| 久久性视频一级片| av福利片在线观看| 欧美在线一区亚洲| 国产免费男女视频| 91成年电影在线观看| 久久人妻av系列| 97人妻精品一区二区三区麻豆| 国产熟女午夜一区二区三区| 日本成人三级电影网站| 我要搜黄色片| 国产亚洲精品综合一区在线观看 | 中文资源天堂在线| 久久久国产精品麻豆| 国产aⅴ精品一区二区三区波| 亚洲自拍偷在线| 别揉我奶头~嗯~啊~动态视频| 亚洲精品一卡2卡三卡4卡5卡| 久久中文字幕人妻熟女| 看免费av毛片| 日韩欧美三级三区| 中文字幕最新亚洲高清| 美女午夜性视频免费| 狂野欧美白嫩少妇大欣赏| 亚洲国产高清在线一区二区三| 一夜夜www| 看片在线看免费视频| 宅男免费午夜| 亚洲欧洲精品一区二区精品久久久| 无遮挡黄片免费观看| 别揉我奶头~嗯~啊~动态视频| 三级男女做爰猛烈吃奶摸视频| 91av网站免费观看| 亚洲av成人一区二区三| 亚洲18禁久久av| 午夜激情福利司机影院| 白带黄色成豆腐渣| 香蕉丝袜av| 一本精品99久久精品77| av欧美777| 91国产中文字幕| 国语自产精品视频在线第100页| 听说在线观看完整版免费高清| 观看免费一级毛片| 宅男免费午夜| 欧美大码av| 国产精品香港三级国产av潘金莲| 午夜老司机福利片| 亚洲 国产 在线| 青草久久国产| 在线观看www视频免费| 老鸭窝网址在线观看| 一本大道久久a久久精品| 免费av毛片视频| 狂野欧美白嫩少妇大欣赏| 亚洲国产日韩欧美精品在线观看 | 伦理电影免费视频| 人妻夜夜爽99麻豆av| 久久久久久大精品| 女警被强在线播放| 亚洲精品粉嫩美女一区| 午夜免费激情av| 久久人妻av系列| 夜夜看夜夜爽夜夜摸| 老汉色av国产亚洲站长工具| 毛片女人毛片| 一级毛片精品| 亚洲自拍偷在线| 一个人观看的视频www高清免费观看 | 看黄色毛片网站| 亚洲精品一区av在线观看| 我要搜黄色片| 高清在线国产一区| 少妇粗大呻吟视频| 床上黄色一级片| 国产v大片淫在线免费观看| 老司机深夜福利视频在线观看| 国内揄拍国产精品人妻在线| 国模一区二区三区四区视频 | 亚洲 欧美 日韩 在线 免费| 一本久久中文字幕| 国产亚洲精品久久久久久毛片| 一个人观看的视频www高清免费观看 | 久久久久久久久中文| 久久久水蜜桃国产精品网| 久久国产精品人妻蜜桃| 国产成年人精品一区二区| 日韩欧美国产在线观看| 日韩欧美三级三区| 亚洲欧美日韩高清在线视频| www.999成人在线观看| 国产成人aa在线观看| 悠悠久久av| 欧美成狂野欧美在线观看| 熟女电影av网| 欧美成狂野欧美在线观看| 亚洲熟女毛片儿| 麻豆av在线久日| 国语自产精品视频在线第100页| 亚洲精品中文字幕一二三四区| 亚洲专区字幕在线| 搡老妇女老女人老熟妇| 少妇熟女aⅴ在线视频| 俄罗斯特黄特色一大片| 欧美绝顶高潮抽搐喷水| 国产亚洲欧美在线一区二区| 亚洲国产欧美网| 色综合站精品国产| 久久精品国产综合久久久| a级毛片a级免费在线| 中文字幕高清在线视频| 极品教师在线免费播放| 超碰成人久久| 夜夜爽天天搞| 国产精品亚洲一级av第二区| 两性夫妻黄色片| 欧美精品亚洲一区二区| a在线观看视频网站| 在线观看66精品国产| 夜夜爽天天搞| 熟女电影av网| 熟女电影av网| 好男人在线观看高清免费视频| 亚洲中文字幕一区二区三区有码在线看 | 久久九九热精品免费| 在线观看免费视频日本深夜| 久久九九热精品免费| 国产午夜福利久久久久久| 欧美乱妇无乱码| 国产午夜福利久久久久久| 一本一本综合久久| 天天躁夜夜躁狠狠躁躁| 中文字幕精品亚洲无线码一区| 日韩免费av在线播放| 一进一出抽搐动态| 最好的美女福利视频网| 欧美成人免费av一区二区三区| 美女午夜性视频免费| 美女扒开内裤让男人捅视频| 在线永久观看黄色视频| 淫妇啪啪啪对白视频| 男女床上黄色一级片免费看| 国产1区2区3区精品| 亚洲国产欧美人成| 国产亚洲精品一区二区www| 制服丝袜大香蕉在线| 欧美性长视频在线观看| 亚洲五月天丁香| 精品不卡国产一区二区三区| 后天国语完整版免费观看| 不卡一级毛片| 美女 人体艺术 gogo| 又粗又爽又猛毛片免费看| 国产成人一区二区三区免费视频网站| 村上凉子中文字幕在线| 精品欧美一区二区三区在线| 欧美黑人精品巨大| 中文字幕最新亚洲高清| 久久久水蜜桃国产精品网| 国产激情偷乱视频一区二区| 久久中文字幕一级| 老熟妇乱子伦视频在线观看| a在线观看视频网站| 999久久久精品免费观看国产| 男人的好看免费观看在线视频 | 久久久久久亚洲精品国产蜜桃av| av在线播放免费不卡| 99热这里只有精品一区 | 久久久国产精品麻豆| 看黄色毛片网站| 亚洲av美国av| 丰满人妻熟妇乱又伦精品不卡| 男女做爰动态图高潮gif福利片| 日本三级黄在线观看| 色综合亚洲欧美另类图片| 午夜久久久久精精品| 国产成人系列免费观看| 午夜激情福利司机影院| 亚洲美女黄片视频| 欧美极品一区二区三区四区| 欧美日本亚洲视频在线播放| 黄色片一级片一级黄色片| 久久欧美精品欧美久久欧美| 国内精品久久久久久久电影| 巨乳人妻的诱惑在线观看| 亚洲国产精品久久男人天堂| 天堂动漫精品| netflix在线观看网站| 性色av乱码一区二区三区2| 亚洲国产精品合色在线| 99riav亚洲国产免费| 精品一区二区三区视频在线观看免费| 黄色片一级片一级黄色片| 村上凉子中文字幕在线| svipshipincom国产片| 精品乱码久久久久久99久播| 一级毛片高清免费大全| 听说在线观看完整版免费高清| 免费在线观看视频国产中文字幕亚洲| 亚洲黑人精品在线| 中文亚洲av片在线观看爽| 99在线视频只有这里精品首页| 好男人电影高清在线观看| 国产69精品久久久久777片 | 18禁裸乳无遮挡免费网站照片| 色哟哟哟哟哟哟| 日韩有码中文字幕| 国产野战对白在线观看| 成人18禁在线播放| 19禁男女啪啪无遮挡网站| 91在线观看av| 亚洲国产欧美一区二区综合| 香蕉国产在线看| 一夜夜www| 性色av乱码一区二区三区2| 久久久久精品国产欧美久久久| 18禁观看日本| 两个人视频免费观看高清| 日韩欧美在线乱码| 一区福利在线观看| 欧美中文日本在线观看视频| 波多野结衣高清无吗| 国产av又大| 成人特级黄色片久久久久久久| 欧美zozozo另类| 国产亚洲精品一区二区www| 久久久久久大精品| 窝窝影院91人妻| 日本精品一区二区三区蜜桃| 久久热在线av| 香蕉av资源在线| √禁漫天堂资源中文www| 一a级毛片在线观看| www国产在线视频色| 午夜福利18| 久久久久性生活片| ponron亚洲| 午夜精品久久久久久毛片777| 全区人妻精品视频| 国产成人影院久久av| 小说图片视频综合网站| 久久精品91无色码中文字幕| 超碰成人久久| 亚洲五月天丁香| 国产亚洲av高清不卡| 别揉我奶头~嗯~啊~动态视频| 亚洲真实伦在线观看| 日本三级黄在线观看| 99在线人妻在线中文字幕| 久久精品亚洲精品国产色婷小说| 亚洲精品久久国产高清桃花| 男人舔女人的私密视频| 免费看日本二区| 丰满的人妻完整版| 一级毛片高清免费大全| 黄频高清免费视频| 动漫黄色视频在线观看| 狂野欧美白嫩少妇大欣赏| 国产不卡一卡二| 99热这里只有是精品50| 91av网站免费观看| 亚洲成人免费电影在线观看| 亚洲人与动物交配视频| 亚洲片人在线观看| 亚洲av中文字字幕乱码综合| 最近最新中文字幕大全免费视频| 亚洲 欧美 日韩 在线 免费| 在线看三级毛片| 久久精品国产亚洲av香蕉五月| 日韩大尺度精品在线看网址| а√天堂www在线а√下载| 国产99白浆流出| 国产精品亚洲av一区麻豆| 黄色女人牲交| 亚洲国产欧美网| 久久精品91无色码中文字幕| 大型黄色视频在线免费观看| 99热这里只有精品一区 | 在线永久观看黄色视频| 一级毛片高清免费大全| 亚洲国产精品成人综合色| 每晚都被弄得嗷嗷叫到高潮| 91麻豆av在线| 久久久水蜜桃国产精品网| 少妇熟女aⅴ在线视频| 制服诱惑二区| 亚洲国产精品sss在线观看| 精品久久久久久久久久免费视频| 一个人免费在线观看电影 | 国产精品一及| 久久久久免费精品人妻一区二区| www.www免费av| 日韩欧美国产一区二区入口| 午夜成年电影在线免费观看| 国产亚洲欧美在线一区二区| 欧美性长视频在线观看| 久久九九热精品免费| 亚洲国产精品sss在线观看| 亚洲在线自拍视频| 好男人在线观看高清免费视频| 国产在线观看jvid| 欧美最黄视频在线播放免费| 欧美在线一区亚洲| 又粗又爽又猛毛片免费看| 色尼玛亚洲综合影院| 欧美性猛交黑人性爽| 黄色成人免费大全| 老司机午夜十八禁免费视频| 丰满的人妻完整版| 国产亚洲精品一区二区www| 国产精华一区二区三区| 搡老岳熟女国产| 人妻丰满熟妇av一区二区三区| 最近最新中文字幕大全免费视频| 老熟妇仑乱视频hdxx| 法律面前人人平等表现在哪些方面| 久久精品综合一区二区三区| 又紧又爽又黄一区二区| 在线看三级毛片| 久久久久久久久久黄片| 少妇熟女aⅴ在线视频| 国产熟女午夜一区二区三区| 亚洲片人在线观看| 在线国产一区二区在线| 丝袜人妻中文字幕| 精品人妻1区二区| 久久午夜综合久久蜜桃| 男女之事视频高清在线观看| 少妇粗大呻吟视频| 成人18禁在线播放| 国产精品av久久久久免费| 日韩中文字幕欧美一区二区| 欧美一级毛片孕妇| 国产精品九九99| 亚洲专区国产一区二区| 亚洲成a人片在线一区二区| 中文字幕人妻丝袜一区二区| 亚洲一区中文字幕在线| 久久久精品国产亚洲av高清涩受| 婷婷六月久久综合丁香| 香蕉久久夜色| 制服诱惑二区| 一区福利在线观看| 亚洲国产日韩欧美精品在线观看 | 精品久久久久久久末码| 日韩免费av在线播放| 国产成人啪精品午夜网站| 国产精品自产拍在线观看55亚洲| 日本黄色视频三级网站网址| 成年免费大片在线观看| АⅤ资源中文在线天堂| 在线观看66精品国产| 国产成人欧美在线观看| 亚洲熟妇熟女久久| 国产精品1区2区在线观看.| 国内精品久久久久久久电影| 男女床上黄色一级片免费看| 老熟妇乱子伦视频在线观看| 欧美色欧美亚洲另类二区| 十八禁网站免费在线| 手机成人av网站| 欧美色欧美亚洲另类二区| 十八禁网站免费在线| av福利片在线观看| 国产成人欧美在线观看| 亚洲国产欧美人成| 亚洲国产高清在线一区二区三| 一卡2卡三卡四卡精品乱码亚洲| 中文字幕熟女人妻在线| 日韩中文字幕欧美一区二区| 校园春色视频在线观看| 精品久久久久久久毛片微露脸| 亚洲一区二区三区色噜噜|