Planet Formation around Supermassive Black Holes in the Active Galactic Nuclei

As a natural consequence of the elementary processes of dust growth, we discovered that a new class of planets can be formed around supermassive black holes (SMBHs). We investigated a growth path from submicron sized icy dust monomers to Earth-sized bodies outside the “snow line,” located several parsecs from SMBHs in low luminosity active galactic nuclei (AGNs). In contrast to protoplanetary disks, the “radial drift barrier” does not prevent the formation of planetesimals. In the early phase of the evolution, low collision velocity between dust particles promotes sticking; therefore, the internal density of the dust aggregates decreases with growth. When the porous aggregate’s size reaches 0.1–1 cm, the collisional compression becomes effective, and the decrease in internal density stops. Once 10–100 m sized aggregates are formed, they are decoupled from gas turbulence, and the aggregate layer becomes gravitationally unstable, leading to the formation of planets by the fragmentation of the layer, with 10 times the mass of the Earth. The growth timescale depends on the turbulent strength of the circumnuclear disk and the black hole mass MBH, and it is comparable to the AGN’s lifetime (∼108 yr) for low mass (MBH ∼ 106M⊙) SMBHs.


Introduction
Planetary systems are ubiquitous-more than 4000 exoplanets have been discovered thus far. 5 However, protoplanetary disks around stars may not be the only sites for planet formation in the universe. Here we propose a new site of "planet" formation: the circumnuclear disk around supermassive black holes (SMBHs).
Most galaxies host SMBHs at their centers, with masses ranging from a few million to billion solar masses. Gas disks around SMBHs emit large amounts of energy owing to mass accretion onto the SMBHs, which are known as the "central engine" of active galactic nuclei (AGNs). It is believed that the mass of the SMBH in a galaxy depends on its host galaxy's bulge mass (Marconi & Hunt 2003). Researchers are more convinced of the presence of SMBHs since the discovery of the "black hole shadow" in M87 (Event Horizon Telescope Collaboration et al. 2019). In the "unified model" of AGNs (Antonucci 1993;Netzer 2015), the gas and dust form a geometrically and optically thick "torus," and it obscures the broad emission line (line width is several 1000 km s −1 ) region around the central accretion disk. This hypothesis successfully explains the type-1 and type-2 dichotomy of Seyfert galaxies' spectra, depending on the viewing angle of the tori. The real structure of the tori has been unclear for many years. Recently, the Atacama Large Millimeter/submillimeter Array (ALMA) spatially resolved the molecular tori in nearby AGNs (García-Burillo et al. 2016;Imanishi et al. 2018;Izumi et al. 2018;Combes et al. 2019). Their internal structure is still not well resolved; however, recent 3D radiation-hydrodynamic simulations suggested a dynamic structure energized by a radiation-driven fountain flow to sustain their geometrical thickness (Wada 2012;Wada 2015;Wada et al. 2018, see also Figure 1). Notably, even in this situation, cold, dense gas forms a geometrically thin disk (Schartmann et al. 2014;Wada et al. 2016), and this stratified structure is also consistent with recent X-ray surveys (Buchner et al. 2014).
The remainder of this paper is organized as follows. In Section 2, we describe dust and its environment around SMBHs, and their differences from the standard situation, i.e., in the circumstellar disks. In Section 3, we show a typical evolutional track of representative dust particles from a monomer to a planet-sized body. Four stages of the dust coagulation based on the recent theoretical model proposed for the protoplanet disks are described in detail in the Appendix. We also discuss how the evolutional timescales depend on parameters in Section 4.

The Snow Line in AGNs and the Major difference from the Protoplanetary Disks
The dusty gas around SMBHs extends beyond the dust sublimation radius r sub , where the dust temperature is higher than the sublimation temperature of the dust grains (T sub ∼1500 K). The radius depends on the AGN luminosity: where L UV is the ultraviolet luminosity of the AGN and a d is the dust size (Barvainis 1987). The temperature of the gas and dust beyond r sub , especially at the midplane of the dusty torus, should be cold 100 K (Schartmann et al. 2014), because the radiation originated from the accretion disk is weaker in the direction of the disk plane, and it is further attenuated by the dense dusty gas. Interestingly, even for X-rays, a large fraction of AGNs are Compton-thick, i.e., the hydrogen column density is N H > 10 24 cm −2 (Buchner et al. 2014). Although near-infrared and mid-infrared interferometer observations of AGNs show the presence of hot dusts (several 100 K) around AGNs , colder dust particles are also present in this dense media around AGNs. The total amount of dust in the central 6-27 pc around SMBHs estimated from recent molecular lines (e.g., CO) observations of nearby AGNs by ALMA (Combes et al. 2019) is enormous, e.g., ∼(0.7-3.9)×10 5 M e for the dust-to-gas mass ratio of ∼0.01 (Draine 2011). This number could be even larger for the high metallicity environment around AGNs (Groves et al. 2006;Rémy-Ruyer et al. 2014). The internal dynamics and structure of the molecular tori are still observationally unclear. However, since the mass feeding to the AGN through the circumnuclear disk is necessary during their lifetime (∼10 7 -10 8 yr), the turbulent viscosity works in the dusty gas disk. Here, we model the turbulent disk based on the α-viscosity formalism (Shakura & Sunyaev 1973;Shlosman & Begelman 1987).
The dust grains in the circumnuclear disks around SMBHs are in a qualitatively similar situation as the ones in protoplanetary disks. The major differences between the two systems are summarized in Table 1. The "snow line" for a d =0.1 μm dust irradiated by X-ray around an AGN with an SMBH of 10 7 M e for the Eddington ratio (γ Edd ) of 0.1 is  (Barvainis 1987). 6 Moreover, AGNs are often heavily obscured (Compton-thick) even for hard-X-rays (Buchner et al. 2014), suggesting that a cold dusty layer exits around SMBHs. Therefore, it is expected that the dust in most of the circumnuclear disk is icy.

Outline of Evolution of "fluffy" Dust Aggregates
We then apply recent models of coagulation of dust particles and their aggregates outside the snow line in the protoplanetary disk (Okuzumi et al. 2012;Suyama et al. 2012;Kataoka et al. 2013;Michikoshi & Kokubo 2016 to the dust around AGNs. The coagulation of "fluffy icy dust" is one of the plausible solutions to avoid the theoretical obstacles that prevent the growth of dust grains (monomers) to "planetesimal" such as the "radial drift barrier" (Okuzumi et al. 2012). In this scenario, the evolution of the dust can be divided into four stages: (1) the "hit-and-stick" phase, (2) the collisional or gas pressure compression phase, (3) the gravitationally compression phase, and (4) the gravitational instability (GI) phase (e.g., Goldreich & Ward 1973). We investigated each phase in the circumnuclear disk as discussed below (see also the Appendix).
We track the evolution of icy monomers, whose size and density are a 0 =0.1 μm and ρ 0 =1 g cm −1 , and their aggregates. Here, we investigate the evolution of a representative dust particle size, using the single-size approximation (Sato et al. 2016). In the hit-and-stick phase, when two monomers/ aggregates collide, the internal structure of the aggregates becomes porous (i.e., average internal density is smaller than ρ 0 ) with internal voids (Suyama et al. 2012). This "fluffy dust" formation is also examined by numerical experiments (Dominik & Tielens 1997;Wada et al. 2008). The internal density and size of the aggregates are ρ int ∼(m d /m 0 ) −1/2 ρ 0 and a d ∼(m d /m 0 ) 1/2 a 0 , where m d and m 0 are the masses of the aggregate and monomers, for the fractal dimension of 2. When the collision energy exceeds a critical value, the porous Figure 1. Schematic picture of the active galactic nucleus (AGN) and the circumnuclear disk. A supermassive black hole (the mass is 10 6 -10 9 M e ) is surrounded by an accretion disk, which radiates enormous energy (∼10 42 -10 45 erg s −1 ) mostly in the ultraviolet and X-ray. The dust particles in the central r<r sub ∼0.1-a few parsecs are sublimated owing to the heating by the central radiation. The radiation forms conical ionized gas (Narrow emission-line region) and also contributes to producing outflows of the dusty gas and torus (Wada 2012;Izumi et al. 2018;Wada et al. 2018). In the midplane of the torus, cold, dense gas forms a thin disk, where icy dust particles can present beyond the snow line r snow . The dust aggregates evolve by collisions to form planetesimals, and eventually "circum-black-hole planets" by the gravitational instability of the aggregate disk.  6 The approximate proportionalitya d 1 2 comes from the absorption efficiency of a dust grain being roughly proportional to its radius at a certain wavelength in the near-IR (Draine & Lee 1984). aggregates start to get compressed, and the evolution of the internal density changes beyond this point (Suyama et al. 2012). During this compression phase, the aggregates' mass rapidly increase; however, their internal densities gradually increase as well, from ρ int ∼10 −6 g cm −3 to ∼10 −5 g cm −3 .
In contrast to the dust coagulation process in the protoplanetary disks (Weidenschilling 1977), the drag between dust particles and gas obeys the Epstein law only. The aggregate's size (a d ) is always much smaller than the mean free path of the gas ( , where σ mol and n mol are the collisional cross-section and number density of the gas). At all times the radial drift velocity of the dust is negligibly small compared to the Kepler velocity v K (i.e., 10 −4 -10 −5 v K ).
In both the protostellar and the circumnuclear disks, the dustgas coupling is characterized by the normalized stopping time, i.e., the Stokes number, S t ≡Ω K t stop , where t stop is the timescale of the dust particles to reach the terminal velocity due to the gas drag. In the Epstein law, where Q g is the Toomre's Q-value for a gas disk and Q g ≡c s Ω K /(πGΣ g ), with the surface density of the gas disk Σ g , and a d,0.1 ≡a d /(0.1 μm), the sound velocity of the gas c s,1 =c s /(1 km s −1 ) and ρ int,1 ≡ρ int /(1 g cm −3 ). Q g <1 is a necessary condition for the ring-mode GI. The radial velocity of the dust v r,d relative to the gas (Weidenschilling 1977;Tsukamoto et al. 2017) is where η is a parameter that determines the sub-Keplerian motion of the gas, where M BH,7 ≡M BH /10 7 M e . Therefore, both in the early stage of the dust evolution (S t =1) and in the late phase (S t ∼ 1) in the circumnuclear disk, v r,d is much smaller than v K ; thus, we can ignore the radial drift of the aggregates in the circumnuclear disk during the whole evolution. The "radial drift barrier," i.e., the dust growth is limited by infalling to the central stars before dust particles obtain large enough mass as planetesimals, is not a serious problem in the circumnuclear disk.

The Growth Time and Destruction by Collisions
The growth time of the aggregate during the hit-and-stick phase can be estimated as in Tsukamoto et al. (2017): 1 2 3 2 in the circumnuclear disk. The Reynolds number R e and the relative velocity of the dust particles Δv are given in the Appendix. This growth time is comparable to the AGN lifetime (see Section 3 and the Appendix in more details).
During the dust compression phase due to collisions, the kinematics of the dust aggregates are dominated by eddies of turbulence in the gas disk. Therefore, relative velocity of the aggregates, Δv, which is important for both growth and destruction of them, is determined by the property of the turbulence, i.e., R e , and S t (Ormel & Cuzzi 2007). The dimensionless parameter α is a parameter to determine the kinematic viscosity in the turbulent disk (Shakura & Sunyaev 1973): Edd,6 ,1 3 whereṀ g is the radial mass accretion rate in the disk, and γ Edd,6 is the Eddington rate for the BH mass with M BH =10 6 M e for the energy conversion efficiency of 0.1.
Here, we assume α is a constant, smaller than unity throughout the disk. In this phase, S t gradually increases, and eventually the phase ends when S t ∼1, then the aggregates decouple from the turbulent gas. At this time, the mass and size of the aggregates become m d ∼10 5 g and a d ∼100 m, respectively. There internal density is still very low (i.e., "fluffy").
In the next gravitational compression phase, the aggregates are compressed owing to their self-gravity, and their internal density increases as r µ m d int 0.4 (Kataoka et al. 2013). The relative velocity of the aggregates in this phase is determined by the energy balance among gravitational scattering, collisional energy loss, turbulent stirring, turbulent scattering, and gas drag (Michikoshi & Kokubo 2016; see also the Appendix). The aggregates finally grow to ∼kilometer-sized bodies (i.e., planetesimals).
The value of the critical velocity for collisional destruction of centimeter-to kilometer-sized dust aggregates is not clear. Numerical experiments of collisions between aggregates (Wada et al. 2009) showed that the critical velocity is 50-100 m s −1 for ∼10 4 monomers (m d ∼ 10 −11 g), and this scales with the mass of the aggregates as µm d 1 4 . If the critical velocity simply scales, it should exceed 1 km s −1 for kilometer-sized "planetesimals." In the regime with S t >1, if the Toomore Q-value for the disk of aggregates becomes smaller than ∼2, the GI takes place (e.g., Goldreich & Ward 1973), and spiral density enhancements are formed, and it leads to rapid growth of planetesimals in a rotational period t K ≡2π/Ω K ;9.5×10 4 yr (M BH / 10 6 M e ) −1/2 (r/1 pc) 3/2 (see also Michikoshi & Kokubo 2017). In fact, we found that the aggregate disks become gravitationally unstable soon after S t reaches unity in most cases.

A Typical Evolution Track toward Planets
We investigated the evolution of icy dust particles based on the processes explained in Section 2 to see if the four evolution stages of the dust aggregates are completed, by changing the parameters, such as the black hole mass M BH , and turbulent viscous parameter α. The circumnuclear cold gas disk embedded in the geometrically thick torus (see Figure 1) is assumed to be gravitationally stable. We assign a constant Toomre's Q-value, i.e., Q g =2 with the gas temperature T g =100 K in the disk. 7 The hydrogen column density of the gas disk is therefore N H ;6.2× 10 23 cm −2 (M BH /10 6 M e ) 1/2 (r/1 pc) −3/2 (T g /100 K) 1/2 . The gas mass between r=0.1 and 10 pc in the thin disk is M g ;5.7×10 4 M e (M BH /10 6 M e ) 1/2 (T g /100 K) 1/2 . The total gas mass in the whole torus system of tens of parsecs could be comparable to M BH , as suggested by recent ALMA observations Combes et al. 2019). Figure 2(a) shows a typical evolution of a dust aggregate at r=5.5 pc, just outside of the snow line (r snow =4.7 pc) around a low luminosity AGN with M BH =10 7 M e and the Eddington ratio, i.e., the luminosity ratio to the bolometric luminosity, γ Edd =0.01. The internal density of the aggregate ρ int is plotted as a function of its mass m d . Here, we assume α=0.1. The internal density decreases monotonically from the monomer's density, ρ 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 ∼0.1 cm. After this hit-and-stick phase, the fluffy dust aggregates keep growing by collisions in the turbulent gas motion until S t ;1. During this phase, the aggregates are compressed by the collisions, and therefore ρ int gradually increase during this phase (m d = 10 −6 g to 10 5 g). At the end of this phase, their size becomes a d ∼1 km. After this stage, the aggregates are compressed by their self-gravity, therefore ρ int increases quite rapidly as seen below. Figure 2(b) plots collisional velocity Δv of the aggregate [cm s −1 ] and its size a d [cm] as a function of m d . a d monotonically increases for S t <1. Initially Δv is 50 cm s −1 and it decreases in the hit-and-stick phase, after which it increases from 10 cm s −1 to ∼500 m s −1 around S t =1 in the collisional compression phase. It is far below the limit of the collisional destruction velocity of aggregates extrapolated from the numerical experiments of collisions between porous aggregates (Wada et al. 2008), which scales with the mass as m d 1 4 . Figure 3 shows time evolution of ρ int , a d , and Δv for the same model in Figure 2. The hit-and-stick phase and the collisional compression phase take ∼3.8×10 8 yr. Soon after the aggregates are decoupled from gaseous turbulence, where S t ;1, their evolution is determined by various heating and dissipation (i.e., cooling) processes in the N-body system of the aggregates (Michikoshi & Kokubo 2017). Note that Figure 3(b) depicts a d grows up to ∼1000 km, but this does not happen because of the GI after S t =1. Figure 3(c) shows that the collisional velocity Δv increases rapidly around S t ∼1, resulting in the rapid growth of the aggregates in the mass and density. The collisional velocity then dramatically decreases , temperature of the gas T g =100 K. The dust sublimation radius is located at r sub =0.3 pc and the snow line is r snow =4.7 pc. This plot for the dust at r=5.5 pc. The aggregates grow by the hit-and-stick process, where the internal density of the aggregates monotonically decreases, that means the aggregates are porous in this phase. After the collisional energy exceeds critical energy, the aggregates start to get compressed (m d > 10 −5 g). The color bar represents the Stokes number. For S t ∼1 the dust aggregates are decoupled from the gaseous turbulence. (b) Collision velocity of the aggregates Δv and size of the aggregate a d as a function of m d . The dashed line shows the limit for the collisional destruction of the aggregates estimated from numerical experiments (Wada et al. 2009). After S t =1 is attained, Δv drops and the disk of the aggregates becomes gravitationally unstable to form more massive "planets," shown by the vertical blue dotted line with "GI" (gravitational instability). 7 Although Q g <1 is the necessary condition for the GI for the m=0 mode perturbation in a thin, uniform disk, the nonaxisymmetric modes can be unstable for Q g 1.5 (e.g., Laughlin & Bodenheimer 1994). Therefore it is safe to assume Q g ∼2 for a gravitationally stable disk (see also Wada & Norman 1999;Wada et al. 2002, where the effective Q-value is 2 in the multiphase, quasi-stable gas disk).
to ∼4 m s −1 due to collisional loss of their kinematic energy. This reduces the Toomre's Q-value of the dust disk, and as a result Q d <2 is attained, therefore the system of kilometer-sized aggregates becomes gravitationally unstable (denoted by the dotted lines with "GI" in Figures 2(a) and (b)). This leads to the formation of spiral instabilities and their fragmentation followed by collapsing massive "planets" consisting of dust. This final unstable phase occurs very quickly with a few rotational periods (i.e., 4 × 10 5 yr at r = 5.5 pc for M BH = 10 7 M e ). The mass of "planets" then would be l Sm M 10 pl d E GI 2 , where λ GI is the most unstable wavelength for the GI, and M E is the Earth mass. In this model, the total number of the "planets" outside the snow line (r=4.7 pc) to r=7 pc is about 8.5×10 4 and its number density is Σ planet ∼10 3 pc −2 .

Discussion
We then explored the evolution of the aggregates by changing α and M BH . In Figure 4, we plot time for which S t becomes unity as a function of α and M BH . The behavior of the dust growth is basically the same as the model with M BH =10 7 M e and α=0.1 (Figures 2 and 3), but the timescale to reach the state with S t =1 depends on α and M BH . For example, for M BH = 10 6 M e and γ Edd =0.01, the snow line is located r=1.4 pc. At r=2 pc, it takes 2.6×10 8 yr when S t exceeds unity. As shown in Equation (5) In other words, formation of "planets" can be expected around the circumnuclear gas disks in low luminosity Seyfert-type AGNs rather than quasar-type high luminosity ones with massive SMBHs.   (Shakura & Sunyaev 1973). After S t ∼1 is attained, the kilometer-sized aggregate system is decoupled from the gaseous turbulence, and it becomes gravitationally unstable, leading to the formation of "planets" within ∼10 6 yr in this parameter range.
The growth time of an aggregate is proportional to α −1/2 . Therefore, if α in the circumnuclear disk is as small as in the protoplanetary disks, i.e., α∼10 −3 -10 −4 , the growth time would be 6×10 8 -2×10 9 yr, which is still smaller than the cosmological time. However, because α is proportional to the mass accretion rate (Equation (6)), the Eddington rate of the central source as a result of the accretion would be γ Edd ∼10 −4 -10 −5 , which corresponds to very low luminosity AGNs (e.g., Ricci et al. 2017). These imply that the planets could be also formed within a few gigayears around very faint AGNs.
The dust-to-gas mass ratio is assumed to be a standard value, i.e., 0.01 (Draine 2011); however, this could be larger by a few factors in the high metallicity environment (e.g., ∼4Z e ) around AGNs (Groves et al. 2006;Rémy-Ruyer et al. 2014). In such a case, the timescale of the dust evolution can be smaller by a few more factors than shown in Figure 4.
Observing planets around SMBHs should be challenging. The standard techniques to detect exoplanets around stars, i.e., Doppler spectroscopy, transit photometry, gravitational microlensing, or direct imaging are hopeless. Photometry by a hard X-ray interferometer in space might be a possible solution, but the occultation of the accretion disk by the "planets" would be hard to distinguish from the intrinsic time variability of AGNs. The other, indirect way is detecting spectral changes in the millimeter-wavelength due to opacity variation associated with the dust growth as used in the protoplanetary disk. The opacity is roughly proportional to ρ int a d , and this increases by more than two orders of magnitude around S t ∼1.
We appreciate the anonymous referee's valuable comments. This work was supported by JSPS KAKENHI grant No. 18K18774. The authors thank Akio Inoue and Tohru Nagao for suggestions on metallicity and the dust-to-gas ratio in AGNs.

2
where H g =c s /Ω K is the scale height of the gas disk.
The relative velocity between aggregates Dv 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) t η =t s =Ω −1 . where Q g is the Toomre's Q-value for the gas disk. The eddy turnover time t L is~Wt L K 1 , and t η ∼t L for the smallest eddy. For the hit-and-stick phase, S t =R e , then Δv obeys the regime I, and where S t,1 and S t,2 are Stokes numbers of two colliding particles. We here assume S t,2 ∼S t,1 /2 (Sato et al. 2016 where v L is velocity of the largest eddy. The size of dust aggregates determines how they interact with the gas (e.g., the Stokes parameter is proportional to ρ int a d for the Epstein law). Dynamics of the aggregates is affected by their cross sections, which depend on the internal inhomogeneous structure. The radius of BCCA cluster a BCCA for large numbers of monomers N is given as a BCCA ;N 0.5 a 0 (Mukai et al. 1992;Wada et al. 2008Wada et al. , 2009, and this was also confirmed by N-body simulations (Suyama et al. 2012).

A.2. Collisional and Gravitational Compression Phases
The hit-and-stick phase ends, when 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 between the porous aggregates. Beyond this point, the

A.2. Collisional and Gravitational Compression Phases
In Equations (15) and (16), the ranges of S t are divided by -Re 1 2 , not Re.

A.3. Evolution of Dust Aggregates as an N-body System
There is a typo in the left-hand side of Equation (17)