Gravitational stability and fragmentation condition for discs around accreting supermassive stars

TL;DR

This study investigates the gravitational stability of accretion discs around supermassive stars, identifying conditions under which disc fragmentation occurs, which impacts the formation of supermassive black holes in the early universe.

Contribution

It provides a detailed analysis of the fragmentation conditions for accretion discs around supermassive stars using a steady thin disc model considering chemical and thermal processes.

Findings

Atomic flows require accretion rates above 10^{-1} M_sun/yr for fragmentation.

Molecular flows have a critical radius beyond which discs become unstable.

Disc fragmentation is likely during supermassive star formation.

Abstract

Supermassive stars (SMSs) with mass $\sim10^{5}~\rm{M}_{\odot}$ are promising candidates for the origin of supermassive black holes observed at redshift $\gtrsim6$. They are supposed to form as a result of rapid accretion of primordial gas, although it can be obstructed by the time variation caused by circum-stellar disc fragmentation due to gravitational instability. To assess the occurrence of fragmentation, we study the structure of marginally gravitationally unstable accretion discs, by using a steady one-dimensional thin disc model with detailed treatment of chemical and thermal processes. Motivated by two SMS formation scenarios, i.e., those with strong ultraviolet radiation background or with large velocity difference between the baryon and the dark matter, we consider two types of flows, i.e., atomic and molecular flows, respectively, for a wide range of the central stellar mass $10-10^5~\rm{M}_{\odot}$ and the accretion rate $10^{-3}-1~\rm{M}_{\odot}~\rm{yr}^{-1}$. In the case of a mostly atomic gas flowing to the disc outer boundary, the fragmentation condition is expressed as the accretion rate being higher than the critical value of $10^{-1}~\rm{M}_{\odot}~\rm{yr}^{-1}$ regardless of the central stellar mass. On the other hand, in the case of molecular flows, there is a critical disc radius outside of which the disc becomes unstable. Those conditions appears to be marginally satisfied according to numerical simulations, suggesting that disc fragmentation can be common during SMS formation.

Tables

Table 1. Table 1: Continuum absorption processes

Number	Name	Process	Reference
1	$H$ bound-free	$H+\gamma \to H^{+}+e$	Rybicki & Lightman (1979)
2	$H^{-}$ bound-free	$H^{-}+\gamma \to H+e$	John (1988)
3	$H^{-}$ free-free	$H+e+\gamma \to H+e$	John (1988)
4	$H$ free-free	$H^{+}+e+\gamma \to H^{+}+e$	John (1988)
5	$H_2$-$H_2$ CIA	$H_2$(v,J)$+H_2+\gamma \to H_2$(v’,J’)$+H_2$	Borysow et al. (1997)
6	$H_2$-$\mathrm{He}$ CIA	$H_2$(v,J)$+\mathrm{He}+\gamma \to H_2$(v’,J’)$+\mathrm{He}$	Borysow et al. (1997)
7	$H$ Rayleigh	$H+\gamma \to H+{\gamma }^{\prime }$	Kurucz (1970)
8	$H_2$ Rayleigh	$H_2+\gamma \to H_2+{\gamma }^{\prime }$	Dalgarno & Williams (1962)
9	Thomson	$e+\gamma \to e+{\gamma }^{\prime }$	Rybicki & Lightman (1979)

Table 2. Table 2: Coefficients for fitting function ( 30 )

process	$a_0$	$a_1$	$a_2$	$a_3$	$a_4$
$H$  free-bound	-22.410	-0.51833	2.3058d-3	0.0000	0.0000
$H^{-}$ free-bound	-31.211	1.1678	-3.0473d-2	0.0000	0.0000
$H^{-}$ free-free	-38.564	3.4681	-0.17880	0.0000	0.0000
$H$ free-free	-28.030	0.54235	-5.3291d-3	0.0000	0.0000
$H_2$-$H_2$ CIE	75.232	-1.4487d+2	65.126	-12.304	0.84166
$H_2$-$\mathrm{He}$ CIE	87.458	-1.6412d+2	75.157	-14.482	1.0120

Table 3. Table 3: Planck mean opacity log  κ P subscript 𝜅 P \kappa_{\mathrm{P}}

						$\log \rho$
$\log T$	-17.0	-16.0	-15.0	-14.0	-13.0	-12.0	-11.0	-10.0	-9.0	-8.0	-7.0	-6.0
2.7	-14.928	-13.898	-12.898	-11.898	-10.898	-9.898	-8.898	-7.898	-6.898	-5.898	-4.898	-3.898
2.8	-14.934	-13.905	-12.905	-11.905	-10.905	-9.905	-8.905	-7.905	-6.905	-5.905	-4.905	-3.905
2.9	-14.939	-13.904	-12.904	-11.904	-10.904	-9.904	-8.904	-7.904	-6.904	-5.904	-4.904	-3.904
3.0	-14.856	-13.815	-12.814	-11.814	-10.814	-9.814	-8.814	-7.814	-6.814	-5.814	-4.814	-3.814
3.1	-15.181	-13.858	-12.776	-11.750	-10.742	-9.739	-8.738	-7.738	-6.738	-5.738	-4.738	-3.738
3.2	-20.250	-16.670	-14.552	-12.483	-10.944	-9.758	-8.698	-7.679	-6.673	-5.672	-4.671	-3.671
3.3	-16.156	-15.656	-15.155	-14.611	-13.486	-11.469	-9.422	-7.903	-6.724	-5.667	-4.649	-3.643
3.4	-12.847	-12.347	-11.847	-11.347	-10.847	-10.347	-9.846	-9.262	-7.760	-6.078	-4.815	-3.730
3.5	-10.161	-9.727	-9.251	-8.759	-8.261	-7.762	-7.263	-6.763	-6.266	-5.793	-5.292	-4.094
3.6	-6.682	-6.646	-6.576	-6.413	-6.115	-5.705	-5.238	-4.749	-4.252	-3.755	-3.266	-2.851
3.7	-4.209	-3.707	-3.521	-3.456	-3.421	-3.370	-3.248	-3.002	-2.628	-2.177	-1.694	-1.207
3.8	-4.495	-3.494	-2.515	-1.665	-1.182	-1.002	-0.937	-0.896	-0.825	-0.664	-0.368	0.039
3.9	-4.942	-3.940	-2.940	-1.941	-0.949	-0.023	0.617	0.888	0.983	1.033	1.104	1.258
4.0	-5.388	-4.387	-3.387	-2.387	-1.387	-0.388	0.602	1.516	2.126	2.377	2.470	2.529
4.1	-5.835	-4.833	-3.833	-2.833	-1.833	-0.833	0.166	1.164	2.137	2.950	3.383	3.544
4.2	-6.280	-5.279	-4.279	-3.279	-2.279	-1.279	-0.279	0.721	1.719	2.705	3.587	4.128
4.3	-6.724	-5.724	-4.724	-3.724	-2.724	-1.724	-0.724	0.276	1.276	2.275	3.261	4.152

Table 4. Table 4: Rosseland mean opacity log  κ R subscript 𝜅 R \kappa_{\mathrm{R}}

						$\log \rho$
$\log T$	-17.0	-16.0	-15.0	-14.0	-13.0	-12.0	-11.0	-10.0	-9.0	-8.0	-7.0	-6.0
2.7	-11.152	-11.122	-10.938	-10.465	-9.823	-9.115	-8.370	-7.688	-7.113	-6.605	-6.041	-5.222
2.8	-10.753	-10.740	-10.638	-10.259	-9.636	-8.925	-8.189	-7.531	-7.002	-6.513	-5.839	-4.918
2.9	-10.354	-10.348	-10.292	-10.012	-9.436	-8.734	-8.009	-7.367	-6.861	-6.340	-5.552	-4.608
3.0	-9.954	-9.951	-9.922	-9.733	-9.225	-8.536	-7.816	-7.174	-6.666	-6.066	-5.220	-4.364
3.1	-9.491	-9.529	-9.529	-9.412	-8.984	-8.321	-7.611	-6.976	-6.458	-5.805	-5.022	-4.349
3.2	-8.764	-8.863	-8.944	-9.008	-8.798	-8.142	-7.410	-6.761	-6.216	-5.555	-4.917	-4.394
3.3	-7.320	-7.454	-7.587	-7.715	-7.822	-7.837	-7.419	-6.627	-5.992	-5.354	-4.848	-4.417
3.4	-5.907	-6.095	-6.267	-6.427	-6.576	-6.689	-6.696	-6.596	-6.118	-5.319	-4.796	-4.408
3.5	-4.143	-4.495	-4.803	-5.071	-5.305	-5.506	-5.635	-5.613	-5.491	-5.329	-4.812	-4.325
3.6	-2.061	-2.538	-2.997	-3.425	-3.811	-4.143	-4.389	-4.440	-4.211	-3.822	-3.385	-2.984
3.7	-0.606	-0.856	-1.260	-1.721	-2.189	-2.630	-2.959	-2.969	-2.657	-2.220	-1.751	-1.273
3.8	-0.526	-0.527	-0.537	-0.611	-0.865	-1.259	-1.649	-1.777	-1.504	-1.066	-0.589	-0.099
3.9	-0.526	-0.526	-0.526	-0.526	-0.529	-0.556	-0.686	-0.832	-0.653	-0.227	0.253	0.745
4.0	-0.526	-0.526	-0.526	-0.526	-0.525	-0.521	-0.496	-0.407	-0.163	0.321	0.845	1.353
4.1	-0.526	-0.526	-0.525	-0.524	-0.524	-0.521	-0.505	-0.426	-0.164	0.421	1.147	1.754
4.2	-0.526	-0.526	-0.526	-0.524	-0.518	-0.513	-0.504	-0.455	-0.265	0.228	1.064	1.911
4.3	-0.526	-0.526	-0.526	-0.525	-0.521	-0.501	-0.480	-0.447	-0.312	0.094	0.878	1.832

Table 5. Table 5: line cooling parameters for Li I

					Temperature  $T$ [K]
log $\tilde{N}$(Li)		300	500	800	1000	2000	3000	5000	8000	10000
				log $L_0$ [erg cm3 s-1]
		-46.86	-36.02	-29.48	-27.22	-22.60	-21.03	-19.75	-18.99	-18.72
				log ${\overline{\beta }}_{\mathrm{esc}}{ℒ}_{\mathrm{LTE}}$ [erg s-1]
6.0		-34.78	-22.36	-15.38	-13.05	-8.40	-6.85	-5.60	-4.89	-4.66
7.0		-35.64	-23.23	-16.25	-13.92	-9.26	-7.71	-6.44	-5.64	-5.36
8.0		-36.64	-24.23	-17.25	-14.92	-10.26	-8.71	-7.40	-6.47	-6.13
9.0		-37.60	-25.19	-18.21	-15.88	-11.23	-9.67	-8.31	-7.27	-6.89
10.0		-38.28	-25.87	-18.88	-16.56	-11.90	-10.35	-9.05	-8.13	-7.78
11.0		-38.68	-26.27	-19.29	-16.96	-12.31	-10.75	-9.50	-8.73	-8.46
12.0		-38.73	-26.31	-19.33	-17.00	-12.35	-10.80	-9.55	-8.85	-8.63
13.0		-38.73	-26.32	-19.34	-17.01	-12.35	-10.80	-9.56	-8.87	-8.65
14.0		-38.73	-26.32	-19.34	-17.01	-12.35	-10.80	-9.56	-8.87	-8.65
15.0		-38.73	-26.32	-19.34	-17.01	-12.35	-10.80	-9.56	-8.87	-8.65

Table 6. Table 6: Chemical reactions

Number	Reaction	Rate coefficient of forward reaction (${\mathrm{cm}}^3{s}^{-1}$)	Reference
1	$H+e\rightleftharpoons H^{+}+2e$	$k_1=\exp [-3.271396786\times {10}^1$	Janev et al. (1987)
		$+1.35365560\times {10}^1\ln T_e$
		$-5.73932875\times {10}^0{\left(\ln T_e\right)}^2$
		$+1.56315498\times {10}^0{\left(\ln T_e\right)}^3$
		$-2.87705600\times {10}^{-1}{\left(\ln T_e\right)}^4$
		$+3.48255977\times {10}^{-2}{\left(\ln T_e\right)}^5$
		$-2.63197617\times {10}^{-3}{\left(\ln T_e\right)}^6$
		$+1.11954395\times {10}^{-4}{\left(\ln T_e\right)}^7$
		$-2.03914985\times {10}^{-6}{\left(\ln T_e\right)}^8]$
2	$H^{-}+H\rightleftharpoons H_2+e$	$k_2=1.3500\times {10}^{-9}(T^{9.8493\times {10}^{-2}}+3.2852\times {10}^{-1}{T}^{5.5610\times {10}^{-1}}$	Kreckel et al. (2010)
		$+2.7710\times {10}^{-7}{T}^{2.1826})/(1.0+6.1910\times {10}^{-3}{T}^{1.0461}$
		$+8.9712\times {10}^{-11}{T}^{3.0424})+3.2576\times {10}^{-14}{T}^{3.7741})$
3	$2H\rightleftharpoons H^{+}+e+H$	$k_3=1.2\times {10}^{-17}{T}^{1.2}\exp \left(-\frac{157800}{T}\right)$	Lenzuni et al. (1991)
4	$H_2+e\rightleftharpoons 2H+e$	$k_4=k_{4,H}^{1-a}{k}_{4,L}^a$
		$k_{4,H}=1.91\times {10}^{-9}{T}^{0.136}\exp \left(-53407.1/T\right)$	Trevisan & Tennyson (2002)
		$k_{4,L}=4.49\times {10}^{-9}{T}^{0.11}\exp \left(-101858/T\right)$	Trevisan & Tennyson (2002)
		$a={\left(1+n_H/n_{\mathrm{crit}}\right)}^{-1}$
		$n_{\mathrm{crit}}=y\left(H\right)/n_{\mathrm{crit}}\left(H\right)+2y\left(H_2\right)/n_{\mathrm{crit}}\left(H_2\right)+y\left(\mathrm{He}\right)/n_{\mathrm{crit}}\left(\mathrm{He}\right)$
		$\log \left(n_{\mathrm{crit}}\left(H\right)\right)=3-0.416\log \left(T/{10}^4\right)-0.372{\left[\log \left(T/{10}^4\right)\right]}^2$
		$\log \left(n_{\mathrm{crit}}\left(H_2\right)\right)=4.845-1.3\log \left(T/{10}^4\right)+1.62{\left[\log \left(T/{10}^4\right)\right]}^2$
		$\log \left(n_{\mathrm{crit}}\left(\mathrm{He}\right)\right)=5.0792\left[1-1.23\times {10}^{-5}\left(T-2000\right)\right]$
5	$H_2+H\rightleftharpoons 3H$	see the reference	Martin et al. (1996)
6	$H^{-}+H^{+}\rightleftharpoons 2H$	$k_6=2.4\times {10}^{-6}{T}^{-0.5}\left(1.0+T/20000\right)$	Croft et al. (1999)
7	$2H_2\rightleftharpoons 2H+H_2$	$k_7=k_{7,H}^{1-a}{k}_{7,L}^a$
		$k_{7,H}=1.3\times {10}^{-9}\exp \left(-53300/T\right)$	Shapiro & Kang (1987)
		$k_{7,L}=\frac{5.996\times {10}^{-30}{T}^{4.1881}}{{\left(1.0+6.761\times {10}^{-6}T\right)}^{5.6881}}\exp \left(-54657.4/T\right)$	Martin et al. (1998)
8	$H_2+\mathrm{He}\rightleftharpoons 2H+\mathrm{He}$	$k_8=k_{8,H}^{1-a}{k}_{8,L}^a$
		$k_{8,H}=\mathrm{dex}\left[-1.75\log T-2.729-23474/T\right]$	Dove et al. (1987)
		$k_{8,L}=\mathrm{dex}\left[3.801\log T-27.029-29487/T\right]$	Dove et al. (1987)
9	$H^{+}+e\rightleftharpoons H+\gamma$	$k_9=2.753\times {10}^{-14}{\left(315614/T\right)}^{1.5}$	Ferland et al. (1992)
		$\times {\left[1.0+{\left(115188/T\right)}^{0.407}\right]}^{-2.242}$
10	$H+e\rightleftharpoons H^{-}+\gamma$	$k_{10}=\mathrm{dex}[-17.845+0.762\log T$	Wishart (1979)
		$+0.1523{\left(\log T\right)}^2$
		$-0.03274{\left(\log T\right)}^3]$                 ()
		$=\mathrm{dex}[-16.4199+0.1998{\left(\log T\right)}^2$
		$-5.447\times {10}^{-3}{\left(\log T\right)}^4$
		$+4.0415\times {10}^{-5}{\left(\log T\right)}^6]$   ($T>6000K$)
Note.— The temperature $T_e$ is in eV.

Equations

$$\Omega =\sqrt{\frac{G{M}_{\ast }}{r^3}},$$

$$Q=\frac{c_{\mathrm{s}}\Omega }{\pi G\Sigma }\sim 1,$$

$$\Sigma =\frac{c_{\mathrm{s}}\Omega }{\pi G}.$$

$$c_{\mathrm{s}}=\sqrt{\frac{\gamma k_{\mathrm{B}}T}{\mu m_{\mathrm{H}}}},$$

$$\rho =\frac{\Sigma }{2H},$$

$$H=\frac{c_{\mathrm{s}}}{\Omega }.$$

$$v_r=-\frac{\dot{M}}{2\pi r\Sigma }.$$

$$\frac{\mathrm{d}{e}_{\mathrm{in}}}{\mathrm{d}t}=-p\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{1}{\rho }\right)+\Gamma -\Lambda ,$$

$$e_{\mathrm{in}}=\frac{1}{\gamma -1}\frac{k_{\mathrm{B}}T}{\mu m_{\mathrm{H}}}.$$

$$p=\frac{k_{\mathrm{B}}T}{\mu m_{\mathrm{H}}}\rho .$$

$${\Gamma }_{\mathrm{vis}}=\frac{9}{4}\nu {\Omega }^2,$$

$$\nu =\frac{\dot{M}}{3\pi \Sigma },$$

$${\Lambda }_{\mathrm{cont}}=\frac{{\Lambda }_{\mathrm{cont},\mathrm{thin}}}{1+{\tau }_{\mathrm{P}}+\frac{3}{4}{\tau }_{\mathrm{R}}{\tau }_{\mathrm{P}}},$$

$${\tau }_{\mathrm{P}(\mathrm{or}\mathrm{R})}=\frac{1}{2}\Sigma {\kappa }_{\mathrm{P}(\mathrm{or}\mathrm{R})},$$

$${\Lambda }_{{\mathrm{H}}_2}={\overline{\beta }}_{\mathrm{esc}}{\Lambda }_{{\mathrm{H}}_2,\mathrm{thin}}{\mathrm{e}}^{-\tau },$$

$${\Lambda }_{\mathrm{Li}}={\overline{\beta }}_{\mathrm{esc}}{\Lambda }_{\mathrm{Li},\mathrm{thin}}{\mathrm{e}}^{-\tau }.$$

$${\Lambda }_{\mathrm{chem}}=\left({\chi }_{\mathrm{H}}\frac{\mathrm{d}y({\mathrm{H}}^{+})}{\mathrm{d}t}-{\chi }_{{\mathrm{H}}_2}\frac{\mathrm{d}y({\mathrm{H}}_2)}{\mathrm{d}t}\right)\frac{n_{\mathrm{H}}}{\rho },$$

$$y(i)=\frac{n(i)}{n_{\mathrm{H}}},$$

$$n_{\mathrm{H}}=\frac{\rho }{(1+4y_{\mathrm{He}})m_{\mathrm{H}}},$$

$$\alpha =\frac{\nu \Omega }{c_{\mathrm{s}}^2}=\frac{G}{3\pi }\frac{\dot{M}}{c_{\mathrm{s}}^3},$$

$${\Gamma }_{\mathrm{vis}}=\frac{3}{8\pi }\frac{{\Omega }^3}{\rho c_{\mathrm{s}}}\dot{M}.$$

$${\Lambda }_{{\mathrm{H}}^{-}\mathrm{fb}}$$

$$=\frac{4\pi {\Omega }^3{f}_{{\mathrm{H}}^{-}\mathrm{fb}}(T)K(T{)}^{-\frac{1}{2}}}{{\left[2\pi (1+4y_{\mathrm{He}})m_{\mathrm{H}}G\right]}^{\frac{3}{2}}\rho },$$

$$y(\mathrm{e})=1/\sqrt{n_{\mathrm{H}}K(T)},$$

$$f_{{\mathrm{H}}^{-}\mathrm{fb}}(T)K(T{)}^{-\frac{1}{2}}{c}_{\mathrm{s}}\propto \dot{M},$$

$${\dot{M}}_{\mathrm{crit}}=0.1{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}{\left(\frac{T}{3000\mathrm{K}}\right)}^{\frac{3}{2}}.$$

$$\dot{M}=\varphi \frac{c_{\mathrm{s},\mathrm{col}}^3}{G}=0.1\varphi {\left(\frac{c_{\mathrm{s},\mathrm{col}}}{8\mathrm{km}{\mathrm{s}}^{-1}}\right)}^3{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}.$$

$${\eta }_{\nu }=\frac{2h{\nu }^3}{c^2}\frac{z_{\mathrm{B}}}{z_{\mathrm{A}}{z}_{\mathrm{e}}}$$

$$\times \exp \left[-\frac{h\nu -{\chi }_B}{k_{\mathrm{B}}T}\right]n(A)n(\mathrm{e});$$

$${\eta }_{\nu }$$

$${\Lambda }_{i,\mathrm{thin}}=\frac{4\pi }{\rho }\int {\eta }_{i,\nu }\mathrm{d}\nu =\frac{4\pi }{\rho }{f}_i(T)n(\mathrm{A})n(\mathrm{B}),$$

Full text

††pubyear: 2018††pagerange: Gravitational stability and fragmentation condition for discs around accreting supermassive stars–C

1 Introduction

Discoveries of dozens of quasars at high redshift $z\gtrsim 6.5$ demonstrate the very existence of supermassive black holes (SMBHs) as massive as $\sim 10^{9}~{}\mathrm{M}_{\odot}$ in less than a billion years after the Big Bang (e.g., Venemans et al. 2013; Matsuoka et al. 2018; see also Gallerani et al. 2017 for a review), including the most distant one with $8\times 10^{8}~{}\mathrm{M}_{\odot}$ at redshift $z=7.5$ (Bañados et al., 2018). Although the origin of those BHs is still a mystery, short available time for their growth seems to favor massive seed BHs (see, e.g., Volonteri 2012; Haiman 2013 for a review).

Remnant BHs of Population III stars with $\sim 10^{2}~{}\mathrm{M}_{\odot}$ are one of the candidates for such seeds (e.g., Madau & Rees, 2001). They can grow SMBHs in the available time if either with continuous accretion at the Eddington limit or with short episode of super-Eddington growth. In reality, however, the BH’s growth can easily be hindered by its own radiative feedback (Milosavljević et al., 2009; Park & Ricotti, 2011; Sugimura et al., 2018; Orofino et al., 2018). An alternative and attractive pathway for seed BH formation is via the so-called direct collapse (e.g., Bromm & Loeb 2003), where a supermassive star (SMS) with $\sim 10^{5}~{}\mathrm{M}_{\odot}$ collapses into a BH with a similar mass by general relativistic instability (Umeda et al., 2016; Woods et al., 2017; Haemmerlé et al., 2018). Here, SMSs are supposed to form from a primordial gas in some peculiar sites. Unlike in ordinary first star formation, which is driven by $\mathrm{H}_{2}$ cooling (see, e.g., Bromm & Larson 2004; Glover 2013 for a review), $\mathrm{H}_{2}$ is dissociated by strong external far-ultraviolet (UV) radiation from nearby galaxies and the contraction of clouds is solely caused by $\mathrm{H}$ atomic cooling in the most intensively studied channel for SMS formation (Omukai, 2001; Omukai et al., 2008; Shang et al., 2010; Regan et al., 2014; Sugimura et al., 2014). Such clouds contract almost isothermally at $\sim 10^{4}$ K without experiencing vigorous fragmentation. The protostar formed at the centre subsequently accretes the gas at a high rate of $10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ due to high temperature (Latif et al., 2013; Inayoshi et al., 2014; Becerra et al., 2015; Chon et al., 2018). Such rapidly accreting protostar inflates greatly in radius with effective temperature of several 1000 K and grows supermassive with $\gtrsim 10^{5}~{}\mathrm{M}_{\odot}$ avoiding ionization radiation feedback on the accretion flow (Hosokawa et al., 2012; Hosokawa et al., 2013), before collapsing by general relativistic instability.

In addition to the above scenario for SMS formation, Hirano et al. (2017) recently proposed another channel where the suppression of $\mathrm{H}_{2}$ cooling is not always required. They considered a pristine cloud in a halo with 3$\sigma$ higher than the average velocity difference between the dark matter and the baryon generated before the recombination epoch at $z\gtrsim 1100$ (Tseliakhovich & Hirata, 2010; Fialkov et al., 2012; Schauer et al., 2017). As a result of enhanced effective Jeans mass by such streaming motions, the formed protostar in the halo accretes a gas with a high rate of $\sim 10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$. The star eventually grows to $4\times 10^{4}~{}\mathrm{M}_{\odot}$ in their calculation.

Protostellar ionization feedback is suppressed as long as the accretion rate higher than $10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ is maintained (Hosokawa et al., 2012; Hosokawa et al., 2013). Once the accretion rate falls below this value for sometime, however, the radiation feedback may become significant and eventually terminate the stellar growth. Adopting a simple sporadic accretion history with repeating burst and quiescent phases, Sakurai et al. (2015) showed that the protostar contracts to the main sequence and ionization feedback becomes active if the quiescent phase exceeds $\sim 10^{3}$ yrs. Such time variations would be caused by clump formation in the circum-(proto)stellar disc, which is formed due to the finite angular momentum of the cloud. In fact, numerical simulations of SMS formation (Chon et al., 2018; Hirano et al., 2017) found a sign of disc fragmentation by gravitational instability and resulting accretion variation although their spatial resolution and chemical/thermal modeling are not good enough to draw definitive conclusions.

Inayoshi & Haiman (2014) and Latif & Schleicher (2015) studied the gravitational stability of circum-stellar discs in the process of SMS formation by way of the steady thin disc model (Shakura & Sunyaev, 1973), and claimed that the disc becomes unstable and fragmentation occurs for a disc accreting at the typical rate of $10^{-1}\ \mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$. In their calculation, however, simplified treatment of thermal modeling is adopted without solving detailed chemical reactions: Latif & Schleicher (2015) assumed that the accretion flow is composed of $\mathrm{H}_{2}$ at the outer boundary and becomes atomic $\mathrm{H}$ once the temperature exceeds 4000 K, while Inayoshi & Haiman (2014) assumed that the gas is always in the atomic state and the temperature is determined by the balance between the viscous heating and the $\mathrm{H}^{-}$ free-bound emission cooling. Their conclusion needs to be justified by more detailed modeling.

In this paper, we construct models for one-dimensional steady accretion discs by solving non-equilibrium chemical and thermal evolution and investigate their gravitational stability. While the works by Inayoshi & Haiman (2014) and Latif & Schleicher (2015) are limited to a typical accretion rate of $10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$, here we also study the cases with different values of the accretion rate and central stellar mass.

The paper is organized as follows. We describe our disc model and stability criterion in Section 2. We then present the obtained disc structures and the result for the stability analysis in Section 3. Discussion and conclusion are given in Section 4. Appendices are for the details of radiative and chemical processes.

2 Model

2.1 Basic Framework

We construct the one-dimensional axisymmetric and steady model for an accretion disc feeding the central star under a given central stellar mass and accretion rate. The disc is assumed to rotate at the Keplerian angular velocity

[TABLE]

where $r$, $G$ and $M_{\ast}$ are the radius, the gravitational constant and the central stellar mass, respectively. We assume the disc is marginally gravitationally unstable and the angular momentum is transported by the spiral arms. Thermal and chemical evolution in the flow is calculated as in Sections 2.2 and 2.3. We describe the angular momentum transport by the $\alpha$-viscosity prescription (Shakura & Sunyaev, 1973) and judge disc fragmentation from the required value of $\alpha$-parameter (see Section 2.4).

For a marginally unstable disc, Toomre’s $Q$-value (Toomre, 1964) is

[TABLE]

where $c_{\mathrm{s}}$ is the sound velocity and $\Sigma$ is the surface density, as a result of the angular momentum transfer by the spiral arms (Vorobyov & Basu, 2007; Tsukamoto et al., 2015). Another possible mechanism of the angular momentum transfer is via the turbulent viscosity. According to the results shown in Section 3, however, the actual $\alpha$-value required is always larger than the typical value for the turbulent viscosity ($\alpha\sim 10^{-2}$ for magnetorotational instability; e.g., Bai & Stone 2013). This means that the turbulent viscosity alone cannot maintain the disc in a gravitationally stable state and the disc eventually settles in a marginally unstable state where the angular momentum is transferred due to the spiral arms. From the condition $Q=1$, the surface density is obtained as

[TABLE]

The sound velocity is

[TABLE]

where $\gamma$ is the specific heat ratio, $k_{\mathrm{B}}$ is the Boltzmann constant, $T$ is the gas temperature, $\mu$ is the mean molecular weight and $m_{\mathrm{H}}$ is the mass of a hydrogen nucleus.

The gas mass density is given by

[TABLE]

where the scale height $H$ characterizing the thickness of the disc is estimated from vertical hydrostatic balance:

[TABLE]

The radial velocity is related to the mass accretion rate $\dot{M}$ as:

[TABLE]

2.2 Thermal Evolution

We calculate the time evolution of internal energy $e_{\mathrm{in}}$ (per unit mass) in the flow by solving the energy equation,

[TABLE]

where $p$ is the pressure, $\Gamma$ is the heating rate per unit mass and $\Lambda$ is the cooling rate per unit mass. The internal energy and the temperature are related as

[TABLE]

The pressure is given by the ideal-gas equation of state,

[TABLE]

We consider viscous heating, continuum cooling, $\mathrm{H}_{2}$ molecular and Li atomic line cooling and chemical cooling as the heating and cooling processes.

The viscous heating rate per unit mass is

[TABLE]

where $\nu$ is the kinematic viscosity. For consistency with the assumption that the disc is steady and rotating at the Keplerian velocity, the kinematic viscosity must satisfy the relation

[TABLE]

from the angular momentum conservation. In our model, the kinematic viscosity is determined so that the surface density satisfies Equation (3) for a given accretion rate. The continuum radiation processes considered are $\mathrm{H}$ free-bound emission, $\mathrm{H}^{-}$ free-bound emission, $\mathrm{H}^{-}$ free-free emission, $\mathrm{H}$ free-free emission, $\mathrm{H}_{2}$-$\mathrm{H}_{2}$ collision induced emission (CIE) and $\mathrm{H}_{2}$-$\mathrm{He}$ CIE (reverse processes of 1-6 in Table 1). The continuum cooling rate in optically thin regime $\Lambda_{\mathrm{cont,thin}}$ is given by sum of Equation (29) in Appendix A.1 over all the processes. Smoothly connecting the optically thin and thick limits, the continuum cooing rate in both regimes is given as (Tanaka & Omukai, 2014)

[TABLE]

where $\tau_{\mathrm{P(or\ R)}}$ is the Planck (Rosseland) mean optical depth, given by

[TABLE]

with the Planck (Rosseland, respectively) mean opacity $\kappa_{\mathrm{P(R)}}$ calculated as in Appendix A.2. The $\mathrm{H}_{2}$-line cooling rate is given by multiplying the value in the optically thin regime $\Lambda_{\mathrm{H}_{2},\mathrm{thin}}$ with the line-averaged escape probability $\overline{\beta}_{\mathrm{esc}}$:

[TABLE]

where $\tau=\sqrt{\tau_{\mathrm{P}}\tau_{\mathrm{R}}}$ is the effective optical depth for continuum radiation. We use fitting functions for $\overline{\beta}_{\mathrm{esc}}$ from Fukushima et al. (2018) and for $\Lambda_{\mathrm{thin},\mathrm{H}_{2}}$ from Glover (2015), respectively. As chemical cooling/heating processes, we consider those associated with H ionization/recombination and H2 dissociation/formation. The Li cooling rate including line self-absorption effect $\overline{\beta}_{\mathrm{esc}}\Lambda_{\mathrm{Li,thin}}$ is given by Equation 36 (see Appendix B). Considering also the continuum absorption, we obtain the Li cooling rate

[TABLE]

The chemical cooling rate is

[TABLE]

where $\chi_{\mathrm{H}}=13.6$ eV and $\chi_{\mathrm{H}_{2}}=4.48$ eV are the binding energies. The chemical fraction of species $i$, $y(i)$, is defined by the number fraction relative to hydrogen nuclei:

[TABLE]

where $n(i)$ and $n_{\mathrm{H}}$ are the number densities of species $i$ and hydrogen nuclei, respectively. The latter is related to the mass density of the gas $\rho$ as

[TABLE]

where $y_{\mathrm{He}}=8.333\times 10^{-2}$ the number fraction of He relative to hydrogen nuclei.

We do not include radiative heating by the central star in our calculation. We estimated it under the assumption of optically-thin irradiation at the disc midplane and found it is generally lower than the viscous heating rate even at the Eddington luminosity because the coupling between the stellar radiation and disc is very inefficient due to the geometrical dilution of radiation ($\propto r^{-2}$) and transparency of the disc (typically $\tau\sim 10^{-3}$-$10^{-4}$ in the case of atomic flow and $\sim 10^{-6}$-$10^{-10}$ in the case of molecular flow at 1000 au, where the discs tend to be most gravitationally unstable; see Figs. 1c and 3c). The radiative heating does not change the disc structure except in the outer regions in the case of low accretion rate $\lesssim 10^{-2}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$. In this case, the stellar radiation heating can be comparable to the viscous heating and raises the temperature to $\sim$1000 K. Although this somewhat stabilizes this part of the disc, the disc remains unstable (see Section 2.4) as in the case without stellar radiation. We thus expect that the effect of stellar radiation heating does not change our conclusion on disc fragmentation.

2.3 Chemical Reactions

We take into account the 20 chemical reactions 111We include the collisional dissociation of $\mathrm{H}_{2}$ by atomic helium impact (reaction 8 in Table 6) for consistency with the expression used for the critical density of $\mathrm{H}_{2}$, where the helium impact is also considered. among the following five H-bearing species, $\mathrm{H}$, $\mathrm{H}_{2}$, $\mathrm{H}^{+}$, $\mathrm{H}^{-}$ and $\mathrm{e}$, summarized in Table 6 (see Appendix C for details). All the reactions are paired with their reverse reactions. We solve non-equilibrium kinetic equations for $\mathrm{H}$, $\mathrm{H}_{2}$, $\mathrm{H}^{+}$ and $\mathrm{e}$, while $\mathrm{H}^{-}$ fraction is calculated from the equilibrium between the forward and reverse reactions of reactions 2, 6 and 10 in Table 6. We assume $\mathrm{He}$ to be all neutral with fractional abundance $y_{\rm He}=8.333\times 10^{-2}$ as our temperature range is limited below 104 K. We also assume that Li remains neutral with $y_{\rm Li}=4.15\times 10^{-10}$ (Mayer & Duschl, 2005) since they are still dominantly atomic in the calculated range and the small amount of Li ions does not affect thermal evolution at all. $\mathrm{H}_{2}$ photodissociation and $\mathrm{H}^{-}$ photodetachment by external radiation are not included. Because the number density is larger than $10^{8}~{}\mathrm{cm}^{-3}$ in our calculation range and the shielding of far-UV radiation is effective due to the large column density, the collisional dissociation works more efficiently than the dissociation due to external far-UV radiation. The effect of $\mathrm{H}^{-}$ photodetachment is also negligible because $\mathrm{H_{2}}$ formation via the three body reaction is more effective than the $\mathrm{H}^{-}$ channel at such high densities as considered in this work.

2.4 Fragmentation Condition

We judge whether the disc fragments or not from the viscosity parameter

[TABLE]

where we have used Equations (3) and (12) for the last expression. We here adopt $\alpha>1$ as the fragmentation condition following the results of Zhu et al. (2012), who claimed that discs fragment when the $\alpha$-value corresponding to the angular momentum transfer by gravitational torque exceeds unity from two-dimensional numerical simulations of protoplanetary discs. In other words, there is an upper limit on the angular momentum transported by the spiral arms and if more matter is accumulated than the disc can extract its angular momentum, the disc will fragment. In contrast, if $\alpha<1$ throughout, the disc remains stable and the angular momentum is transported by the gravitational torque.

2.5 Model setup

We study cases with a range of the central stellar mass $10~{}\mathrm{M}_{\odot}\leq M_{\ast}\leq 10^{5}~{}\mathrm{M}_{\odot}$ and the accretion rate $10^{-3}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}\leq\dot{M}\leq 1~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$. The disc structures at different central stellar masses can be regarded as temporal evolution around a growing SMS. Recall that the typical accretion rate is $10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ for SMS formation and $10^{-3}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ for ordinary first star formation, respectively. The calculated range of accretion rate encompasses these values.

We set the outer boundary of the disc at 1000 au and assigning the temperature and chemical fractions there. We evaluate the surface density by Equation (3) using the temperature and the central stellar mass, and calculate the radial velocity using Equation (7). Equation (12) gives the kinematic viscosity that is required to maintain such a disc structure. We then calculate the heating, cooling and chemical reaction rate, and update the temperature and the chemical fractions along the inflow. We continue the radial integration until the temperature exceeds $10^{4}~{}\mathrm{K}$, above which $\alpha<0.4$ for $\dot{M}<1~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ (Equation 20) and thus the disc is always stable.

We consider two different conditions in the flow at outer boundary: the atomic flow with temperature $T=3000$ K, molecular fraction $y({\rm H_{2}})=10^{-8}$, and ionization degree $y({\rm H^{+}})=10^{-6}$ and the molecular flow with $T=1000$ K, $y({\rm H_{2}})=0.495$, and $y({\rm H^{+}})=10^{-10}$. Since ${\rm H^{-}}$ fraction is much lower than those of other species, the ${\rm H}$ and ${\rm e}$ fractions are calculated by y(\mathrm{H})=1-2y(\mathrm{H}_{2})-y(\mathrm{H}^{+})\ and $y(\mathrm{H}^{+})=y(\mathrm{e})$, respectively, at the outer boundary. The atomic flow corresponds to an SMS forming in an atomically cooling cloud as a result of $\mathrm{H}_{2}$ photodissociation by external radiation, while the molecular one to that forming in an $\mathrm{H}_{2}$ cloud due, e.g., to baryonic streaming motions. For the former, we take the temperature and ionization degree from Omukai (2001) for the far-UV-irradiated cloud at number density $n_{\mathrm{H}}=10^{10}~{}\mathrm{cm}^{-3}$ with the far-UV intensity $J_{21}=10^{5}$ (in unit of $10^{-21}$ erg $\mathrm{s^{-1}}~{}\mathrm{Hz^{-1}}~{}\mathrm{str^{-1}}~{}\mathrm{cm^{-2}}$) and for the latter we take the self-regulated disc model from Schleicher et al. (2016) at number density $n_{\mathrm{H}}=10^{10}~{}\mathrm{cm}^{-3}$.

3 Result

We present the disc structures and examine their gravitational stability for atomic flows in Section 3.1 and for molecular flows in Section 3.2, respectively.

3.1 Atomic Accretion Flows

3.1.1 Disc Structures

The disc structures for atomic flows accreting at $\dot{M}=10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ are shown in Fig. 1 for the central stellar mass of $M_{\ast}=10$, $10^{2}$, $10^{3}$, $10^{4}$ and $10^{5}~{}\mathrm{M}_{\odot}$. In all cases, the temperature first jumps up to around 4000 K at which the cooling and heating rates balances at the outer boundary at 1000 au, and then gradually decreases as the gas flows inward (Fig. 1a). When the optical depth exceeds unity, the radiative cooling becomes ineffective and the temperature suddenly rises to $>10^{4}$ K (Figs. 1a and c). Just before the optical depth reaches unity, the temperature takes its minimum at 3000 K, which does not depend on the central stellar mass. The density depends on the radius as $n_{\mathrm{H}}\propto r^{-3}$ (Fig. 1b), as can be derived analytically from Equations (1), (3), (5) and (6).

To understand the temperature evolution in Fig. 1a, we show in Fig. 2 the radial distribution of the chemical fractions (Fig. 2a) and the heating and the cooling rates (Fig. 2b) for the case with $M_{\ast}=10^{2}~{}\mathrm{M}_{\odot}$ and $\dot{M}=10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ (corresponding to the blue lines in Fig. 1). Dominant chemical species is atomic $\mathrm{H}$ from 1000 au to 10 au, where the ionization starts and the temperature increases (Fig. 2a), and dominant heating and cooling processes are the viscous heating and the $\mathrm{H}^{-}$ free-bound emission cooling, respectively, and the temperature is determined by their balance (Fig. 2b). This is also true for all the cases shown in Fig. 1. Note that the contribution of lithium line cooling is always negligible although Mayer & Duschl (2005) insisted lithium cooling can be important, if it is at LTE value, at the relevant density and temperature ranges (see their Fig. 4). This is because the actual rate is far smaller than the LTE rate as the electron number density is generally several orders of magnitude below the critical density, even with the photon trapping effect. In the inner region, although the electron density increases, the continuum optical depth exceeds unity before reaching the critical density and the cooling rate decreases exponentially thereafter.

The minimum temperature is set by the balance between the viscous heating (from Equations 6, 5 and 12),

[TABLE]

and the cooling rate by $\mathrm{H}^{-}$ free-bound emission (from Equations 13 and 29)

[TABLE]

where $f_{\mathrm{H}^{-}\mathrm{fb}}(T)$ is a function of temperature (see Equation 29) and we have assumed that mostly atomic gas of $y(\mathrm{H})\sim 1$ (Fig. 2a) with small ionization degree given by the equilibrium value

[TABLE]

due to high enough density ($\sim 10^{16}~{}\mathrm{cm}^{-3}$) at the temperature minimum (Fig. 1a and b), where the equilibrium constant $K(T)$ is again a function only of temperature (see Appendix C). We use the optically thin rate as the optical depth is still small at this epoch (Fig. 1c). Equating Equations (21) and (22), we obtain

[TABLE]

which implicitly gives the temperature $T$ as a function of $\dot{M}$. For example, for $\dot{M}=10^{-1}\ \mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$, the equilibrium temperature is $T=$3100 K, consistent with that in Fig. 1a. The slight outward temperature increase in Fig. 1 is due to the smaller actual cooling rate than that by Equation (22) as a result of smaller electron fraction than the equilibrium value used above. It should be noted that the minimum temperature in the disc given as a solution of Equation (24) depends only on the accretion rate and not on the central stellar mass, as seen in Fig. 1a.

Next we study the dependence on the accretion rate. The disc structures for the central stellar mass $M_{\ast}=10^{2}~{}\mathrm{M}_{\odot}$ are shown in Fig. 3 for four different accretion rates $\dot{M}=10^{-3}$, $10^{-2}$, $10^{-1}$ and $1~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$. The temperature profiles in the high accretion rate cases of $10^{-1}$ and $1~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ behave in a similar way as seen before in Fig. 1: the gradual inward decrease from 4000 K to 3000 K followed by the sudden jump when the disc becomes optically thick. At lower accretion rates, i.e., $10^{-3}$ and $10^{-2}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$, however, the temperature experiences a steep drop to 1000-2000K at some radius before eventually jumping up to $>10^{4}$ K. This is a result of H2 formation and its CIE cooling. For illustration, we show the radial distribution of the chemical fractions (Fig. 4a), and the heating and cooling rates (Fig. 4b) for the case of $\dot{M}=10^{-2}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ (green lines in Fig. 3). At radii outer than 9 au, the gases are mostly atomic and the temperatures are set by the balance between the viscous heating and the cooling by the $\mathrm{H}^{-}$ free-bound emission. Once the temperature becomes as low as 3000 K, transition to the molecular phase occurs, accompanied with a sudden temperature drop to 2000K due to the resultant $\mathrm{H}_{2}$ CIE cooling. The transition radius moves outward to 100 au in the case of $\dot{M}=10^{-3}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ due to lower viscous heating and temperature. Note that this molecular transition in the disc has not been recognized before (Inayoshi & Haiman, 2014; Latif & Schleicher, 2015) and only becomes able to be followed by our detailed thermal and chemical modeling. Similar dependence of disc structures on the accretion rate is also observed for other central stellar masses.

3.1.2 Gravitational Instability of Discs

We now examine the disc instability for the atomic inflow using the $\alpha$ parameters. The radial distributions of the $\alpha$-value are shown for different stellar masses but with the same accretion rate of $\dot{M}=10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ in Fig. 1d, while those for the same mass $M_{\ast}=10^{2}~{}\mathrm{M}_{\odot}$ but with four different accretion rates are shown in Fig. 3d. Since $\alpha\propto\dot{M}/T^{3/2}$ (Equation 20) for constant $\mu$, the radial variation of $\alpha$ follows that of the temperature but in a upside-down way. As seen in Fig. 3d, $\alpha$ tends to have higher value for higher accretion rate. With high accretion rate of $\dot{M}\geq 10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$, $\alpha$ is almost constant radially in the optically thin part of the disc, while with lower rate $\dot{M}\leq 10^{-2}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$, the maximum value of $\alpha$ is reached at the the temperature minimum by the runaway $\mathrm{H}_{2}$ formation.

As seen in Section 3.1.1, the minimum temperature depends only on the accretion rate and not on the central stellar mass, so does the maximum value of $\alpha$, $\alpha_{\mathrm{max}}$ (see Fig. 1d). Fig. 5 shows $\alpha_{\mathrm{max}}$ as a function of the accretion rate. There are two branches where $\alpha_{\mathrm{max}}\propto\dot{M}$ below and above $\dot{M}=2\times 10^{-2}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$, which correspond, respectively, to the case where $\alpha_{\mathrm{max}}$ is attained in the $\mathrm{H}_{2}$ phase (as shown in Fig. 4) and the case where the gas always remains atomic in the optically thin part (Fig. 2). Here we adopt $\alpha>1$ as the condition for disc fragmentation (Section 2.4). As seen in Fig. 5, this occurs for $\dot{M}\ga 10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$. This critical accretion rate for disc fragmentation can be obtained by substituting $\alpha=1$ and the sound velocity for neutral gas with the minimum temperature of 3000 K (Fig. 1a) into Equation (20):

[TABLE]

Since the typical accretion rate for SMSs formation falls on around this critical value, whether fragmentation actually occurs or not would depend on detailed dynamics of collapse and accretion (see Section 4).

3.2 Molecular Accretion Flows

3.2.1 Disc Structures

The disc structure for a molecular flow is shown in Fig. 6 in the case with accretion rate $\dot{M}=10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ for different central stellar masses $M_{\ast}=10$, $10^{2}$, $10^{3}$, $10^{4}$ and $10^{5}~{}\mathrm{M}_{\odot}$. Soon after the beginning of calculation, the temperature near the outer boundary quickly converses to the thermal equilibrium value, which is higher for the higher accretion rate (Fig. 6a). In the highest mass case of $M_{\ast}=10^{5}~{}\mathrm{M}_{\odot}$, this temperature reaches $\sim 3000$ K and the gas becomes atomic immediately as a result of collisional H2 dissociation. At lower stellar masses $M_{\ast}\leq 10^{4}~{}\mathrm{M}_{\odot}$, the initial temperature at the outer boundary is lower than $\sim 2000$ K but increases as the gas flows inward, first gradually up to $\sim 2000$ K and then suddenly to $\sim 3000$ K. As an example, for the case of $M_{\ast}=10^{2}~{}\mathrm{M}_{\odot}$ (blue lines in Fig. 6), radial profiles of the chemical fractions (Fig. 7a) and the heating and cooling rates (Fig. 7b) are presented in Fig. 7. At 100 au, transition from the molecular to atomic phases proceeds due to the H2 collisional dissociation (Fig. 7a) and simultaneously the dominant cooling process shifts from the $\mathrm{H}_{2}$-line emission to the $\mathrm{H}^{-}$ free-bound emission (Fig. 7b), which causes the abrupt temperature increase to 3000 K. After the transition to $\mathrm{H}$, the temperature is set by the balance between the viscous heating and the $\mathrm{H}^{-}$ free-bound emission cooling and slightly decreases inward until the disc becomes optically thick, as already seen for the atomic flows (Fig. 1a). Note that the decrease of $\mathrm{H}_{2}$ fraction with density here looks more abrupt (Figs. 6a and 7a) than in the case of collapsing prestellar cores (e.g., Omukai, 2001). This difference can be attributable to smaller chemical cooling rate in our case: the latent heat associated with $\mathrm{H}_{2}$ dissociation is extracted in a longer timescale for density enhancement, i.e., disc viscous time, than the free-fall timescale of collapsing cores, and thus at a smaller chamical cooling rate. The temperature in the disc increases more easily, resulting in the sudden $\mathrm{H}_{2}$ dissociation.

Flows with different accretion rates $\dot{M}=10^{-3}$, $10^{-2}$, $10^{-1}$ and $1\ \mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ are shown in Fig. 8 for the same stellar mass $M_{\ast}=10^{2}~{}\mathrm{M}_{\odot}$. The dependence of the disc structure on the accretion rate is similar also in the cases with other values of stellar mass. In high accretion-rate cases of $10^{-1}$ and $1~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$, the disc structures are the same as discussed above in Fig. 6: the temperature suddenly increases from 2000 K to 3000 K at atomic phase transition, then stays at 3000 K in the optically thin atomic gas and eventually jumps up to $>10^{4}$ K when the disc becomes optically thick. With lower accretion rate of $\dot{M}=10^{-2}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$, the gas returns to molecular phase again at 10 au and temperature drops to 2000K as seen for the atomic accretion flow with low accretion rates (Fig. 3). In the case of $\dot{M}=10^{-3}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$, the temperature variation is more gradual than in the other cases. The chemical fractions and the heating and cooling rates in this case are shown in Fig. 9. Unlike in higher accretion rate cases, $\mathrm{H}_{2}$ now survives as the dominant species until the flow reaches the innermost part. The $\mathrm{H}$ fraction increases from 70 au. After peaking at 30 au, it then decreases inward as the temperature also decreases (Fig. 8). The temperature decrease at $<30$ au is caused by the cooling by $\mathrm{H}_{2}$ CIE (Fig. 9b). The temperature more inside is set by the thermal balance between this cooling and the viscous heating until the disc becomes optically thick to the $\mathrm{H}_{2}$ collision-induced absorption.

3.2.2 Disc Fragmentation

Next, we examine the disc fragmentation condition for the molecular accretion flows. The radial distributions of $\alpha$-value for such flows are presented in Fig. 6d for $\dot{M}=10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ at different stellar masses. The $\alpha$-value takes its maximum just inside the outer boundary where the temperature is the lowest, except in the case with $M_{\ast}=10^{5}~{}\mathrm{M}_{\odot}$ (see Section 3.2.1 for details). Going inward, the $\alpha$-value decreases along with the increasing temperature. The $\mathrm{H}_{2}$ dissociation and resultant temperature jump to 3000 K causes the corresponding drop of $\alpha$ below unity in all the cases. We thus expect that fragmentation, if any, occurs only in the outer $\mathrm{H}_{2}$ region. The cases for different accretion rates are shown for $M_{\ast}=10^{2}~{}\mathrm{M}_{\odot}$ in Fig. 8d. Again, the maximum $\alpha$ is reached around the outer boundary and its value is almost proportional to the accretion rate (Equation 20) as in the case of atomic flows.

As discussed above, for the molecular flows, the $\alpha$-value tends to be higher at outer radius, suggesting that if the disc is large enough the $\alpha$ exceeds unity at some radius. We can thus define the critical radius here as the radius where $\alpha=1$ is satisfied and plot it in Fig. 10 as a function of the accretion rate for each central stellar mass. If the size of a disc is smaller than this, the disc remains stable. Otherwise, its outer part would become gravitationally unstable and fragment. In all the cases shown in Fig. 10, the critical radius becomes abruptly small at $\sim 10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ because the condition of $\alpha<1$ is attained only when the gas ionizes at innermost radius for higher accretion rate. Also, the critical radius increases with the central stellar mass, reflecting that the temperature profiles shift toward outer radius as seen in Fig. 6a. With increasing stellar mass by accretion, the disc at a given radius becomes more stable and only the more outer part can be gravitationally unstable.

4 Summary and Discussion

Seed black hole (BH) formation by the direct collapse and its subsequent growth is one viable scenario for the supermassive black hole (SMBH) formation in the early Universe. In this framework, a small protostar is first formed by the gravitational collapse of a cloud, and grows to supermassive star (SMS) by rapid accretion, and finally the SMS collapses by general relativistic instability into a seed BH. In this paper, we have investigated gravitational instability of the discs around rapidly accreting protostars with mass $10-10^{5}~{}\mathrm{M}_{\odot}$ at rate $10^{-3}-1~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ by way of the steady $\alpha$-disc model with detailed treatment of chemical and thermal processes. We have considered two possible compositions of the inflow, i.e., mostly atomic and molecular gases, at the disc outer boundary. We have constructed marginally gravitationally stable disc structure. By comparing the required value of viscosity parameter $\alpha$ to maintain such structure and the possible maximum value of unity, we have judged the gravitational instability of discs.

In the case of atomic flows, we have found the following two types of disc structures depending on the accretion rate. At high rate ($\gtrsim 10^{-2}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$), the gas composition remains atomic in the optically thin part with temperature of $\sim 3000$ K, set by the H- free-bound emission. At lower rate, in contrast, the compositional transition to the H2 phase occurs at some radius and $\mathrm{H}_{2}$ collision-induced emission dominates the cooling further inward. In both cases, once the disc becomes optically thick, the temperature jumps up to $>10^{4}$ K and the gas is ionized. The discs are gravitationally unstable if the accretion rate is higher than $10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ (Fig. 5) regardless of the central stellar mass.

If the inflow is mainly composed of a molecular gas, the $\mathrm{H}_{2}$-line cooling keeps the temperature below 2000 K in outer part of the disc although the temperature gradually increases inward. With the accretion rate higher than $10^{-2}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$, the $\mathrm{H}_{2}$ dissociation occurs while the disc is still optically thin and the temperature jumps up to 3000K. With a lower rate, the gas remains in the molecular form until it becomes optically thick to the H2 collision-induced absorption. For a given accretion rate and stellar mass, such disc is always unstable at sufficiently large radius (Fig. 10). This critical radius is $10^{2}-10^{3}$ au for $\lesssim 10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ depending on the stellar mass, but decreases abruptly by an order of magnitude with higher accretion rates.

Based on the analysis above, let us briefly discuss the possibility of disc fragmentation during SMS formation in different formation sites, whose main component can be either atomic or molecular. First, for a star forming from the atomic gas, the accretion rate can be estimated as

[TABLE]

where $c_{\mathrm{s,col}}$ is the effective sound speed in the collapsing cloud including also the effect of turbulence, magnetic field etc., and $\phi$ is a factor that reflects the manner of collapse: $\phi\sim 10$ for dynamically collapsing cores and $\phi\sim 1$ for quasi-statically collapsing cores (Hunter, 1977). According to numerical simulations for the collapse induced by atomic cooling (e.g., Inayoshi et al., 2014), $c_{\mathrm{s,col}}$ is $\sim 10~{}\mathrm{km\ s^{-1}}$, which indicates the occurrence of disc fragmentation largely depends on detailed manner of the collapse. Second, for the molecular-gas case, where SMSs are supposed to form in sites with large streaming velocity between baryon and dark matter, whether the size of the discs is bigger or smaller than the critical one (Fig. 10) matters. In Hirano et al. (2017)’s numerical simulations, disc radii are $\sim$1000 au in most of the cases, while the accretion rates fluctuate between $10^{-1}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ and $10^{-2}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$. Thus, the discs are larger than the critical radius and expected to fragment before the central stellar mass reaches $\sim 10^{4}~{}\mathrm{M}_{\odot}$ (Fig. 10), which is actually the case as found by Hirano et al. (2017).

If a disc fragments into multiple clumps, the accretion onto the central protostar becomes episodic with repeating bursts accompanying falls of the clumps and intervening quiescent phases. With a quiescent phase longer than the Kelvin-Helmholtz timescale, the protostar contracts to a main-sequence star and starts emitting ionization photons, which can terminate further growth of the central star (Sakurai et al., 2015, 2016). For quantitative modeling of the accretion-rate variability, hydrodynamical simulations incorporating with the chemical and thermal networks are awaited.

The outcome of the disc fragmentation studied here might be supermassive binary stars, which can be detectable by future gravitational wave observations (Chon et al., 2018). For example, if an equal mass binary system with $10^{4}~{}\mathrm{M}_{\odot}$ forms and ends up with the merger of the remnant BHs with the same mass, the emitted gravitational waves would be detectable by Deci-hertz Interferometer Gravitational wave Observatory (DECIGO: Kawamura et al. 2011) and Laser Interferometer Space Antenna (LISA: Amaro-Seoane et al. 2012) up to redshift $z\sim 30$ and $\sim 10$, respectively.

In this work, we have adopted the critical $\alpha$-value of unity as the condition for disc fragmentation in light of numerical simulations of protoplanetary discs (Zhu et al., 2012). Fragmentation condition is, however, still in dispute (cf., Gammie, 2001; Meru & Bate, 2012; Rice et al., 2014; Takahashi et al., 2016). Furthermore, it is not clear whether the condition obtained for protoplanetary discs is applicable to the those around SMSs. If the critical value of $\alpha$ becomes smaller, discs will fragment more easily. In particular, for atomic flows with low accretion rate ($\lesssim 10^{-2}~{}\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$), the cold inner regions by $\mathrm{H}_{2}$ cooling (see Figs. 3 and 4) may also become gravitationally unstable if a smaller critical $\alpha$-value is adopted.

Here the disc structure is obtained under the thin disc approximation. To check this, we calculated the aspect ratio $H/r$ for all the cases considered in this paper and found that in most cases it falls in the range $\sim 10^{-1}-10^{-2}$, justifying this approximation. Owing to the dependence $H/r\propto r^{1/2}M_{\ast}^{-1/2}$, however, the thin disc approximation becomes worse for lower $M_{\ast}$ and larger $r$. In fact, the ratio $H/r$ becomes $\sim$ 1 in our smallest mass case considered $M_{\ast}=10~{}\mathrm{M}_{\odot}$ around the outer edge of the disc at $\sim$ 1000 au. In reality, however, the disc size would be smaller than our adopted value of 1000 au in such an early stage since the disc size increases with time accompanying the accretion of large angular momentum material originally located at the outer part of the parental dense core.

Considering uncertainties above, further studies are needed to conclude whether SMSs actually form or not. We plan to perform hydrodynamics simulations that follow the formation of SMSs from the collapse of cores, where the disc fragmentation and subsequent evolution of fragments supposedly play an important role.

Acknowledgments

The authors would like to thank Kohei Inayoshi, Daisuke Nakauchi, Hidekazu Tanaka and Hidenobu Yajima for discussions. We also thank Shingo Hirano and Takashi Hosokawa for providing their simulation results and the anonymous reviewer for useful comments on improving thermal and chemical modelling. RM acknowledges financial support from the Graduate Program on Physics for the Universe of Tohoku University. This work is supported in part by MEXT/JSPS KAKENHI Grant Number 17H06360 (KS and KO) and 17H01102, 17H02869 (KO) and by NAOJ ALMA Scientific Research Grant Numbers 2016-02A (ST).

Appendix A Radiative processes

A.1 Continuum cooling rate in the optically thin regime

For continuum emission cooling, we consider the reverse of absorption processes in Table 1. Here, all the emission processes are in the form either of $\mathrm{A+B\to C+\gamma}$ or $\mathrm{A+B\to C+D+\gamma}$, and their emissivities $\eta_{\nu}$ are proportional to the product of $n(\mathrm{A})$ and $n(\mathrm{B})$. More specifically, the emissivity of photons per unit volume via free-bound emission ($\mathrm{A+e\to B+\gamma}$) is given by (Omukai, 2001):

[TABLE]

that via free-free emission ($\mathrm{A}+\mathrm{e}\to\mathrm{A}+\mathrm{e}+\gamma$) or $\mathrm{H}_{2}$ CIE ($\mathrm{H}_{2}(=A)+\mathrm{B}\to\mathrm{H}_{2}+\mathrm{B}+\gamma$)

[TABLE]

By integrating over frequency, the cooling rate by $i$-th process in the optically thin regime is

[TABLE]

with $f_{i}(T)$ a function only of the temperature.

The coefficients $f_{i}(T)$ are calculated by numerically integrating Equations (27) and (28) in the range $10^{-3}\,\mathrm{eV}<h\nu<10^{2}\,\mathrm{eV}$ and the results are fitted for the temperature range $500\,\mathrm{K}<T<20000\,\mathrm{K}$ except in the case of $\mathrm{H}^{-}$ free-free emission ($900\,\mathrm{K}<T<20000\,\mathrm{K}$ in this case) in the form of

[TABLE]

The fitting constants $a_{0}$, $a_{1}$, $a_{2}$, $a_{3}$ and $a_{4}$ are presented in Table 2 for each emission process. The error in the fitting is generally less than a few percent.

A.2 Planck and Rosseland mean opacities

The cooling rate in the optically thick regime is reduced from the optically thin value due to the photon trapping effect, which is taken into account as in Equation (13) by using the Planck and Rosseland optical depths ($\tau_{\rm P}$ and $\tau_{\rm R}$). Here we present the method for calculation of Planck and Rosseland mean opacities. The continuum absorption processes taken into account are listed in Table 1. In addition to those considered in Omukai (2001), we also consider the $\mathrm{H}_{2}$ Rayleigh scattering following Mayer & Duschl 2005. The opacity $\kappa_{\nu}$ per unit mass of gas via a pure absorption process, i.e., bound-free absorption ($\mathrm{A+\gamma\to B+\mathrm{e}}$), free-free absorption ($\mathrm{A+\mathrm{e}+\gamma\to A+\mathrm{e}}$) and $\mathrm{H}_{2}$ collision-induced absorption (CIA; $\mathrm{H}_{2}(={\rm A})+\mathrm{B}+\gamma\to\mathrm{H}_{2}+\mathrm{B}$), is given by considering stimulated emission:

[TABLE]

with the Planck constant $h$ and the photon frequency $\nu$. That for scattering ($\mathrm{A+\gamma\to A+\gamma^{\prime}}$) is

[TABLE]

Note that the cross-sections $\sigma_{\nu}$ for free-free absorption and $\mathrm{H}_{2}$ CIA are proportional to the density of collision partners, $n(\mathrm{e})$ and $n(\mathrm{B})$, respectively.

The Planck mean opacity can be calculated by

[TABLE]

where subscript $i$ represents an individual process in Table 1 and the sum is over only pure absorption processes, while the Rosseland mean opacity is

[TABLE]

where the sum is over all the absorption and scattering processes.

The photon trapping effect becomes remarkable only when the disc is nearly optically thick. In such condition, the chemical fractions are approximately in the equilibrium values, which are given by Saha-Boltzmann equations. For example, among species A, B, and C for which a reaction $\mathrm{A\leftrightarrow B+C}$ exists, the equilibrium fractions are written as

[TABLE]

where the superscript $*$ means the chemical equilibrium value, $\Delta\chi=\chi_{\mathrm{B}}+\chi_{\mathrm{C}}-\chi_{\mathrm{A}}$ is the difference in chemical binding energy, and $m_{\mathrm{X}}$, $z_{\mathrm{X}}$ and $\chi_{\mathrm{X}}$ are the mass, partition function and chemical potential energy of species $\mathrm{X}$, respectively. We adopt $z_{\mathrm{H}}=z_{\mathrm{e}}=2$, $z_{\mathrm{H}^{+}}=z_{\mathrm{H}^{-}}=1$, $\chi_{\mathrm{H}}=13.6\ \mathrm{eV}$, $\chi_{\mathrm{H}_{2}}=4.48\ \mathrm{eV}$, $\chi_{\mathrm{H}^{-}}=0.755\ \mathrm{eV}$, and the fitting function of $z_{\mathrm{H}_{2}}$ from Irwin (1981).

In Tables 3 and 4, we present the numerical values for the Planck and Rosseland mean opacities with the equilibrium chemical fractions for the temperature and density ranges $2.70<\mathrm{log}\ T\ (\rm K)<4.30$ and $-17.0<\mathrm{log}\ \rho\ (\mathrm{g}\ \mathrm{cm}^{-3})<-6.0$. The frequency integrations in Equations (33) and (34) are carried out in the range of $10^{-2}<h\nu/k_{\mathrm{B}}T<10^{2}$.

Appendix B Lithium cooling rate

We calculate the cooling rate of atomic lithium bound-bound emission by counting level populations from the statistical balance among ten (from 2s up to 5s) levels. Photon trapping effect in large column density cases is taken into account by way of the escape probability formalism (e.g. Omukai 2001). The level energies and the spontaneous radiative decay rate are taken from Klevas et al. (2016). We only consider electron impact as the collisional excitation process and its rate coefficients are taken from Osorio et al. (2011). The cooling rate including the line self-absorption $\overline{\beta}_{\mathrm{esc}}\Lambda_{\mathrm{Li,thin}}=n(\mathrm{Li})n(\mathrm{e})L_{\mathrm{Li}}/\rho~{}\mathrm{(erg~{}g^{-1}~{}s^{-1})}$ is calculated for given sets of the temperature $T$, the electron number density $n({\rm e})$ and the lithium column density $N({\rm Li})$. We found the results are excellently reproduced by the interpolation in the form of

[TABLE]

where $L_{0}$ is the cooling function in the low-density limit and $\overline{\beta}_{\mathrm{esc}}\mathcal{L}_{\mathrm{LTE}}$ is the cooling rate per lithium atom in the LTE case with the photon trapping effect included, which are presented in Table 5 as a function of the temperature $T$ and the column density parameter $\tilde{N}(\mathrm{Li})=N(\mathrm{Li})/\sqrt{2k_{\mathrm{B}}T/m_{\mathrm{Li}}}$, where $N(\mathrm{Li})$ is the column density of lithium and $m_{\mathrm{Li}}$ is the mass of a lithium atom. The error in this approximation is generally less than one percent and at most several percent in the presented range even in the worst case at high temperatures and the electron number densities.

Appendix C Chemical Reactions

Our chemical network consists of reactions 1-10 in Table 6, including both forward and reverse reactions. The forward (i.e., left to right) reaction coefficients are presented in the table. The reverse reaction coefficient $k^{\prime}_{\rm c}$ for collisional reactions (reactions 1-8 in Table 6), ${\rm A}+{\rm B}\rightleftharpoons{\rm C}+{\rm D}$, whose forward reaction coefficient is $k_{\rm c}$, is given from the law of mass action as

[TABLE]

with equilibrium constant $K(T)=\left[\frac{n({\rm C})n({\rm D})}{n({\rm A})n({\rm B})}\right]^{*}$ calculated by the Saha-Boltzmann equation. The reverse reaction coefficients of radiative reactions (i.e., reactions 9 and 10), which are expressed symbolically as ${\rm A+B\rightleftharpoons C+\gamma}$ depend on the radiation field. The radiation intensity $J_{\nu}$ in the disc is proportional to the effective optical depth for the contiuum $\tau=\sqrt{\tau_{\mathrm{P}}\tau_{\mathrm{R}}}$ and saturates to the black-body value $B_{\nu}$ at $\tau\simeq 1$. Recalling that the chemical equilibrium is realized for $J_{\nu}=B_{\nu}$, the reverse reaction coefficient $k^{\prime}_{\rm r}$ is related to the forward reaction coefficient $k_{\rm r}$ as

[TABLE]

where the equilibrium coefficient $K(T)=\left[\frac{n({\rm C})}{n({\rm A})n({\rm B})}\right]^{*}$.