Formation of"Blanets"from Dust Grains around the Supermassive Black Holes in Galaxies

In Wada, Tsukamoto, and Kokubo (2019), we proposed for the first time that a new class of planets,"blanets"(i.e., black hole planets), can be formed around supermassive black holes (SMBHs) in the galactic center. Here, we investigate the dust coagulation processes and physical conditions of the blanet formation outside the snowline ($r_{snow} \sim $ several parsecs) in more detail, especially considering the effect of the radial advection of the dust aggregates. We found that the viscous $\alpha$-parameter in the turbulent circumnuclear disk should be smaller than 0.04, to prevent the destruction of the aggregates due to collisions. The formation timescale of blanets $\tau_{GI}$ at $r_{snow}$ is, $\tau_{GI} \simeq$ 70-80 Myr for $\alpha = 0.01-0.04$ and $M_{BH} = 10^6 M_\odot$. The mass of the blanets ranges from $\sim 20 M_E$ to $3000 M_E$ in $r<4$ pc for $\alpha = 0.02$ ($M_E$ is the Earth mass), which is in contrast with $ 4 M_E-6 M_E$ for the case without the radial advection. Our results suggest that \textit{blanets} could be formed around relatively low-luminosity AGNs ($L_{bol} \sim 10^{42}$ erg s$^{-1}$) during their lifetime ($\lesssim 10^8$ yr).


INTRODUCTION
There is enough evidence suggesting that planets are formed in the protoplanetary disks around stars. However, protoplanetary disks might not be the only site for planet formation. Recently, in Wada, Tsukamoto, and Kokubo (2019) (hereafter Paper I), we claimed a new class of "planets" that orbit around super-massive black holes (SMBHs) in galactic centers. Paper I theoretically investigated the growth processes of planets, from sub-micron-sized icy dust monomers to Earth-sized bodies outside the snowline in a circumnuclear disk around a SMBH, typically located several parsecs from the SMBHs. As is the case in a protostellar disk, in the early phase of the dust evolution, low-velocity collisions between dust particles promote sticking; therefore, the internal density of the dust aggregates decreases with growth (Okuzumi et al. 2012;Kataoka et al. 2013). When the size of porous dust aggregates reaches 0.1-1 cm, the collisional and the gas-drag compression become effective, and as a result, the internal density stops decreasing. Once 10-100 m sized aggregates are formed, they decouple from gas turbulence, and as a result, the aggregate layer becomes gravitationally unstable (Michikoshi, & Kokubo 2016, leading to the formation of "planets" due to the fragmentation of the layer, with ten times the mass of the earth. The objects orbit the SMBHs with an orbital time of 10 5 − 10 6 years. To distinguish them from standard planets, we hereafter call these hypothetical astronomical objects blanets (i.e., black hole planets).
The results reported in Paper I, however, have two major limitations. One is that the collisional velocity between the dust aggregates might become too large (> several 100 m s −1 at the Stokes parameter, S t ∼ 1). And if the collisional velocity is that large, rather than growing, the aggregates might get destroyed. In Paper I, we used the numerical experiments conducted by Wada et al. (2009) 1 on the collisions between the dust aggregates, wherein the critical collisional velocity (v crit ) scales with the mass m d of the dust aggregates, as v crit ∝ m 1/4 d . However, this is correct only for the headon collisions, as stated in the paper. Moreover, Wada et al. (2009Wada et al. ( , 2013 showed that the growth efficiency of the dust aggregates depends on the impact parameter of the collisions, and as a result, v crit does not strongly depend on the mass of the dust aggregates, if off-set collisions are taken into account. They concluded that v crit 80 m s −1 for the icy monomers. This low critical velocity is also one of the obstacles in the planet formation in protoplanetary disks. In this follow-up paper, we adopt v crit 80 m s −1 as a constraint on the growth of the dust aggregates. Another limitation of Paper I is that the size of dust aggregates a d and collisional velocity ∆v show runaway growth in the collisional compression phase around S t ∼ 1. However, this rapid growth would not be realistic if a more natural treatment of the internal density of the dust is considered ( §2.2.1, see also §3).
Moreover, there is a critical process that may promote blanet formation. In paper I, we did not take into account the radial drift of the dust particles as the first approximation. The radial velocity of the dust v r,d relative to the gas (Weidenschilling 1977;Tsukamoto et al. 2017 where c s is the gas sound velocity and v K is the Keplerian rotational velocity. In the circumnuclear disk around a SMBH, initially S t η ∼ 10 −4 − 10 −3 . Then the drift time of the dust particle t drif t ∼ r/v r,d ∼ 5 − 50(M BH /10 7 M ) 1/2 (r/1 pc) 1/2 Myr. This is not negligibly small for the lifetime of the active galactic nucleus (AGN), i.e., 10 7 − 10 8 yr. In this paper, we investigate the effects of the radial advection of the dust particles.
The remainder of this paper is organized as follows. In §2, we describe the models for the dust evolution and its application to the circumnuclear region. In §3, we show the results of the models with and without the radial advection of the dust particles. In §4, we discuss how the maximum collisional velocity and the formation timescale of blanets depend on the parameters α and M BH . We also discuss the expected mass of the blanets and their radial distribution. Finally, we summarize the results in §5. Here we briefly summarize the concept of dust evolution around SMBHs, as discussed in Paper I. Figure 1 shows a schematic of the active galactic nucleus (AGN) and the circumnuclear disk. A SMBH (with a mass of 10 6 − 10 10 M ) is surrounded by an accretion disk, which radiates enormous energy (the bolometric luminosity is ∼ 10 42 − 10 45 erg s −1 ), mostly as ultra-violet and X-rays. The dust particles in the central r < r sub are sublimated by the nuclear radiation. The radius depends on the AGN luminosity: (1) where L U V is the ultra-violet luminosity of the AGN, and a d is the dust size (Barvainis 1987). The radiation forms conical ionized gas (narrow emission-line region) and also contributes to producing the outflows of the dusty gas and torus (Wada 2012;Wada et al. 2018;Izumi et al. 2018). In the mid-plane of the torus, cold, dense gas forms a thin disk, where icy dust particles can be present beyond the snowline r snow (see §2.4).
We model the turbulent disk based on the α-viscosity formalism (Shakura, & Sunyaev 1973) with a dimensionless parameter α 1.0. In the circumnuclear disk in AGNs, the value of α is highly uncertain. Therefore, we here treat α as a free parameter to check how it alters the results, especially the maximum collisional velocity between the dust aggregates and the onset of the gravitational instability of the dust disk.
In contrast to the dust coagulation process in the protoplanetary disks (Weidenschilling 1977), the drag between dust particles and gas obeys the Epstein law. The aggregate's size (a d ) is always much smaller than the mean free path of the gas, λ g ∼ 10 12 cm(σ mol /10 −15 cm 2 ) −1 (n mol /10 3 cm −3 ) −1 , where σ mol and n mol are the collisional cross-section and number density of the gas, respectively.

Evolution of dust aggregates in each stage
The model for the growth of dust particles here is based on the elementary processes found around stars. The evolution of dust particles is divided into four stages as described below.

Hit-and-Stick stage
If the dust aggregates grow through ballistic cluster-cluster aggregation (BCCA), the internal structure of the aggregate should be porous (i.e., the internal density is much smaller than the monomer's density: ρ int ρ 0 ), and its fractal dimension is D 1.9 (Mukai et al. 1992;Okuzumi et al. 2009). This is called the hit-and-stick stage (Okuzumi et al. 2012;Kataoka et al. 2013), and the internal density is given by where m d is the mass of the aggregate, and m 0 is the monomer's mass. We assume that m 0 = 10 −15 g and ρ 0 = 1 g cm −3 . The growth rate for m d is where ∆v is collision velocity between the aggregates, and H d is the scale height of the dust disk as given in (Youdin, & Lithwick 2007;Tsukamoto et al. 2017); where H g = c s /Ω K is the scale height of the gas disk, and the Stokes parameter is defined as The collision velocity between aggregates ∆v for S t < 1 can be divided into two regimes (Ormel, & Cuzzi 2007): regime I) t s t η = t L Re −1/2 , and regime II) Here the Reynolds number, R e ≡ αc 2 s /(ν mol Ω K ) with the molecular viscosity ν mol ∼ 1 2 c s λ g is where Q g is the Toomre's Q-value for the gas disk and γ Edd is the Eddington ratio for the AGN. The eddy turn over time t L is t L ∼ Ω −1 K , and t η ∼ t L for the smallest eddy. For the hit-and-stick stage, S t R −1/2 e , then for the regime I, where S t,1 and S t,2 are Stokes numbers of the two colliding particles, and C I is a constant of the order of unity (Sato et al. 2016). For regime II, on the other hand, where v L is velocity of the largest eddy. We assume that ∆v I = ∆v II at the transition.
The size of dust aggregates determines how they interact with the gas. The dynamics of the aggregates is affected by their cross sections, which depend on their internal inhomogeneous structure. The radius of BCCA cluster a BCCA consisted of N monomers (N = m d /m 0 ) is given as a BCCA N 0.5 a 0 for N 1 (Mukai et al. 1992;Wada et al. 2008Wada et al. , 2009, which was also confirmed through N -body simulations (Suyama et al. 2012). We then assume that

Compression stages
The hit-and-stick stage ends due to collisions between the aggregates (collisional compression), or due to their interaction with the ambient gas (gas-drag compression). In the collisional compression, the rolling energy E roll , which is the energy required to rotate a particle around a connecting point by 90 • , is comparable to the impact energy, E imp = m d ∆v 2 /4, between the two porous dust aggregates of the same mass, m d . Beyond this point, the aggregates start to get compressed due to mutual collisions and interaction with the gas (i.e., the ram pressure).
According to Suyama et al. (2012), the internal density of the aggregated ρ int,f formed by collisions between two equal-mass aggregates, with density ρ int , is calculated for E imp > 0.45E roll : ρ f is the fractal density of the dust aggregate: ρ f ≈ m d /(7.7a 2.5 d ), and E roll = 4.37 × 10 −9 erg.
Moreover, the fluffy dust aggregates can be compressed owing to the ram pressure of the ambient gas (Kataoka et al. 2013). The internal density of the aggregates that are compressed by the gas is given where the ram pressure for a dust aggregate is with the stopping time t s = S t /Ω K (Kataoka et al. 2013).

N -body stage
When S t 1, the aggregates are decoupled from the turbulent gas, and they evolve as N -body system, after which the collision velocity between the aggregates is determined by a balance between heating and cooling processes as the N -body particles. According to Michikoshi, & Kokubo (2016, we solve the following equation to get the equilibrium random velocity of the dust aggregates v d , The first three heating terms represent the gravitational scattering of the aggregates, stirring by the turbulence, and gravitational scattering by density fluctuation of the turbulence, respectively. The two cooling terms in eq.(13) represent the collisional damping and the gas drag. We assume the collision velocity ∆v ≈ v d .

Radial advection of the dust particles
In Paper I, we ignored the radial advection of the dust particles in the disk. However, as mentioned in §1, this is not always obvious. Here, we solve the following governing equations for the dust evolution based on the assumption that the mass distribution of the dust particles at each orbit radius is singly peaked at a mass (Tsukamoto et al. 2017;Sato et al. 2016); Here t coll in eq.(15) is the collision time, and the source term is where n d is the number density of the dust particles. The dust particles have a radial velocity due to the drag with the ambient gas: where v K is the Kepler velocity and η is a parameter that determines the sub-Keplerian motion of the gas, and the radial velocity of the gas v r,g is given with the mass accretion rateṀ : where the mass accretion rateṀ is assumed to be using the Eddington mass accretion rateṀ = γ EddṀEdd with the Eddington ratio γ Edd .

Gravitational instability of the dust disk and formation of blanets
We investigate the gravitational instability (GI) of the disk consisting of dust aggregates with S t > 1 using the Toomre's Q-value defined as For the axi-symmetric mode, Q d < 1 is the necessary condition for GI, but the non-axisymmetric mode can develop for Q d 2. In this case, spiral-like density enhancements are formed followed by fragmentation of the spirals (Michikoshi, & Kokubo 2017), which leads formation of massive objects, i.e., blanets. The mass of blanets is estimated as where the critical wavelength for GI is

Initial and boundary conditions
In all the models, the circumnuclear cold gas disk embedded in the geometrically thick torus (see Fig. 1) is assumed to be gravitationally stable; the Toomre's Q-value, Q g ≡ c s Ω K /πGΣ g = 2. The gas sound velocity is assumed to be c 2 s = k B T g /µ m H with T g = 100 K and µ = 2.0.
The Eddington ratio is assumed to be γ Edd = 0.01. The AGN bolometric luminosity is then L bol = 1.3 × 10 42 erg s −1 M BH /(10 6 M ). The X-ray luminosity of the AGN, which is used to determine the snowline radius, L X = 0.1L bol . This can be attributed to the fact that the UV flux from the accretion disk is attenuated in the dense circumnuclear disk.
The snowline for a d = 0.1 µm, r snow ≈ 1.5 pc L X 1.3 × 10 41 erg s −1 Therefore, it is expected that the dust in the most part of the circumnuclear disk is icy. We assume T ice = 170 K. We solve the governing equations ( §2.2.4) between r = 0.1 pc and 200 pc with 600 grid cells.  (Wada et al. 2009). After S t = 1 is attained, ∆v drops and the disk of the aggregates becomes gravitationally unstable. Figure 2a shows a typical evolution of a dust aggregate at the snowline for M BH = 10 6 M and α = 0.02. The internal density of the aggregate ρ int is plotted as a function of its mass m d . Initially, the internal density decreases monotonically from the monomer's initial density, i.e., ρ 0 = 1 g cm −3 to 4 × 10 −6 g cm −3 , as its mass increases from m d ∼ 10 −15 g to ∼ 10 −5 g. At that instant, the size of the aggregate becomes ∼ 1 cm (see Fig. 2b). After this hit-and-stick phase, the fluffy dust aggregates keep growing due to collisions in the turbulent gas motion until S t 1. During this stage (m d = 10 −5 g to 10 10 g), the aggregates are compressed mainly due to the gas drag ( §2.2.2), and therefore ρ int gradually increases 2 . After S t = 1, the aggregates are decoupled from the turbulence and the system goes to the N -body evolution phase. For m d > 10 10 g and S t > 1, the aggregates are compressed due to their self-gravity. In the model shown in Fig. 2, it becomes gravitationally unstable at t = 75 Myr (see also §4.1). Figure 2b plots the collisional velocity ∆v of the aggregate and its size a d as a function of m d . The size a d monotonically increases. In the compression stage (m d > 10 5 g), the increase of a d slows down (see eq. (9)). Initially ∆v is 20 cm s −1 , and it slightly decreases during the hit-and-stick stage. Then it turns to an increase phase until ∼ 57 m s −1 around S t = 1 during the compression stage. The size of the aggregate becomes a d ∼ 10 4 cm at the end of this stage. In this case, the aggregates are not compressed by their self-gravity (m d < 10 10 g) before the dust disk becomes GI. Note that the increase of ∆v slows down at m d ∼ 0.1 g, which corresponds to the transition between ∆v I and ∆v II (eqs. (7) and (8)).  The growth time of an aggregate for ∆v = ∆v I 1/2 √ αc s R 1/4 e S t can be estimated as (23) Figure 3 shows time evolution of ρ int , a d , S t , and ∆v for the same model shown in Fig. 2. The hit-and-stick phase lasts for ∼ 10 Myr as expected by t grow , and S t becomes unity at t = 60 Myr. At this moment, the dust aggregate's size reaches ∼ 100 m (Fig. 3b). Fig. 3b shows that the growth of a d is exponential, or slower in time, in contrast to the results in Paper I. This is a natural consequence of the evolution of the mass of the dust aggregates. The mass increase rate of the dust is

RESUTS
Here we assume the dust layer is sedimanted, i.e., its thickness H d is scaled as H d ∝ S −1/2 t . For ∆v ∝ S 1/2 t (eq. (8)), dm d /dt ∝ m d , therefore, m d grows exponentially. The runaway growth seen in Paper I is caused by the assumption that the internal density of the dust aggregates stays porous (i.e., the fractal dimension is ∼ 2) through the evolution. In reality, when the compression by the ambient gas works, ρ int is nearly constant (i.e., m d ∝ a 3 d ), as shown in Fig. 2, therefore S t ∝ a d . If the scale height of the dust disk is constant, then dm d /dt ∝ S 1/2 t a 2 d ∝ m 5/6 d ; therefore, the growth of the dust aggregate should be slower than exp(t). Fig. 3c shows that the collisional velocity ∆v gradually decreases during the hitand-stick stage, and it turns to rapid increase during the compression stage until S t becomes unity at t = 56 Myr. In the N -body stage, the collisional velocity decreases from its maximum value, 57 m s −1 , and it becomes GI (i.e., Q d ≤ 2) at t 75 Myr.
For comparison, the evolution of the model without the radial advection is shown in Fig. 4. We found that the dust aggregates before S t = 1 evolve almost the same way as that in the model with the radial advection (Fig. 3). However, the time for GI is 136 Myr, in contrast to 71 Myr for the case with the advection.  Fig. 2. (a) Surface density of dust Σ d and the scale height of the dust H d . The gray line is Σ d at t = 0. Note that the surface density of the dust decreases in the outer disk (r > 30 pc), and the total mass of the dust is conserved. The vertical dashed line is the position of the snowline r snow = 1.5 pc. (b) Same as (a), but for the collision velocity ∆v, the internal density of the aggregate ρ int and the Stokes parameter S t . (c) Same as (a), but for the radial velocity of the dust v r,d and v r,g normalized by the Kepler velocity v K . (d). Same as (a), but for the mass and size of the aggregate, m d and a d . Figure 5 shows the radial distributions of Σ d , H d , ∆v, the radial velocity of the dust and gas (v r,d and v r,g ), ρ int , S t , m d , and a d at 200 Myr in the same model shown in Fig.  2. As Fig. 5a shows, the dust is accumulated around r ∼ 2 − 3 pc, where v r,d v r,g (Fig. 5c) and S t turns from S t < 1 to S t > 1 (Fig. 5b). From Fig. 5d, the dust aggregates evolve more rapidly in the inner region (r 3 pc), and the maximum size is ∼ km. We call these objects as blanetesimals. In the models with the radial advection of the dust aggregates, we investigated how the maximum velocity of the collision ∆v max and the time for GI (τ GI ) depend on α and M BH . In Figure 6, we plot ∆v max and τ GI as a function of α for M BH = 10 6 M and 10 7 M . It shows that ∆v max depends on α, and not on M BH . If ∆v max 80 m s −1 is necessary for collisional growth as numerical experiments suggested (Wada et al. 2009), then α should be ∼ 0.04 or smaller. The behavior of the dust growth (e.g., Fig. 2) does not significantly depend on α and M BH , but the timescale to reach S t = 1 is different as shown in Fig. 6b. For example, for M BH = 10 6 M and α = 0.02, it takes ∼ 60 Myr when S t exceeds unity, whereas it is ∼ 225 Myr for M BH = 10 7 M and α = 0.02. This implies that smaller BHs may preferentially host blanets within a lifetime AGNs ( 10 8 yr). Figure 6b shows that, for M BH = 10 7 M , the blanetesimal disk may not become GI earlier than ∼ 150 yr for α < 0.04. For α > 0.05 or 0.06, GI does not occur at r = r snow in the models with M BH = 10 6 M or 10 7 M , respectively.

Number and mass of blanets
In the final stage of the evolution, the blanetecimal disk can be gravitationally unstable, and it fragments into massive objects, i.e., blanets (see §2.3). Figure 7 shows the radial distribution of the mass and typical separation between blanets, λ bl ≈ λ GI (eq.(21)). Two models with the radial advection for M BH = 10 6 M with α = 0.02 and M BH = 10 7 M with α = 0.04 are shown. For comparison, a model without the radial advection is also shown (M BH = 10 6 M and α = 0.02). The mass of blanets ranges from 20M E at r = r snow to 3000M E at r ∼ 3.5 pc for M BH = 10 6 M , in contrast to the model without advection, which is M bl 3M E − 7M E . For M BH = 10 7 M , M bl 10 4 M E − 10 5 M E outside the snowline. However, this extraordinary massive blanet is unlikely, because it is comparable to the minimum mass of brown dwarfs (∼ 2 × 10 4 M E ).
Therefore, the largest size of the blanets would be maximally ∼ 10× Earth's radius at r ∼ 3 pc for M BH = 10 6 M , if the average internal density is similar to that of the Earth.
It would be interesting to investigate whether blanets can obtain the atmosphere. The gas envelope mass attained by the blanet, M E , can be simply estimated by using the Hill radius r H = (M bl /3M BH ) 1/3 r, provided that the gap caused by the accretion is not refilled, If this is the case, the blanets probably cannot obtain the massive envelope of the gas.

SUMMARY
In this follow-up paper of Wada et al. (2019) (Paper I), we theoretically investigated a process of dust evolution around a SMBH in the galactic center. We proposed that a new class of planets, blanets (i.e., black hole planets) can be formed, provided that the standard scenario of planet formation is present in the circumnuclear disk. Here, we investigated the physical conditions of the blanet formation outside the snowline (r snow ∼ several parsecs) in more detail, especially considering the effect of the radial advection of the dust aggregates. We also improved the dust evolution model in Paper I in terms of the internal density evolution of the dust aggregates. We assumed the maximum collisional velocity for destruction, which was suggested by previous numerical simulations, as one of necessary conditions for blanet formation. We found that the viscous parameter α (Shakura, & Sunyaev 1973) in the turbulent circumnuclear disk should be smaller than 0.04 for the black hole mass M BH = 10 6 M ; otherwise, the dust aggregates could be destroyed due to collisions. The formation timescale of blanets τ GI at r snow is τ GI 70-80 Myr for α = 0.01 − 0.04. The blanets (M bl ) are more massive for larger radii; they range from M bl ∼ 20M E − 3000M E in r < 4 pc, in contrast to M bl = 3 − 7M E for the case without the radial advection. The typical separations between the blanets, estimated from the wavelength of the gravitational instability, would be ∼ 0.01 pc. For M BH ≥ 10 7 M , the formation timescale is longer than ∼ 150 Myr for α ≤ 0.04. Although the gravitational instability of the blanetesimal disk takes place just outside the snowline (r = 4.7 pc), they should not be blanets because they are more massive than the typical brown dwarf mass (∼ 3 × 10 4 M E ). Note that AGNs are often heavily obscured with dense gas even for hard X-rays (Buchner et al. 2014) (N H > 10 23 cm −2 , i.e., Compton-thick). If this is the case, the snowline is located at the inner region (e.g. r ∼ 2 − 3 pc), and as a result, blanets with M bl ∼ 10M E − 100M E around M BH = 10 7 M could be possible.
Our results suggest that blanets could be formed around relatively low-luminosity AGNs during their lifetime ( 10 8 yr). The gaseous envelope of a blanet should be negligibly small compared with the blanet mass. Therefore, the system of blanets are extraordinarily different from the standard Earth-type planets in the exoplanet systems. The dynamical stability of such a system around a SMBH may be an interesting subject for future studies.