Coronal Properties of Low-mass Population III Stars and the Radiative Feedback in the Early Universe

We carried out one-dimensional magnetohydrodynamic (MHD) simulations of the heating of coronal loops in metal-free stars. To find the systematic properties, we perform the simulations with stellar mass $M=0.8,0.7,0.6,0.5M_{\odot}$ and a wide range of the magnetic field strength.

The loop model and simulations are basically similar to those in Washinoue & Suzuki (2019); we consider a single magnetic flux tube that is anchored at the photosphere and take the coordinate $s$ along the loop. We inject velocity perturbations $\delta v$ in three directions from the footpoints to excite MHD waves.

It is assumed that the loop expands towards a higher altitude with an expansion factor, $f(s)$:

where $f_{\rm max}$ is the value of $f(s)$ at the loop top. $R_{\star}$ is the stellar radius and $h$ is the loop height. $a$ is a constant for each flux tube and we assume a functional form of $a=-0.54h/10^{4}+0.6\times\mathrm{log}_{2}(f_{\rm max}/100)+5.3$, where we set $b=0.02$. The magnetic field along the loop is described as $B_{s}=B_{\rm ph}/f(s)$, where $B_{\rm ph}$ is the photospheric magnetic field strength.

We adopt the loop length $l=8\times 10^{4}$km, which gives $h\simeq 2.5\times 10^{4}$ km, for $M=0.8M_{\odot}$. We note that these are roughly comparable to moderate-sized solar coronal loops. The loop lengths for different stellar masses are set to be proportional to the scale height $\propto a_{\rm ph}^{2}/g$ where $a_{\rm ph}$ and $g$ are the sound speed and the gravitational acceleration at the photosphere, respectively (See Table 1).

We solved the ideal MHD equations including thermal conduction and radiative cooling as follows;

and

where $G$ is the gravitational constant and $R$ is the distance from the center of a star. $e$ is the internal energy, which has the relation $e=\frac{P}{(\gamma-1)/\rho}$ where $\gamma=5/3$. $F_{c}$ is the thermal conductive flux which has the form $F_{c}=\kappa T^{5/2}\frac{\partial T}{\partial s}$, where $\kappa$ is the Spitzer conductivity. $\mu$ is the mean molecular weight that we adopt $\mu=0.6$ for fully ionized plasma, $m_{u}$ is the atomic mass unit and $k_{B}$ is the Boltzmann constant. Note that when we solve Equations (2)-(6), it is taken into account the following curvature effects using Equation (1);

where $\bm{X}$ and $X_{s}$ are any vector and its $s$ component.

$q_{R}$ is radiative cooling rate per volume. Since our simulations treat the wide range of atmospheric layers from the photosphere to the corona, we adopt different prescriptions for the chromosphere and the upper regions. In the chromosphere, we utilized the formulation by Suzuki (2018) which is extended from the empirical model for the solar atmosphere (Anderson & Athay, 1989);

Equation (9) is derived based on the observation of the solar chromosphere; 80 % of the chromospheric emission is contributed from heavy elements.

In the optically thin region,

$\Lambda$ is the cooling function, for which we adopt from Sutherland & Dopita (1993) for $Z=0$. $n$ and $n_{e}$ are ion and electron number densities.

In addition, we define the cutoff temperature of radiative cooling, $T_{\rm off}$, such that $q_{R}=0$ when $T<T_{\rm off}$. It controls a minimum temperature of the low atmosphere, $T_{\rm min}$. Although $T_{\rm min}$ is often lower than $T_{\rm off}$ because the adiabatic cooling works even when $q_{\rm R}$ is turned off, $T_{\rm min}$ is generally comparable to $T_{\rm off}$. A classical model of the solar chromosphere gives $T_{\rm min}\approx 4000$ K $\approx 0.7T_{\rm eff}$ (Vernazza et al., 1981). Because $T_{\rm min}$ for low-metallicity stars has not been specified in observations, we simply refer to this classical solar value to set $T_{\rm off}=0.7T_{\rm eff}$.

We adopt the 2nd order Godunov method for compressible waves and the method of characteristics for incompressible waves (Stone & Norman, 1992; Suzuki & Inutsuka, 2005).

In our model, the initial cool gas with $T=T_{\rm eff}$ is heated via the dissipation of Alfvén waves. In our simulations, the main channel of the wave dissipation is the nonlinear mode conversion from incompressible Alfvénic waves to compressible waves (Suzuki & Inutsuka, 2005, 2006)). In addition, fast-mode MHD shocks, which is formed by the direct steepening of transverse (Alfvénic) waves (Hollweg, 1982; Suzuki, 2004), is another channel of the wave dissipation.

Figure 1 shows the time evolution of the loop profiles for the cases of $M=0.8M_{\odot}$ and $B_{\rm c}=0.18,2.9$ and $11.6$ G. As can be seen from the figure, most of our simulations indicate that the simulated coronal loops undergo the dynamical evolution with intermittent rises and drops of the temperature (top panels). Such thermally non-equilibrium behaviors have been actually seen in observations (Auchère et al., 2014; Froment et al., 2015) as well as simulations (Müller et al., 2003, 2004; Antolin et al., 2010; Froment et al., 2018). However, the dynamical behavior in our simulations may be partially because of our limited treatment of the wave heating in 1D flux tubes (see Section 4.2).

The left panels in Figure 1 correspond to the case with weak $B_{\rm c}=0.18$ G. In this case, the corona in excess of $T=10^{6}$K is not formed. It is only heated up to $T\sim(3-5)\times 10^{5}$K and fluctuates below $10^{6}$K . Smaller $B_{\rm c}$ leads to the larger nonlinearity of Alfvénic waves $\propto\delta B/B_{s}$, which enhances the wave dissipation. Therefore, a larger fraction of the Alfvénic waves dissipate at the lower altitude where the density is higher and most input energy is lost through radiative cooling. In contrast, only a smaller fraction of the initial wave energy remains at higher altitudes, and therefore, the upper atmosphere is not heated up to the coronal temperature.

It is also useful to evaluate characteristic timescales. We define the heating timescale $\tau_{\rm heat}$ as the travel time required for the Alfvén wave to damp the Poynting flux $F_{\rm A}$ by a factor of e:

where we evaluate $v_{\rm A}$ and $F_{\rm A}$ at the loop top, which yields $\tau_{\rm heat}=250$ s for this case. The radiative cooling timescale is given by

We note that the dependence of $\tau_{\rm cool}$ on $T$ is different for different temperature ranges; for $T>10^{6}$ K, $\Lambda\propto T^{1/2}$ so that $\tau_{\rm cool}\propto T^{1/2}\rho^{-1}$. On the other hand, for $T<10^{6}$ K, $\Lambda\propto T^{-1}$ in the case of zero metallicity (Sutherland & Dopita, 1993; Washinoue & Suzuki, 2019) and thus $\tau_{\rm cool}\propto T^{2}\rho^{-1}$. For both temperature ranges, $\tau_{\rm cool}$ increases with temperature so that the gas should be thermally unstable (Balbus, 1995). However, for $T>10^{6}$ K, the gas is not so unstable because the dependence on $T$ is relatively weak and the thermal conduction $\propto T^{5/2}$ works to stabilize it. In contrast, for $T<10^{6}$ K, the gas is quite unstable owing to the inefficient thermal conduction and strong dependence of $\tau_{\rm cool}$ on $T$.

If we estimate $\tau_{\rm cool}$ for the time-averaged values of $T=2.5\times 10^{5}$K and $\rho=4\times 10^{-14}$ g cm^-3 at the loop top, we obtain $\tau_{\rm cool}=180$ s. The comparison between $\tau_{\rm heat}$ and $\tau_{\rm cool}$ indicates that the heating is dominated by the radiative cooling in the weak $B_{\rm c}$ case. Therefore, high-temperature corona is not formed in this case.

In the case of intermediate $B_{\rm c}=2.9$ G (the middle panels), the corona is rapidly heated up to $T=10^{6}$K and reaches $T\sim 3\times 10^{6}$K at $t=29.1$hr. In this case, the wave nonlinearity, $\delta B/B_{s}$, is moderately smaller and a larger fraction of the input Alfvénic waves reach the corona to heat the gas. The heating timescale is $\tau_{\rm heat}=310$ s, which is consistent with the slower dissipation compared with the weak $B_{\rm c}$ case. The cooling timescale is $\tau_{\rm cool}=480$ s, where we use the time-averaged loop top values of $T=2.2\times 10^{6}$K and $\rho=2.0\times 10^{-13}$ g cm^-3. $\tau_{\rm cool}\propto T^{1/2}$ is long enough owing to the high temperature, and consequently the cooling is dominated by the heating, $\tau_{\rm cool}>\tau_{\rm heat}$. Therefore, the loop profile is kept almost steady with time once the temperature is in the thermally stable region of $T>10^{6}$K.

The right panels show the case with strong $B_{\rm c}=11.6$ G. The corona is heated up to $T\sim 10^{6}$K, however, it does not settle into the stable state. In this case, $\tau_{\rm heat}=440$ s and $\tau_{\rm cool}=310$ s using $T=7\times 10^{5}$K and $\rho=1.0\times 10^{-13}$ g cm^-3 for the time-averaged values. However, it is found that $\tau_{\rm cool}$ fluctuates widely over time; the gas is efficiently heated at early time and it leads more chromospheric evaporation to increase the coronal density. Then, the cooling time drops to $\tau_{\rm cool}\sim 1$ s and the denser gas suffers more effective radiative cooling so that the temperature drops to the thermally unstable range of $T\lesssim 5\times 10^{5}$K. As the loop structure collapses and the density also decreases, $\tau_{\rm cool}$ increases to $\sim 10^{3}$ s and the corona is reheated to $T\sim 10^{6}$K. The loop structure evolves dynamically with time in the temperature range $10^{5}$K $<T<10^{6}$K in a quasi-periodic manner.

The more dynamic nature of stronger $B_{\rm c}$ cases is also interpreted in terms of the intermittency of the wave heating. The heating in the loops are done by shock waves that are nonlinearly excited from Alfvénic waves. The wavelength of the Alfvénic waves is longer for larger $B_{\rm c}$ owing to the faster Alfvén velocity, which also gives a longer spatial interval between neighboring shocks. The shock heating does not occur in a continuous manner but at localized regions in a stochastic manner in stronger $B_{\rm c}$ cases. Therefore the stable hot corona is hardly kept but the temperature goes up and down in a cyclic manner.

To compare the average properties of loops with different $B_{\rm c}$ and $M$, we show the time-averaged distributions of temperature and density in Figure 2. The profiles for different $B_{\rm c}$ are displayed by different line colors in each panel; the blue, black and red lines correspond to weak ($<1$ G), intermediate ($\geq 1$ G and $\leq 10$ G) and large ($>10$ G) $B_{\rm c}$, respectively. The parameter sets $(B_{\rm ph}/B_{\rm ph*},B_{\rm c})$ for the different cases are shown in Table 2.

The dependence on $B_{\rm c}$ show a similar tendency for different masses: the hottest and densest corona is obtained in the cases with the intermediate $B_{\rm c}$. For the smaller $B_{\rm c}$, the fast dissipation in the lower atmosphere makes it difficult to form the hot corona. For the larger $B_{\rm c}$, the higher density enhances the radiative cooling and the coronal temperature of $T>10^{6}$K is hard to be maintained for a long time. Thus, the loop dynamically evolves with time and the time-averaged temperature and density are lower than those of the intermediate cases. However, we note that the behavior of the cases with strong coronal field, $B_{\rm c}>10$ G may be modified if we consider additional heating processes, such as the turbulent cascade of Alfvénic waves, which is not taken into account in our simulations. We will discuss this issue in Section 4.2.

When we focus on the same line colors in each panel, it can be seen that the density is higher for lower-mass stars. The gas pressure at the photosphere is larger for lower-mass stars. In addition, the loop height, which is proportional to the pressure scale height (Section 2), is shorter for lower-mass stars. Therefore, the coronal density of lower-mass stars is higher even though the difference between the photospheric density and the coronal density is higher.

We estimated the radiative fluxes in extreme ultraviolet (EUV) and X-ray wavelengths to understand the correlation between the coronal magnetic field and the high-energy radiation. We here assume that the star is covered by the same simulated loops and calculate the luminosity $L$ as follows;

We label $L_{\rm EUV}$ and $L_{\rm X}$ for $L$ in the temperature ranges of $1.5\times 10^{5}\mathrm{K}\leq T\leq 1.1\times 10^{6}$K and $T>1.1\times 10^{6}$K, respectively. We note that $T=1.5\times 10^{5}$K and $=1.1\times 10^{6}$K respectively correspond to 13.6 eV and 0.1 keV.

The left panel of Figure 3 shows the time-averaged EUV luminosities $L_{\rm EUV}$ of all the simulated data against $B_{\rm c}$. We note that multiple cases with different sets of the input parameters of $B_{\rm ph}$ and $f_{\rm max}$ (Section 2.2.2) give the same $B_{\rm c}$ for each $M$. $L_{\rm EUV}$ is peaked around $B\approx$ a few G, which is consistent with Figure 2 and the discussion in Section 3.2; the cases with intermediate $B_{\rm c}=2-6$ G give the hot and dense coronae. Toward weaker $B_{\rm c}<1$ G and stronger $B_{\rm c}>10$ G, $L_{\rm EUV}$ gradually decreases. In the top panels of Figure 2, the horizontal dashed lines indicate $T=1.5\times 10^{5}$K and $T=1.1\times 10^{6}$K, and we define that the EUV radiations are from the plasma between these temperatures. The filled boxes and the open circles in the bottom panels point the location of these temperatures. Figure 2 illustrates that, although the hot corona is not formed in small and large $B_{\rm c}$ cases, EUV radiations are emitted from broad regions of the loops in these cases.

Some cases occasionally emit quite strong radiation although the average temperature around the loop top is low $\lesssim 5\times 10^{5}$ K. We show the upper limits of these cases in the circles with arrows in Figure 3. Some simulated cases with weak $B_{\rm c}$ give quite small $L_{\rm EUV}$, which is below the displayed range of $L_{\rm EUV}\geq 10^{26}$erg s^-1. These cases, which are indicated by the single arrows in the bottom left side of the panel, imply that there is a condition for the coronal magnetic field strength to emit sufficiently large $L_{\rm EUV}$.

The right panel of Figure 3 shows the EUV luminosity normalised by the bolometric luminosity, $L_{\rm EUV}/L_{\rm bol}$, with $B_{\rm c}$. We derived $L_{\rm bol}$ from $L_{\rm bol}=4\pi R_{\star}^{2}\sigma T_{\rm eff}^{4}$, where $\sigma$ is the Stefan-Boltzmann constant and we use the values of $R_{\star}$ and $T_{\rm eff}$ from Table 1. In this panel, the dependence on stellar mass is clearly visible; lower-mass stars give higher $L_{\rm EUV}/L_{\rm bol}$ because the gas density in the loops are higher. The peak value, $L_{\rm EUV}/L_{\rm bol}>10^{-4}$ for $M=0.5M_{\odot}$, while $L_{\rm EUV}/L_{\rm bol}\sim 8\times 10^{-6}$ for $M=0.8M_{\odot}$.

Figure 4 shows the time-averaged X-ray luminosities $L_{\rm X}$ and $L_{\rm X}/L_{\rm bol}$ with $B_{\rm c}$. The trend is almost similar to that of $L_{\rm EUV}$, while it has a strong dependence on $B_{\rm c}$ compared with $L_{\rm EUV}$. $L_{\rm X}$ decreases toward both weaker and stronger sides of $B_{\rm c}$. In the weak $B_{\rm c}$, the X-ray corona with $T\geq 1.1\times 10^{6}$K is hardly heated up (Figures 1 & 2). In the strong $B_{\rm c}$, on the other hand, although the X-ray corona is formed transiently, the temperature drops in a quasi-periodic manner. Therefore, the time-averaged $L_{\rm X}$ is lower than those of intermediate $B_{\rm c}$ cases.

$L_{\rm X}/L_{\rm bol}$ is higher for lower mass as with $L_{\rm EUV}/L_{\rm bol}$. At the peak, $L_{\rm X}/L_{\rm bol}\sim 5\times 10^{-5}$ for $M=0.5M_{\odot}$, while $L_{\rm X}/L_{\rm bol}\sim 10^{-5}$ for $M=0.8M_{\odot}$. Some cases with $M\geq 0.7M_{\odot}$ emit strong X-rays, $L_{\rm X}/L_{\rm bol}\gtrsim L_{\rm EUV}/L{\rm bol}$.

We additionally show the dependence of $L_{\rm X}/L_{\rm bol}$ on the other magnetic parameters of $B_{\rm ph}$ and $f_{\rm max}$, in Figure 5. It can be seen that both are related to $L_{\rm X}/L_{\rm bol}$; the value of $L_{\rm X}/L_{\rm bol}$ tends to decrease at smaller and larger sides of $B_{\rm ph}$ and $f_{\rm max}$. However, the scatters on the same $B_{\rm ph}$ or $f_{\rm max}$ are larger than those on the same $B_{\rm c}$ (Figure 4). Therefore, the dependence on $B_{\rm ph}$ and $f_{\rm max}$ is relatively weak, which indicates that $B_{\rm c}$ is a critical parameter that determines the radiative property of metal-free coronae.

Finally, we fitted the dependence of $L_{\rm EUV}/L_{\rm bol}$ and $L_{\rm X}/L_{\rm bol}$ on $B_{\rm c}$;

where $M$ is the stellar mass. The dashed lines in the right panels of Figures 3 and 4 indicate the fitted curves for the different masses. $L$ is scattered on the same $B_{\rm c}$ because of the difference of the field strength, $B_{\rm ph}$, at the photosphere, however, the fitting formulae using only $B_{\rm c}$ roughly explain the tendency. Again, we mention that the coronal properties at $B_{\rm c}>10$ G may be modified when other heating processes are taken into account. Therefore, the trends of $L$ in the strong $B_{\rm c}$ are not definite, which will be discussed in the next section.

Low-mass main sequence stars with surface convection show a large variety of magnetic activity involving coronae and flares. Young solar-type stars generally exhibit brighter and hotter X-ray emissions (Güdel et al., 1997; Ribas et al., 2005). For example, EK Dra (HD129333) is a young solar analogue with age $\lesssim 0.1$ Gyr and gives the X-ray luminocity, $L_{\rm X}\simeq 10^{30}$erg s^-1, or log$\left(L_{\rm X}/L_{\rm bol}\right)\simeq-3.6$, which are about $10^{3}$ times more luminous than those of the Sun (Güdel et al., 1995; Telleschi et al., 2005).

The X-ray luminosity is correlated with stellar rotation period (Pallavicini et al., 1981; Pizzolato et al., 2003; Wright et al., 2011), which generally decreases with the stellar age (e.g., Skumanich, 1972). The solar-type stars follow the relation, $L_{\rm X}\approx(3\pm 1)\times 10^{28}t^{-1.5\pm 0.3}$erg/s, where $t$ is the stellar age in Gyr with a saturated luminosity $L_{\rm X}\leq 4\times 10^{30}$ erg/s (Güdel, 2007). The strong X-rays from these (near) saturated stars probably originate from explosive flares that involve large loops (Shibata & Yokoyama, 1999), in addition to the steady quiet corona as modeled in the present paper. Indeed, a large flaring loop whose height is approximately the stellar diameter has been directly observed in the Algol system (Peterson et al., 2010).

In our calculations, we fix the loop length for each stellar mass (Table 1) by extrapolating that of typical solar coronal loops. Also, our model cannot treat magnetic reconnection that is believed to trigger flares (Tsuneta, 1996; Shibata, 1999; Takasao et al., 2012). In more realistic treatments, $L_{\rm X}$ from low-mass Pop III stars is expected to show a larger variety. In addition, observational analysis with the large number of samples reports that the stellar spin evolution depends on metallicity and lower-metallicity stars rotate faster on average than higher-metallicity stars (Amard et al., 2020). It supports that young low-mass Pop III stars have a potential to emit much brighter X-rays for a long duration than the simulated values shown in the previous section.

In our simulations Alfvénic waves dissipate only by the excitation of compressible-mode waves, which is one of the limitations of our one dimensional treatment. In reality, however, other processes of the wave dissipation are also expected to play a role in the heating of the plasma, such as Alfvénic turbulence (Matthaeus et al., 1999; Cranmer et al., 2007; Shoda et al., 2019), phase mixing (Heyvaerts & Priest, 1983; Ofman & Aschwanden, 2002; Matsumoto & Suzuki, 2012), and resonant absorption (Ionson, 1978; Antolin et al., 2015). In order to incorporate these processes we have to handle 2D/3D simulations or adopt additional models, which we plan to study in our future work.

We here discuss Alfvén wave turbulence, which is one of the reliable mechanisms that heat up the corona. Upgoing Alfvénic waves interact with reflected backward propagating waves, which triggers cascading Alfvén wave turbulence. This is an efficient mechanism of the dissipation of Alfvénic waves. Recent studies report that both the shock and turbulent heating processes are essential in driving the solar wind (Cranmer et al., 2007; Matsumoto & Suzuki, 2014; Shoda et al., 2018; Matsumoto, 2020) and coronal loops (Matsumoto, 2016); the shock heating is dominant in the lower atmosphere, while the turbulent heating is dominant above the transition region.

In our simulation with $B_{\rm c}>10$ G, a larger fraction of the input energy of the Alfvénic waves travels into the upper atmosphere owing to the smaller nonlinearity $\delta B/B_{s}$. Therefore, the hot and dense corona is expected to be maintained in the cases with larger $B_{\rm c}>10$ G by the contribution from the turbulent cascade of the Alfvénic waves; we might underestimate $L_{\rm EUV}$ and $L_{\rm X}$ at $B_{\rm c}>10$ G. In contrast, the trend in the regime of $B_{\rm c}<1$ G would not be affected, even if the turbulent heating was active, because most of the input Alfvénic waves are damped before reaching the upper atmosphere. To summarise, we speculate that $L$ does not decrease at $B_{\rm c}<10$G and it may be flat or possibly increases as presented by the pink shaded regions in Figure 3 and 4. It is required to incorporate the turbulent heating process to verify the tendency for our future work.

The chromosphere is an interface that transfers the energy from the photosphere to the corona. The temperature in the chromosphere determines the distribution of the density, which controls the propagation and dissipation of waves. The thermal properties of the chromosphere affects the heating of the upper layers. We adopt the simplified treatment for the radiative cooling in the chromosphere of metal-free stars (Equation (9)). The cooling rate is a simple extension from the solar value (Anderson & Athay, 1989), which is modeled by the empirical temperature profile of the VAL model C under the assumption of the static chromosphere with spatially uniform radiative emissions (Vernazza et al., 1981). However, the solar chromosphere has been known to be quite dynamic (Katsukawa et al., 2007; De Pontieu et al., 2014; Skogsrud et al., 2015) and the local thermodynamic equilibrium (LTE) is not satisfied in general. In such circumstances, it is required to solve radiative transfer and detailed rate equations in time-dependently to determine the radiative cooling rate.

The dynamic chromosphere has been studied by theoretical models and numerical simulations (Carlsson & Stein, 2002; Rammacher & Ulmschneider, 2003) and it is claimed that dynamic models show quite different radiative properties from a static model (Carlsson & Stein, 2002). Carlsson & Leenaarts (2012) introduced an approximated recipe to estimate the chromospheric emission from their radiative transfer calculation. A future direction for our study is to extend their recipe to chromosphers with different metallicities.

The cutoff temperature of the radiative cooling, $T_{\rm off}=0.7\times T_{\rm eff}$, is also an important parameter because it controls the minimum temperature of the atmosphere. Smaller $T_{\rm min}$ leads to the more rapid decrease of the density, which enhances the reflection of Alfvénic waves. The modified treatment may give better understandings of $T_{\rm min}$.

We additionally mention the effect of the magnetic diffusion. The chromosphre is partially ionized and the collision between neutrals and charged particles causes the ambipolar diffusion. It facilitates damping the Alfvénic waves and heating the gas particularly in the upper chromosphere (Singh et al., 2011; Khomenko & Collados, 2012; Martínez-Sykora et al., 2017). The ambipolar diffusivity is described as $\eta_{\rm ambi}=(|B|\rho_{n}/\rho)^{2}/\rho_{i}\nu_{in}$, where $\rho_{n},\rho_{i},\rho$ and $\nu_{in}$ are neutral density, ion density, total density and ion-neutral collision frequency (Martínez-Sykora et al., 2015). Although recent numerical simulations of the solar atmosphere indicate that ambipolar diffusion has a small effect on the chromospheric heating compared with the shock heating which is resulted from the nonlinear mode coupling (Arber et al., 2016; Brady & Arber, 2016), it may be important under the metal-deficient environment; the metal-free condition lowers $\rho_{i}$ because of the lack of the ionization sources and raises $\eta_{\rm ambi}$. Therefore, the effect of amibipolar diffusion would be large and needs to be added in the induction equation for proper modeling of metal-free chromospheres.

The role of low-mass Pop III stars has not been investigated in detail in the context of the cosmic reionization to date, partly because the existence of low-mass Pop III stars is unclear and the high-energy radiation from them is quite uncertain. In this subsection, we discuss the effect of the EUV and X-ray radiations from low-mass Pop III stars on the cosmic reionization in comparison to other ionizing sources.

One of the unique characteristics of low-mass Pop III stars is that they also emit soft X-rays, in addition to the EUV radiation. Our simulations show that their X-ray luminosity is roughly comparable to the EUV luminosity (Section 3.3), which is in distinctive contrast to massive metal-free stars. The X-ray photons also play a role in partially ionizing the neutral gas. Previous studies report that the X-ray sources can contribute a few per cent of the ionization (Oh, 2001; Mirabel et al., 2011). Furthermore they can heat up intergalactic medium (IGM hereafter) in a more spatially uniform manner due to their long mean path (Madau et al., 2004; Ricotti & Ostriker, 2004; McQuinn, 2012). X-ray binaries (XRBs) and mini quasars have been considered to be likely sources of the X-ray reionization. In addition to these candidates, we consider whether the X-rays from low-mass Pop III stars can contribute to the process of the reionization.

We estimate the X-ray luminosity per unit star-formation rate (SFR) to directly compare to previous works that discuss other X-ray sources. We define the total X-ray emission from low-mass Pop III stars per galaxy as follows;

where $N$ is the number of low-mass Pop III stars in a galaxy and $L_{\rm X}$ is the typical X-ray luminosity of a single star. We adopt $L_{\rm X}=1.0\times 10^{28}$ erg/s as the typical luminosity for $B_{\rm c}\sim(2-6)$G that gives the peak in Figure 4.

For estimating $N$ in Equation (18), we have to put extra assumptions about the formation of the Pop III stars, which are quite uncertain. Assuming the mass fraction $f_{\rm cor}$ of the stars with corona whose luminosity is $\approx$ $L_{\rm X}$, we obtain

where $\tau=0.5$ Gyr is the typical duration of the cosmic reionization (Fan et al., 2006); we adopt a possible upper limit for the choice of $\tau$ as the duration of Pop III star formation throughout the reionization. In the later time of the duration of $\tau=0.5$ Gyr, second-generation low-mass Pop II stars contribute to the reionization in addition to the Pop III component. Our previous works show that the X-ray luminosity from coronae with $Z\leq 10^{-3}Z_{\odot}$ is comparable to that with $Z=0$ (Suzuki, 2018; Washinoue & Suzuki, 2019). Therefore, $\tau=0.5$ Gyr is a reasonable estimate. We assume that all the Pop III stars with $0.3M_{\odot}<M<0.9M_{\odot}$ emit $L_{\rm X}=1.0\times 10^{28}$erg s^-1. We adopt the normalization, $f_{\rm cor}=0.3$, from the Kroupa IMF (Kroupa, 2001; Weidner & Kroupa, 2004). $0.6M_{\odot}$ in Equation (19) is chosen as the median stellar mass in this range. In addition, SFR is here assumed to be $1M_{\odot}$yr^-1 per galaxy during $\tau$. We note that SFR and the mass fraction we used are those in the local universe for simplicity because there are little constraints on the formation of low-mass Pop III stars. Therefore, $f_{\rm cor}$ should be smaller if the Pop III IMF has a top-heavy shape (e.g., Bromm & Larson, 2004), and accordingly the value of SFR would be suppressed.

We then obtain the ratio of X-ray emission to the SFR,

Referred to previous studies, the expected ratio for the HMXBs is $\sim 10^{40}f_{\rm X}\rm{erg/s}(M_{\odot}\rm{yr}^{-1})$ (Furlanetto, 2006; Fragos et al., 2013), where $f_{\rm X}$ is the X-ray efficiency and its standard value is $f_{\rm X}=1$. This relation is based on the observation of nearby galaxies and we here assume that it can be extrapolated to higher redshifts. However, $f_{\rm X}$ in higher redshifts is still under investigation and it is sometimes considered to be higher values (Mesinger et al., 2013; Xu et al., 2014). In addition, the average contribution from mini-quasars is about 10% of that from XRBs, $\sim 10^{39}f_{\rm X}$ erg s^-1 (Fialkov et al., 2014). If we adopt the normalized value, $L_{\rm X}=10^{28}$erg s^-1, based on our results, the contribution from the low-mass Pop III stars is negligible in comparison to those from XRBs and mini-quasars.

Next, we evaluate the effect on IGM heating by the soft X-rays. We first consider the energy conservation in the expanding universe;

where $T_{\rm K}$ is the kinetic temperature of the gas and $n_{\rm H}$ is the number density of neutral hydrogen. $H(z)$ is the Hubble parameter and $\epsilon_{i}$ is the energy injected to the IGM through process $i$. We here focus on X-ray heating by low-mass stars and describe the energy as follows;

where we adopt the $L_{\rm X}$-SFR relation from Equation (20). $f_{\rm heat}$ is the fraction of the X-ray energy that is used to heat the IGM. The value of $f_{\rm heat}$ depends on the ionized fraction of the gas and the X-ray photon energy and we adopt $f_{\rm heat}=0.25$, which is slightly higher than that of other X-ray sources owing to the softer spectra of low-mass stars (Ricotti et al., 2002; Furlanetto & Stoever, 2010). $\dot{\rho}_{\rm SFR}$ is the star formation rate density and we use $\dot{\rho}_{\rm SFR}=10^{-2}M_{\odot}\rm{yr}^{-1}\rm{Mpc}^{-3}$ at $z=9$ from the parameterized model of the evolving $\dot{\rho}_{\rm SFR}$ in Robertson et al. (2015). Then, the contribution to the gas temperature due to X-ray heating is

which is below the cosmic microwave background temperature $T_{\rm CMB}=2.73(1+z)$ K.

Therefore, if we take our simulations at face value, the low-mass Pop III stars do not make a critical contribution to the IGM heating. However, we mention a possibility that they produce stronger radiation which is not considered in our simulations. As we noted in Section 4.1, young solar-type stars are magnetically active and emit strong X-rays. The low-mass Pop III stars that could contribute to the reionization should be inevitably young ones so that it is important to take their property into account. While the saturated X-ray luminosity, $L_{\rm X}\approx 4\times 10^{30}$ erg/s, of solar-type stars is 400 times larger than our results $L_{\rm X}\approx 10^{28}$ erg/s, such strong emission would be also possible for Pop III stars with long loops and strong magnetic environment. If $L_{\rm X}$ of young Pop III stars is comparable to the saturated luminosity of solar-type stars, Equation (20) then yields $L_{\rm X,gal}/\rm{SFR}\approx 1.0\times 10^{39}$ erg/s ($M_{\odot}\rm{yr}^{-1}$).

In addition, we showed that the X-ray radiation from moderate-sized loops are stronger for smaller metallicity (Washinoue & Suzuki, 2019). Therefore, the saturated luminosity of low-mass Pop III stars is possible to exceed $L_{\rm X}\geq 4\times 10^{30}$ erg/s. Accordingly, $L_{\rm X,gal}/\rm{SFR}$ would be even larger than the above value and it could be comparable to that of other X-ray objects. It also provides a larger effect on the IGM heating; the gas temperature from the X-ray heating could be more than $\Delta T_{\rm K}\approx 50$ K, which exceeds the CMB temperature at high redshifts. While the heating by other sources such as HMXBs is still more effective (e.g., $\sim 10^{3}f_{\rm X}\frac{1+z}{10}$ K in Furlanetto (2006)), low-mass stars may have a role to promote the preheating due to their softer energy spectra and bring forward the moment of the heating transition at which $T_{\rm gas}=T_{\rm CMB}$ to appear in the 21cm signals. Given the strong X-ray activities which we have not considered in this study, low-mass Pop III stars could be a reasonable contributor to the cosmic reionization. To examine the coronal properties of young Pop III stars, we plan to perform the simulations with strong magnetic fields and long loops for our future work.

The low-mass Pop III stars with the simulated masses emit the EUV ($\geq 13.6$eV) radiation with $L_{\rm EUV}\sim 10^{28}$ erg/s in a wide range of $B_{\rm c}$ (Figure 3). On the other hand, the massive metal-free stars have the effective temperature $T_{\rm eff}\sim 10^{5}$K which has a weak dependence on stellar mass (Schaerer, 2002; Bromm et al., 2001) and emit a large amount of ionizing radiation. For instance, the stars with $M_{\star}=120M_{\odot}$ emit $1.4\times 10^{50}$ hydrogen ionizing photons per second during their lifetime $\tau_{\rm life}=2.5$Myr (Schaerer, 2002). Assuming that all the photons have an energy of $h\nu\sim 13.6$ eV, the emission from a massive Pop III star roughly corresponds to $L_{\rm EUV}\approx 10^{39}$ erg/s which is much larger than the typical $L_{\rm EUV}$ from a low-mass star; we here note that, since the EUV radiation from low-mass stars is weighted on the higher-energy side than that from massive stars, the efficiency of the ionization is higher owing to multiple times of ionization by an ionizing photon.

We next estimate how much gas the emitted EUV photons from low-mass Pop III stars can ionize. We determine the Strömgren radius with the expected number of EUV photons $N_{\rm EUV}$ as follows;

where $\alpha=2.6\times 10^{-13}$ cm³s^-1 is the hydrogen recombination coefficient at $T=10^{4}$ K. Using our typical results of $L_{\rm EUV}\sim 10^{28}$ erg/s and assuming that all the photons have the same energy of $h\nu=13.6$ eV,

Incidentally, a similar calculation with a massive star ($N_{\rm EUV}=1.4\times 10^{59}$ /s) gives $R_{\rm s}\approx 5.0\times 10^{20}$ cm. Because the radiation from a low-mass Pop III star is much smaller than that from a higher-mass counterpart, the ionized region produced by a single stars is also small.

Finally, we estimate the total contribution of low-mass Pop III stars to the cosmic reionization. We can determine the total mass of the ionized region in a galaxy using Equation (19) and (25);

where $n_{\rm p}$ and $m_{\rm p}$ are the proton number density and mass. If we consider the ionized fraction inside a Milky Way-like galaxy, the ionized mass by low-mass stars accounts for only $\sim 10^{-9}-10^{-8}$ of the total baryon mass of a galaxy $\sim(0.8-1)\times 10^{11}M_{\odot}$ (Bland-Hawthorn & Gerhard, 2016), while for massive stars $M_{\rm ion}=1.1\times 10^{7}M_{\odot}$ which accounts for $\sim 10^{-4}$ of the total mass.