Thermal model

metal-poor disk evolution by performing hydrodynamic simulations, but the duration of their simulations is too short (∼100 yr after the protostar formation) to discuss how the accreting history affects the stellar evolution. In order to investigate SMS formation in low-metallicity environments in the future, we first need to examine the relation between the disk and stellar evolution in the case of low-mass star formation, for which previous studies exist for the solar metallicity cases. A knowledge from this examination will help us understand phenomena in the disk in the case of metal-enriched SMS formation, such as the disk fragmentation, thermal processes, and time fluctuation of accretion rate. In this chapter, we carry out two-dimensional hydrodynamic simulations of low-mass star formation to observe the disk evolution for longer term than in the previous studies and to evaluate the effect of thermal processes on the disk evolution in each metallicity range.

We also discuss the stellar growth based on the obtained accretion history. Finally, we infer the possible changes in the picture of SMS formation in low-metallicity environments when we apply our results to SMS formation. In Section 3.2, we describe our simulation model and the initial conditions. We then present the time evolution of the disk and envelope structures, as well as the impact of clump formation on protostar evolution, in Section3.3. Summary and discussion are given in Section3.4.

3.2 Method | 53 following expression

Λ_net = Λ_cont+ Λ_H2 + Λ_HD+ Λ_metal+ Λ_chem , (3.1) where Λ_cont is the gas and dust continuum cooling rate, Λ_H2 the H₂-line cooling rate, Λ_HD the HD-line cooling rate, Λ_metal the metal line cooling rate, and Λ_chem the chemical

cooling rate.

The net rate of continuum cooling is

Λ_cont= 4π(η−ξ_aJ) , (3.2)

where η is the emission coefficient, ξ_a is the absorption coefficient, and J is the mean intensity. The emission coefficient is given by

η= σ_SBρ

π κ_P,g(T_g)T_g⁴+κ_P,d(T_d)T_d⁴

, (3.3)

where σ_SB is the Stefan-Boltzmann constant, ρ is the mass density related to the surface density as shown in Equation (2.5), κ_P,g is the Planck mean opacity of gas, κ_P,d is the Planck mean opacity of dust,T_g is the gas temperature, andT_d is the dust temperature.

We calculate the Planck mean opacities using the tables fromMayer & Duschl (2005) for the gas andSemenov et al.(2003) for the dust. The absorption coefficient is written as ξ_a = (κ_P,g(T_g) +κ_P,d(T_d))ρ . (3.4) The mean intensity is (Omukai 2001)

J = 1

1 +x(B(T_irr) +xS) , (3.5)

whereB(T) is the Planck function, T_irr is the irradiation temperature, S is the source function, andx is a function that smoothly connects the optically thin and thick limits, written as (Tanaka & Omukai 2014)

x=τ_P+3

4τ_Pτ_R , (3.6)

with the Planck and Rosseland mean optical depths τP and τR, respectively. Those optical depths are calculated as

τ_{P(or R)} = 1

2 κ_{P(or R),g}(T_g) +κ_{P(or R),d}(T_d)

Σ , (3.7)

The Rosseland mean opacities of gas κ_R,g and dust κ_R,d are obtained from the same references of the Planck mean ones. The Planck function is

B(T) = σ_SB

π T⁴ , (3.8)

where T is the temperature of irradiation or dust (see Equation 3.10). The source function is defined as the ratio of the emission coefficient to the absorption coefficient:

S = η

ξ_a . (3.9)

In order to get the dust temperature, we consider the energy balance on dust grains due to the thermal emission, absorption, and gas-dust collision (Omukai et al. 2010):

κP,d(Td)B(Td) =κP,dJ+ Γcoll , (3.10) where Γ_coll is the dust heating rate by collisions with gas particles and its rate is (Hollenbach & McKee 1979)

Γ_coll = 4.4×10⁻⁶(f /ρ)_dustn_H

T_g 1000 K

1/2

(T_g−T_d) , (3.11) where (f /ρ)_dust is total volume of the dust per unit gas mass taken fromPollack et al.

(1994) andn_H is the number density of hydrogen nuclei, given by Equation (2.19).

The H₂ line cooling rate is given by using Equation (2.12). The HD-line cooling rate is calculated in the similar way as the H₂-line cooling rate. We multiply the cooling rate in the optically thin regime Λ_HD,thin (Flower et al. 2000) by the line-averaged escape probability β¯_esc,HD (Vorobyov et al. 2020b):

Λ_HD = ¯β_esc,HDΛ_HD,thine^−τ , (3.12)

where τ is the effective optical depth for continuum radiation which is given by Equa-tion (2.13).

In this thesis, the metal line cooling is due to the atomic fine-structure line emission of CII and OI, and its rate Λ_metal is given by summing the two contributions:

Λ_metal = Λ_CII+ Λ_OI . (3.13)

In order to calculate the line cooling rate of CII and OI, we model CII as a two level system and OI as a three level system and calculate the level populations from the

3.2 Method | 55 statistical balance among the levels. The CII and OI line cooling rates are given by

Λ_CII(OI) =y_CII(OI)n_HZ

×X

ul

hν_ulβ_esc,ulA_ulf_uS(ν_ul)−B_νul(T_irr)

S(ν_ul) , (3.14)

where the chemical fractions of CII and OI arey_CII = 9.27×10⁻⁵ andy_OI = 3.568×10⁻⁴, Z is the metallicity relative to the solar one, hν_ul is the energy difference between the upper level u and lower level l,β_esc,ul the escape probability of the line photon from the upper to the lower levels, A_ul the spontaneous radiative decay rate, f_u the occupancy of the upper level, S(ν_ul) the source function, andB_νul(T_irr) the Planck function. We take Z ranging from 0 (no metal) to 1Z_⊙(solar metal). The level energy and the spontaneous radiative decay rate are referred from Hollenbach & McKee (1989). The line escape probability is

β_esc,ul=

1−e^−τul τ_ul

e^−τ , (3.15)

where τ_ul is the optical depth for line emission from upper level to lower level, given by (see AppendixC)

τ_ul = c³ 8π^3/2ν_ul³ A_ul

g_u

g_lf_u−f_l

N_CII(OI)

v_th , (3.16)

where g_u,l is the statistical weight taken from Hollenbach & McKee(1989),N_CII is the column density of CII (OI), and v_th is the thermal velocity. The column density is

N_CII(OI) = 2H_gn_Hy_CII(OI) . (3.17)

The thermal velocity is

v_th= s

2k_BT_g

µm_H , (3.18)

where k_B is the Boltzmann constant, µ is the mean molecular weight, and m_H is the mass of a hydrogen nucleus. The source function for the line emission is

S(ν_ul) = 2hν_ul³ c²

g_u g_l

f_l f_u −1

−1

. (3.19)

The Planck function is

B_νul(T) = 2hν_ul³ c²

exp

hν_ul kBT

−1 ⁻¹

. (3.20)

The chemical cooling/heating processes include H ionization/recombination and H₂ dissociation/formation. Although we also include H⁻ detachment/attachment in Section 2.2.2, we do not take into account the cooling/heating associated with H⁻ reactions. The cooling/heating rate associated with H⁻ reactions are negligible because the chemical fraction of H⁻ is always smaller than those of other species. In this section, we calculate the chemical cooling rate as follows:

Λ_chem=

dy(H⁺)

dt χ_H−dy(H₂) dt χ_H2

n_H . (3.21)

When we calculate the irradiation temperature T_irr, we use Equation (2.24). In the simulation of the SMS formation, we compute the stellar luminosity L∗ by using the analytical expression in Hosokawa et al. (2012). In this study, we need to consider the stellar evolution to match the value of metallicity. The stellar luminosityL∗ is the sum of the accretion and stellar photospheric luminosities. The stellar photospheric luminosity is calculated by using tabulated stellar evolution tracks (Fukushima et al. 2020). The stellar luminosity and radius are provided from the tracks as functions of the stellar mass and accretion rate. Fukushima et al. (2020) have utilized the tracks in the case of the metallicities Z = 0,10⁻³, and 10⁻² Z_⊙. The tracks in other metallicity cases are obtained through the private communication with Hajime Fukushima.