Dynamical evidence of the sub-parsec counter-rotating disc for a close binary of supermassive black holes in the nucleus of NGC 1068

It arises a puzzle in \NGC\, how to secularly maintain the counter-rotating disc from $0.2$ to $7\,$pc unambiguously detected by recent ALMA observations of molecular gas. Upon further analysis of disc dynamics, we find that the Kelvin-Helmholtz (KH) instability (KHI) results in an unavoidable catastrophe of the disc developed at the interface between the reversely rotating parts, and demonstrate that a close binary of supermassive black holes provides tidal torques as the unique external sources to prevent the disc from the KH catastrophe. We are led to the inescapable conclusion that there must be a binary black hole at the center of NGC 1068, to prevent it from the KH catastrophe. The binary is composed of black holes with a separation of $0.1\,$pc from GRAVITY/VLTI observations, a total mass of $1.3\times 10^{7}\:M_{\odot}$ and a mass ratio of $\sim 0.3$ estimated from the angular momentum budge of the global system. The KHI gives rise to forming a gap without cold gas at the velocity interface which overlaps with the observed gap of hot and cold dust regions. Releases of kinematic energies from the KHI of the disc are in agreement with observed emissions in radio and $\gamma$-rays. Such a binary is shrinking with a timescale much longer than the local Hubble time via gravitational waves, however, the KHI leads to an efficient annihilation of the orbital angular momentum and speed up merge of the binary, providing a new paradigm of solving the long term issue of"final parsec problem". Future observations of GRAVITY+/VLTI are expected to be able to spatially resolve the CB-SMBHs suggested in this paper.


INTRODUCTION
For many decades, the Seyfert 2 galaxy NGC 1068 serving as an archetype for unification model of active galactic nuclei (Antonucci 1993) has received the most intensive observations from radio, infrared, optical, ultraviolet, X-ray and γ−ray bands, which had been individually reported by several hundreds of papers. It is brighter than 9.1 mag with a redshift of z = 0.003793 and a distance of 14.4 Mpc). However, the nature of the NGC 1068 nucleus hidden by the dusty torus still remains open so far.
ALMA targeted CO (6 − 5) molecular line and resolved a 7 − 10 pc geometrically thick-disc structure in NGC 1068 (Garcia-Burillo et al. 2016;Gallimore et al. 2016). The dynamics shows non-circular motions and enhanced turbulence superposed on a surprisingly slow rotation pattern of the disc (Garcia-Burillo et al. 2016;Gallimore et al. 2016). However, ALMA recent observations of HCN gas spatially resolved a counter-rotating disc (CRD) from 0.5 − 1.2 pc and beyond reversely extended to 7 pc (Imanishi et al. E-mail: wangjm@ihep.ac.cn 2018; Impellizzeri et al. 2019). Though the interface of gas density map between the prograde and retrograde parts is not distinguished with ALMA's current spatial resolutions (1.4 pc for NGC 1068), but the velocity interface is quite clear in the velocity map of HCN gas (Impellizzeri et al. 2019). In particularly, it has been found that the outer disc shows Keplerian rotation consistent with extrapolation from the inner disc (Impellizzeri et al. 2019).
It is very motivated that the CRD formation could be related with the two streamers at ∼ 10 pc scale developed from the south and the north tongues by SINFONI/VLT (Müller Sánchez et al. 2009;Gravity et al. 2020). Random collisions of molecular clouds from large scale distances drive them into nuclear regions and trigger Seyfert activity (Sanders 1981). Two tongues around the central black hole are formed by this way, and later could lead to formation of the CRD in NGC 1068. It is obvious that the CRD in NGC 1068 is unavoidably suffering from the well-known Kelvin-Helmholtz instability (KHI) so that it will be destroyed within one Keplerian orbit. This directly conflicts with the presence of the observed CRD, what mechanism is secularly maintaining it at least as long as the viscosity timescale remains a puzzle.  (R 1.5 pc) and the retrograde (3 pc R 7 pc) parts in NGC 1068 observed by ALMA (Impellizzeri et al. 2019) and MIDI/VLTI (Raban et al. 2009). R CBD is the inner edge of the CBD measured by GRAVITY (Gravity et al. 2020), which is consistent with dust sublimation radius. The prograde part is the NIR emission regions where water maser clouds (green circles) are co-spaced. The retrograde part is the FIR emission regions. Color of the red represents the dusty molecular disc (HCN, HCO + and CO) radiating in infrared. The interface regions in purple between the prograde and the retrograde parts are undergoing the Kelvin-Helmholtz instability (KHI) to form shocks and turbulences, driving formation of a gap with a width of δR ≈ 0. 82M H 0 . The KHI dissipates the kinematic energy of the counter rotating disc giving rise to bremsstrahlung emissions observable in radio with morphology tightly related to the disc shapes. Moreover, shocks in the KHI layer accelerate some electrons to be relativistic to radiate γ-ray emissions observed by Fermi-LAT.
In this paper, we demonstrate that the CRD could be only secularly maintained by a close binary of supermassive black holes (CB-SMBHs) supplying angular momentum to the CRD avoiding the KH catastrophe. A model based on observations is provided for the binary system with a counter-rotating circumnuclear disc. Such a system is showing evidence for an efficient way of removing the angular momentum through the retrograde accretion to harden the binary.

THE CIRCUMNUCLEAR DISC SUPPORTED BY
BINARY BLACK HOLES Figure 1 shows the geometric structure of the CRD obtained from ALMA, VLTI/MIDI, GRAVITY, VLBI and BLVI (see details in Appendix A). The CRD of HCN molecular gas in NGC 1068 is composed of the inner counterclockwise (prograde: 0.5−1.2 pc) and the outer clockwise rotation (retrograde: < 7 pc) parts distinguished by the interface radius around R * ≈ 2.5 pc revealed from the velocity map (Impellizzeri et al. 2019). The black hole mass has been measured by water maser dynamics (Greenhill et al. 1996(Greenhill et al. , 1997, but corrected by some effects of massive disc gravity (Hure et al. 2002(Hure et al. , 2011Lodato & Bertin 2003). However, the mass of disc within a few parsec is still dominated by the central black hole (Imanishi et al. 2018), questioning the validity of the corrections. Neglecting the uncertainties, we take the black hole mass of ∼ 10 7 M in this paper, which should be robust. On the other hand, NGC 1068 is a prototypical oval galaxy with an unusually massive pseudobulge (Kormendy & Ho 2013) of M bulge = 10 10.9±0.1 M and stellar dispersion velocity of σ e = 151 ± 7 km s −1 . Usually, pseudobulges are located below the well-known M • − σ relation (Tremaine et al. 2002), the upper limits of the black holes are M • < 3 × 10 7 M from this relation.
The KHI risen by shear collisions is unavoidable at interface regions, making the CRD lifetime just comparable with the Keplerian rotation timescale of t K = Ω −1 K = 1.3 × 10 4 M −1/2 7 R 3/2 2pc yrs through efficient annihilation of angular momentum (AM) simply estimated by the linear analysis (Quach et al. 2015), where R 2pc = R/2 pc is the radius of disc and M 7 = M • /10 7 M is the black hole mass. This is the KH catastrophe of the CRD in NGC 1068. For a thin disc keeping a vertical equilibrium, it secularly evolves with the viscosity timescale (Pringle 1981(Pringle , 1991 2pc yrs, where α = 0.1α 0.1 is the viscosity parameter, h = H 0 /R = 0.1 h 0.1 is the relative height of the disc. Observations show that the disc of HCN gas follows Keplerian rotations in the both prograde and retrograde parts (Impellizzeri et al. 2019), indicating that the two parts are in dynamical equilibrium and therefore already exist at least as long as a lifetime of t vis .

Kelvin-Helmholtz instability
The CRD is supersonically rotating. The Mach number of disc rotation is where v K = RΩ K is the rotation velocity and c 0 s is the sound speed outside the KHI layer. In this paper, the subscripts of ⊕ and symbols denote the prograde and retrograde rotations, respectively. In NGC 1068, GRAVITY observations show h 0.14 in the dust sublimation ring, but h ≈ 0.3 for the R 1.2 pc regions (near infrared regions with temperatures of 800−1500 K), namely, M ⊕ ∼ 3. In the mid-infrared regions ( 3 pc regions with temperatures 300 K) (Raban et al. 2009), we have h ≈ 0.6, as a geometrically slim, and M ≈ 1.7. Being different from the well-understood normal KHI modes (Drazin 2002), both the prograde and the retrograde parts as supersonic flows lead to an extremely complicated interface from AM cancellations subject to black hole gravity and fast cooling of the hot plasma heated by the formed shocks. We take the simplest approach using the Rankine-Hugoniot conditions of shocks for the supersonic collisions of two streams in a slab geometry (Lee et al. 1996). Both the prograde (v 1 = v K ) and retrograde (v 2 = −v 1 ) parts have the same Keplerian velocity (v K ) but they have densities of ρ 0 and χρ 0 ( χ ≤ 1) at the interface, respectively, where χ is the density ratio depending on propterties of the CRD. The velocity, density and pressures of the post shock regions are where Γ is adiabatic index. For the case with identical densities ( χ = 1), we have the most efficient dissipation of the kinematic energy into the thermal since v = 0. This gives rise to multiwavelength emissions tightly related to the dynamical structure of this region.
Supposing the interface radius (R KHI ) with a thickness of δR, the AM cancellation driven by the KHI is happening with a timescale of t KHI ≈ Ω −1 K and δR ≈ t KHI v R , where v R is the velocity of gas entering the KHI zone. In a linear analysis, V R is given by the viscosity stress, which is much smaller than the sound speed (c s ) in the KHI layer formed due to shear collisions. In the present context, however, we have v R ≈ c s in the layer because the efficient dissipation of the kinematic energy of the interface layer heats the medium and the waves propagate with c s . With the shock conditions, we have c 2 s = p/ρ = 2 χ(1 + √ χ) −1 (Γ − 1)v 2 K = (2/3)v 2 K , for χ = 1 and Γ = 5/3, which results in the KHI gas to have two thirds of the virial temperatures. Considering H 0 = c 0 s Ω −1 K , we have the layer width approximated by δR ≈ t KHI c s ≈ 0.82 Ω −1 K v K = 0.82 MH 0 . For the current case of NGC 1068, the interface of the density map is not resolved by the current ALMA observations (Impellizzeri et al. 2019), but the velocity interface is quite clear from the velocity map (see Figure 2 in Impellizzeri et al. (2019)). The projected radius is about R * ≈ 2.5 pc, which exactly overlaps with the geometrical gap between the NIR and MIR component (Raban et al. 2009). We thus take the interface radius R KHI = R * . For the simplest estimation, we just take M KHI ≈ 2.5 before shear collisions, which is the averaged value between the NIR and MIR regions. The KHI layer therefore has a width of δR ≈ 2 H 0 ≈ 1.6 pc, where H 0 ≈ 0.8 pc is the averaged height of the NIR and MIR regions in light of observations (Raban et al. 2009). This gap should be observed in density map, however, it is not spatially resolved (1.4 pc for NGC 1068) of the current ALMA (Impellizzeri et al. 2019). However, we note that the width of the gap is affected by the prograde part supported by the tidal torques of one CB-SMBH.
The KH catastrophe can be evaluated by AM distributions of the CRD. The total AM of the retrograde part (within 7 pc) can be estimated by J ≈ −M gas v KR ≈ −325.3 M 1/2 7 M pc 2 yr −1 , whereR = 5 pc is the averaged radius of the retrograde regions (3-7 pc), and v K ≈ 10 7 cm s −1 at this radius. For the prograde part (< 1.5 pc), we have the total AM of J ⊕ ≈ 190.9 M 1/2 7 M pc 2 yr −1 , where we take M ⊕ gas = 9 × 10 5 M and the averaged radius R ⊕ = 1 pc. Owing to |J ⊕ | < |J |, the prograde part will be destroyed by the retrograde via the KHI action. Given an AM distribution of a disc as J ∝ R β (β > 0), the cancelled AM can be estimated by δJ /J = β (δR/R KHI ). From δR ≈ 0.82MH 0 , we have δJ /J ≈ 0.82β. Though the gaseous dynamics has been measured in NGC 1068, the spatial distribution of mass density is poorly done. In the 14 pc × 10 pc regions (Imanishi et al. 2018), the mass is about 2 × 10 6 M corresponding to an averaged surface density of Σ gas ≈ 1.4 × 10 4 M pc −2 over this region. We have Σ gas ∝ R −0.6 from the prograde (R ⊕ , Σ ⊕ given in SM) to retrograde (R , Σ ) parts, giving rise to J ∝ R 1.9 . With β ≈ 1.9, we find that the AM of the CRD can be cancelled completely by once action of the KHI within a dynamical timescale of t KHI ∼ Ω −1 K = 1.9 × 10 4 M −1/2 7 (R KHI /R * ) 3/2 yrs. This is a very short timescale compared with the typical AGN lifetime, even the flickering lifetime (∼ 10 5 yrs) (Schawinski et al. 2015). It would be worse for the prograde part if the retrograde part is able to obtain more AM from the outer boundary. The KH catastrophe cannot be avoided unless there is an extra supplement from external torque sources. Obviously, the only way to maintain the prograde part is the gain AM from a CB-SMBH.

Close binary of supermassive black holes
As seen in a brief review in Appendix B, evidence for a CB-SMBH is still quite elusive in observations though there are a few candidates. Starting from one CB-SMBH with a circular orbit of a separation (A) composed of the primary (M p ) and the secondary (M s ), we have its period of P orb ≈ 877 (1 + q) −1/2 M 7 a 3/2 5 yrs, and orbital AM of G is the gravitational constant and c is the speed of light. The tidal torques of one CB-SMBH strongly interact with its circumbinary disc (CBD) to form a cavity with an inner radius of R CBD = 2A in light of analytical calculations and simulations (Artymowicz et al. 1994;Armitage et al. 2002). In Appendix C, we demonstrate that 2µm-imaged hole discovered by GRAVITY (Gravity et al. 2020) is approximate to R CBD = 0.24 pc with help of the gas dynamics of ALMA-detected CRD. The timescale of gravitational waves is given by t GW = 2.0 × 10 14 M 7 a 4 5 /q(1 + q) yrs for circular orbits (Peters 1964) much longer than the Hubble time, which is the wellknown "final parsec problem". To maintain the prograde part, we need J B J ⊕ , yielding q 0.3 (for M 7 = a 5 = 1) in order to make up for the AM loss due to the KHI action.
The CBD regions governed by the tidal torques of the CB-SMBH can be estimated in a simple way. Tidal forces can be approximated by F tid ≈ GM tot δm/R 2 (A/R) for a point mass of δm, which is about perpendicular to the direction to the center, and the corresponding torques are T tid ≈ F tid R = δm(1 + q)c 2 ar −2 /4, where r = R/R Sch . For the mass with angular momentum of Considering the rapid decays of the tidal torques with radius (Papaloizou & Lin 1984; Pringle 1991), we recast the tidal timescales of t tid = ε 0 (1 + q) −1/2 (a/r) −4 t K , where ε 0 ≈ 0.1. In the tidal-governed regions t tid ≤ t vis , we have the critical radius of R tid ≤ 10 (ε 0. where ε 0.1 = ε 0 /0.1, beyond which the viscosity torques becomes dominant. For NGC 1068, we have R tid ≈ 1.2 pc for typical values (ε 0.1 , α 0.1 , h 0.1 ) = 1, which is in agreement with the sizes of the prograde part (Impellizzeri et al. 2019). Such a CB-SMBH guarantees the stability of the prograde part and prevents it from the KH catastrophe.
The gas disc beyond R tid is jointly governed by viscosity stress and the tidal torques, and its outer boundary conditions. Since R > R tid part is retrograde with respect to the CB-SMBH orbits, it loses AM due to the tidal torques, enhancing accretion onto the center (Nixon et al. 2011;Roedig & Sesana 2014;Bankert et al. 2015).
In the context of the inner part as a decretion disc (Pringle 1991), its outer boundary moves outward with the viscosity timescale but encounters the inward flows from the retrograde part, leading to an equilibrium radius (R KHI ) between the inward and the outward flows. Dynamical evolution of the CRD could compress the width of the gap. Numerical simulations are expected for details of such a complicated case. For NGC 1068, J B > |J ⊕ +J | indicates that the CB-SMBH with a total mass of 1.3 × 10 7 M , mass ratio of q = 0.3 and a separation of 0.1 pc can support this observed CRD in the galactic center. On the other hand, the CB-SMBH in NGC 1068, as the first observed object, faces the well-known "final parsec problem" (Begelman et al. 1980;Goicovic et al. 2018), that the binary black holes are not able to merge within local Hubble time. Successive series of random accretion of molecular clouds onto binary black holes have been suggested by simulations (Goicovic et al. 2018), on the other hand, from observational side, two tongues from south and north direction around the center appear in NGC 1068 ∼ 10 pc regions (Müller Sánchez et al. 2009). The retrograde part can efficiently cancel the binary orbital AM through the KHI of random and episodic accretion. The episodic numbers can be roughly estimated by N = J B /|J + J ⊕ |, and accretion with low AM (i.e., |J + J ⊕ | ≈ 0) can efficiently shrink the CB-SMBH in light of t GW ∝ M −3 p for the case given A. For NGC 1068, we have N ≈ 4.8 from ALMA observations (Imanishi et al. 2018;Impellizzeri et al. 2019). This indicates that the retrograde accretion onto CB-SMBHs is an efficient way of removing orbital AM of the binaries.

Multi-λ emissions
During the AM cancellation through the KHI, a huge amount of kinematic energies of the CRD is released, which can be estimated by ∆E K ≈ ∆M KHI v 2 K within the interval of Ω −1 K . The radiative luminosity from the dissipation is given by .1 (taking r KHI = 1 and r * = R * /R Sch = 5 × 10 6 ) and M torus is the accretion rates of the torus. For the typical values of the parameters, we have L KHI ≈ 1.1 × 10 44 M 1 ergs s −1 at the peak frequency of peak ∼ 0.1 r −1 KHI keV from the hot KHI layer during the interval of the KHI timescale, where M 1 = M torus /10 M yr −1 . The infall with a typical rate of 10 M yr −1 at a few parsec scales in NGC 1068 has been observed by SINFONI (Müller Sánchez et al. 2009;) (most of this rate will be channeled into outflows). Such a powerful emission can be detectable in type 1 AGNs, however, for NGC 1068 as a Compton thick AGN (Zaino et al. 2020), the major emissions from the KHI layer will be undetectable because of absorption of the torus along the line of sight. The KHI provides a new explanation of soft X-ray excess if there is a CRD like in NGC 1068. The excess should be relatively stationary because of its origins from a large scale.
On the other hand, fortunately, radio emission from this HKI dissipation is detectable, which can be estimated by the bremsstrahlung emission from the KHI layer. The shear collisions of the two streamers make the post-shocked gas with temperatures close to the virial one as T gas ≈ 1.2 × 10 6 r −1 KHI K, where r KHI = R KHI /R * . Since the bremsstrahlung cooling timescale of hot plasma is (Rybicki & Lightman 1979) t cool ≈ 3.0 n −1 6 T 1/2 6 yrs, the KHI layer can be efficiently cooled, where n 6 = n e /10 6 cm −3 is the density (about the CRD gas density) and T 6 = T/10 6 K is the layer temperature. The fast cooling (t cool Ω −1 K ) prevents the shocks from propagation outside the layer. Owing to the black hole gravity, the cooled layer will rapidly collapse onto the prograde part with the free-fall timescale of t ff ≈ 1.8 × 10 4 r 3/2 KHI yrs. On the other hand, gas is supplied beyond the KHI layer with a timescale of t vis ≈ 10 3 α −1 0.1 h −2 0.1 t K by the viscous stress from the retrograde part. The comparison of t cool t KHI ∼ t ff t vis implies that the gas supply to the layer will be halted. The cooled gas in the layer will be emptied by the gravity of the central black holes, forming a gap with a width of δR. Fortunately, such a gap has been revealed by MIDI/VLTI, showing spatial separations between the NIR (hot dust with temperature > 800 K, and located < 1.4 pc) and MIR (< 300 K, and located > 3 pc) emissions (Raban et al. 2009).
With the KHI layer's volume of V KHI = 2πR KHI δRH 0 = 1.6πR 3 KHI M (H 0 /R) 2 , we have radio emissions of L 256GHz = 1.0 × 10 39 n 2 6 T −1/2 6 R 3 2.5 h 2 0.1 ergs s −1 at ν = 256 GHz according to the thermal emissivity [Eq.5.14b in Rybicki & Lightman (1979)]. This radio emission agrees with 12.7 mJy observed by ALMA (Impellizzeri et al. 2019) at 256 GHz. As emitted from the layer, the radio morphology should follow the velocity interface. Indeed, this agrees with the observations by comparing the radio 256 GHz map with the HCN velocity map from Figure 2 in Impellizzeri et al. (2019), showing the north-east part of the HCN disc and southwest part where the velocity interfaces matches the radio morphology. Moreover, the layer has a luminosity of L 5 GHz ≈ 2.0×10 37 ergs s −1 at 5GHz from bremsstrahlung emissions, which roughly agrees with the compact radio source S 1 in the nuclear center.
Additionally, the strong shocks unavoidably accelerate some electrons to be relativistic in the KHI layer, leading to nonthermal high energy emissions from the layer. Here the dissipation of the KHI layer could provide an alternative explanation. Magnetic fields are of B = 10 −3 Gauss in this region (Krolik & Begelman 1987). The shock velocity is V sh ∼ 3 × 10 7 T 1/2 6 cm s −1 , if Γ = 5/3, yielding the acceleration timescale of t acc = R L c/V 2 sh = 5.8 × 10 6 γ 5 B −1 −3 T 6 sec, where γ 5 = γ/10 5 is the Lorentz factor of relativistic electrons, R L is the Lamore gyration radius (Blandford & Eichler 1987). The energy density of infrared photons is about U IR = 1.8 × 10 −6 L 43 R −2 pc ergs cm −3 dominating over the magnetic fields, where L 43 = L IR /10 43 ergs s −1 is the infrared luminosity and R pc = R/1 pc is the distance of the accelerating regions to the center. In such a context, the cooling of the electrons is governed by inverse Compton scatter and given by t IC = 3m e c/4σ T γU IR = 1.7 × 10 8 γ −1 5 L −1 43 R 2 pc sec (Rybicki & Lightman 1979), and the relativistic electrons contribute less in radio band. The maximum Lorentz factor is γ max = 5.5 × 10 5 L −1/2 The maximum energetic photons are max ≈ γ 2 max IR = 74 (γ max /γ 0 ) 2 0.25 GeV, where γ 0 = 5.5 × 10 5 , 0.25 = IR /0.25 eV is IR photon energy in units of λ = 5µm estimated from the spectral energy distribution (Gravity et al. 2020). If the total budget of γ-rays is conservatively 1% [see Blandford & Eichler (1987)] of L KHI , about 10 42 ergs s −1 is produced that is enough to explain the observed γ-rays.

Magnetic fields
We first consider the role of a strong azimuthal magnetic field to support the CRD. The Alfvén velocity is given by where B φ is the azimuthal velocity, µ 0 is permeability of free space and ρ 0 is the density of gas. The plasma couples with neutral part tightly so that the magnetic fields controlling plasma govern the neutral part of gas. If the azimuthal magnetic field supports the CRD, the Alfvén waves provide enough energy to secularly support the CRD, namely v 2 A ≥ (t vis /t K ) (∆V) 2 , where ∆V ≈ 2v K is the different velocity of the CRD. The magnetic field is given by where ρ −18 = ρ 0 /10 −18 g cm −3 [the ALMA-observed density of gas from the HCN line (Imanishi et al. 2018)] and v 2 = v K /10 2 km s −1 . Recent polarization observations show B ≈ 0.7 mG in the few parsec regions of its circumnuclear region (Lopez-Rodriguez et al. 2020) in NGC 1068, which is much lower than the lower limit of the magnetic fields to support the CRD. This directly rules out the magnetic fields to support the unusual CRD.

A CRD around a single SMBH
Since nuclear regions are much smaller than the circumnuclear disk, random accretion onto black holes is naturally taking place (Sanders 1981;. In a cosmic timescale, such a kind of the random accretion leads to spin-down evolution of SMBHs as a result of cancellations of spin AM from different episodes (Kuznetsov et al. 1999;King et al. 2008), but speeds up cosmic growth of SMBHs due to low radiation efficiency. According to the cosmic evolution of duty cycle of SMBH activity episodes , the spin-down behaviors of SMBH evolution are found from application of cosmic equation of radiative efficiency to the survey data (Wang et al. 2009;Li et al. 2012), providing strong evidence for the random accretion. In particular, cold clumps are universal in galaxies (Tremblay et al. 2016;Temi et al. 2018). On the other hand, this indicates that there were a series of CRDs across cosmic time. The CRD will give rise to very fast accretion onto SMBHs [also see Ref (Kuznetsov et al. 1999;Quach et al. 2015;] and show some imprints at scale of a few parsec. Is there any independent evidence for this? Thanks are given to recent ALMA observations of nearby Seyfert galaxies for a new clue to understanding this issue. Among this sample, three Seyfert galaxies (NGC 1365, NGC 1566 and NGC 1672) show a gas hole or depletion in the center within a few parsec (Combes et al. 2019), of which two are more nearby than NGC 1068. It is very interesting to note that NGC 1365 and NGC 1566 are changing-look AGNs. NGC 1365 is changing its column density from Compton-thick to -thin (Winter et al. 2009) whereas NGC 1566 is changing spectral type (Oknyansky et al. 2019). We preliminarily speculate that the gas hole could be caused by the KHI, which leads to AM annihilation of CRDs depleting the parsec scale regions. If this happened, a fast collapse onto very central regions results in high accretion rates of AGNs. The parsec scale hole of gas could be an indicator of the remanent of a past CRD around a single SMBH. Considering the short lifetime of a CRD around single SMBH as low as t KHI /t vis ∼ 10 −3 , we think that the CRD will be detected in an extremely low opportunity. Future MUSE observations of ALMA-observed targets should explore detailed structure at a few parsec scales to test if there are something related to "two tongues" in NGC 1068 (Müller Sánchez et al. 2009). Additionally, we should check their polarized spectra for optical Fe emission lines for accretion status of the central black holes (Du & Wang 2019). On the other hand, candidate AGNs harboring CB-SMBHs with CRDs can be selected from Seyfert 2 galaxies whose polarized spectra show asymmetric profiles of broad emission lines.

CONCLUSIONS AND REMARKS
The puzzle of a counter-rotating disc from 0.2 − 7 parsec regions discovered by ALMA observations arises a serious issue how to secularly maintain it. The Kelvin-Helmholtz instability will destroy the disc unless there is an extra supply of angular momentum. A close binary of supermassive black holes is expected to support the unusual disc for a viscosity timescale. The binary black holes are of a total mass of 1.3×10 7 M and mass ratio of 0.3 and a separation of 0.1 parsec. With annihilation of angular momentum due to the KHI, the disc efficiently dissipates its kinematic energy into radiation peaking at soft X-rays and radio bands from bremsstrahlung emissions. Though soft X-rays are absorbed by the torus, the radio emissions are in agreement with observations. We also find significant γ-rays from the shocks developed from the KHI, which agrees with Fermi and MAGIC observations.
The polarized spectra of NGC 1068 show asymmetric profiles of broad Hβ (with polarization degree of ∼ 3%) (Miller et al. 1991) indicating complicated structure of the broad-line regions. In the NIR, the core has a m K ≈ 8.7 mag (Gravity et al. 2020), the polarized magnitudes could be of ∼ 12.5 mag in K−band. The polarized fluxes of broad Brackettγ line (Martins et al. 2010) could be too weak to detect through the GRAVITY/VLTI. Fortunately, it is strong enough for GRAVITY+ as the next generation of GRAVITY/VLTI to make efforts to spatially resolve the < 0.1 parsec regions for the appearance of the CB-SMBHs through differential phase curves (Songsheng et al. 2019).
Given the black hole mass, mass ratio and separations, the CB-SMBH is radiating gravitational waves (GW) with frequencies of f GW = 0.1 (1 + q) 1/2 M −1 7 a −3/2 5 nHz and intrinsic strain of amplitudes is h s = 5 × 10 −19 qM 7 a −1 5 d −1 14.4 , where d 14.4 = d/14.4 Mpc is the distance of NGC 1068 to observers. GW background composed of NGC 1068-like AGNs can be detected by Square Kilometer Array (SKA) (Barack et al. 2018). However, the retrograde part eventually shrink the CB-SMBH through cancelling its orbital AM to merge, providing a paradigm of solving the "final parsec problem", which is related with many concomitant phenomenon risen from the KHI. The CB-SMBH in NGC 1068, as one of long-sought-after candidates, has an angular size of 2 mas, which is larger by a factor of 100 than the spatial resolution of the Event Horizon Telescopes (EHT), and is thus worth observing via EHT.

APPENDIX A: OBSERVATIONAL PROPERTIES
Polarized spectra. The polarized spectra of NGC 1068 show typical spectra of Seyfert 1 galaxies (Antonucci 1985;Miller et al. 1991) and hence demonstrate the existence of an obscured broadline region (BLR) by a dusty torus. In light of the scenario, an orientation-based unification scheme was suggested for all Seyfert galaxies as well as for quasars (Antonucci 1993). Broad Hβ line appears in the optical polarized spectra along with strong and broad Fe emission lines. Using the polarized spectra, the black hole mass can be estimated by the virialized motion of BLR. We find the black hole mass of (8.4 ± 0.4) × 10 6 M from the 5100 Å luminosity estimated from the [O ] line (Du et al. 2017), which is consistent with results estimated from 2-10 keV luminosity (Zaino et al. 2020).
It is important to note that the polarized Hβ profile is asymmetric (Miller et al. 1991), implying complicated structures of the hidden BLR, at least existence of sub-structures of ionized gas (Du et al. 2018). Since the polarized lights to observers are from scatters of an electron screen vertically located at a distance of 100 pc from the galactic center, it will be very hard to observe variations of the polarized Hβ line through a feasible length campaign. Actually, the Brackett γ line profiles are double-peaked, see Figure 4 in Martins et al. (2010). NGC 1068 is bright enough in near infrared for for GRAVITY spectroastrometry to measure the differential phase curves. Signatures of CB-SMBHs could appear in the differential phase curves (Songsheng et al. 2019).
Dusty discs. The most important information for infrared continuum is the extinction, which can be estimated by near infrared emission lines. The Brackett α broad line component observed by Infrared Space Observatory (ISO) constrains a lower limit of the extinction A Brα,4.05 ≥ 2.4, implying the extinction A K ≥ 6 in the K-band (Lutz et al. 2000). Observations of mid-infrared interferometric instrument (MIDI) for NGC 1068 reveal two spatially resolved regions: 1) inner and hot (> 800 K) zone with 1.35 pc long and 0.45 pc thick at a PA = −42 • , and 2) outer and warm (∼ 300 K) zone with 3 × 4 pc (Jaffe et al. 2004;Raban et al. 2009) but extended to 7 pc regions along the north-south axis (López-Gonzaga et al. 2014), and somehow may correspond to polar dust emission in the ionization cone. The gap between the NIR and MIR dust regions can be approximated by R * ≈ 2.5 pc, actually overlap with the velocity interface of the HCN gas map (Impellizzeri et al. 2019). We thus adopt it as the interface radius where the KHI is taking place.
GRAVITY onboard Very Large Telescope Interferometer (VLTI) observations find the dust sublimation radius R sub = 0.24 pc, where hot dust particles (∼ 1500 K) are responsible for NIR and MIR emissions with extinctions of foreground dust (A K = 5.5) (Gravity et al. 2020). A structure with a hole has been found around 0.24 pc (Gravity et al. 2020). This radius agrees with the empirical relation established by NIR dust reverberation mapping of AGNs (Minezaki et al. 2019). However, the innermost part of the dusty disc is geometrically thin as H/R 0.1, which is still much thicker than the Shakura-Sunyaev disc.
Maser disc. Water (and OH) masers observed by VLBI trace inner edge of the torus from a radius of 0.4 to ∼ 1 pc with a PA = −45 • completely different from radio jet (Greenhill et al. 1996(Greenhill et al. , 1997Gallimore et al. 1996b,c). The inclination of the maser disc is ∼ 80 • whereas an inclination of ∼ 70 • is obtained by fitting infrared continuum (Lopez-Rodriguez et al. 2018). Generally the maser disc aligns with the dusty disc though slightly orientated differently.
Radio morphology. A radio core source, called S 1 , associated with the disc of H 2 O vapor megamaser emission has been found (Greenhill et al. 1996;Gallimore et al. 2001). It is interesting to find that S 1 overlaps with the inner prograde HCN disc (Impellizzeri et al. 2019). ALMA observations of 256 GHz continuum show very interesting features, not discussed in Ref (Impellizzeri et al. 2019), that the radio flux contours overlap with the velocity interface between the two counter rotating parts of the HCN disc, see Figure  2 in Ref (Impellizzeri et al. 2019).

APPENDIX B: CLOSE BINARY OF SUPERMASSIVE BLACK HOLES
As expected for several decades from the theory of the merger tree of galaxy evolution, CB-SMBHs (defined as those with separations less than 0.1 pc) evolved from dual galactic cores must be located in centers of some galaxies sometimes (Begelman et al. 1980;Volonteri et al. 2003). Most phenomenon of AGN activities can be explained by accretion onto a single SMBH (Rees 1984;Osterbrock & Mathews 1986), implying that most CB-SMBHs have finished the final mergers. Single epoch spectrum with double-peaked profiles cannot be used as a direct diagnostic of CB-SMBHs (Popovic 2012; Gaskell 2010). Though dual AGNs are quite common among galaxies (Comerford et al. 2013;Fu et al. 2015;Liu et al. 2018), however, observational evidence for CB-SMBHs is very elusive (Popovic 2012; Dotti et al. 2012). Though the radio map of NGC 7624 (z = 0.0029) shows two radio cores with a separation of 0.35 pc in the nucleus and plausibly implies a candidate binary black holes (Kharb et al. 2017), none of CB-SMBHs are identified so far. A very brief summary is given below about the progress of searching for CB-SMBHs.
The final parsec problem. The formation of massive black hole binaries is unavoidable as a outcome of successive mergers of galaxies, however, it turns out that the decay of the binary orbits in a real galaxy would be expected to stall at separations of parsec scales unless some additional mechanism is able to efficiently remove the orbital AM (Milosavljević & Merritt 2001). At large separations, dynamical friction of the binaries with background stars controls orbital evolution (Milosavljević & Merritt 2001;Yu 2002;Wang & Yuan 2012), but becomes inefficient when the mass enclosed in the orbits are comparable in gas-poor environment and hardly evolve into stages of radiating gravitational waves. This is the socalled "final parsec problem" (Begelman et al. 1980;Milosavljević & Merritt 2001).
Numerous simulations show that the gas-rich environment can efficiently absorb and transfer the orbital AM through the tidal interaction with the prograde CBD (Escala et al. 2005;MacFadyen & Milosavljevi 2008;Cuadra et al. 2009;Lodato et al. 2009), or the retrograde (Nixon et al. 2011;Roedig & Sesana 2014;Bankert et al. 2015). However, the real situations are complicated by successive random accretion of clumpy clouds formed in circumnuclear regions (Tremblay et al. 2016;Temi et al. 2018), which have been found as an efficient clue to harden binary black holes from simulations (Goicovic et al. 2017(Goicovic et al. , 2018. Fortunately, such a configuration of accretion onto binary black holes has been found in NGC 1068 by ALMA (Impellizzeri et al. 2019), which provides unique opportunity to investigate the orbital evolution predicted by simulations.
Signatures of CB-SMBHs. The suggested signatures so far are: 1) periodicity of long term radio or optical variations of AGNs from a few years to a few tens of years, which is regarded as results of modulation of the orbital motion of CB-SMBHs; 2) ultraviolet deficit of continuum rising from the central cavity of the CBD governed by the tidal torques of the binary; 3) shifts of red and blue peaks of broad emission lines in long term variations similar to normal binary stars; 4) deficit profile of brightness of galactic centers as a result of ejection of stars through interaction of stars with the binary black holes; 5) 2-dimensional kinematic maps of broad emission lines in AGNs showing signatures of two point potentials generated through reverberation mapping campaigns; 6) GRAVITY (or Extremely Large Telescope: ELT) with spectroastrometry may reveal some signals of double BLRs. Additionally, X-shaped radio jet was regarded as a signature of CB-SMBHs (Merritt & Ekers 2002), but it hardly applies to measure orbital parameters of the binaries.
Great efforts have been made to search for periodicity of AGN long-term variations over several decades. Nowadays, a sample composed of about 100 AGNs was built up through Catalina Real-time Transient (Graham et al. 2015a), among which PG 1302-102 is the best one with a period of 5 yrs (Graham et al. 2015b;D'Orazio et al. 2015). Other three objects with periodical variations are OJ 287 in radio (Valtonen et al. 2008), NGC 5548 with a period of about ∼ 13 yrs over 40 yrs (Li et al. 2016), and Akn 120 with a period of ∼ 20 yrs over the last 40 yrs . Progress in this way is quite slow so far unless the future LSST can efficiently discover some AGNs with short periods from the fast large scale of time domain surveys.
Continuum features of CB-SMBHs radiated from accretion result from the cavity of the CBD governed by tidal interaction (Gültekin & Miller 2012). However, the shapes of UV continuum are easily affected by dust extinctions (Leighly et al. 2016), making that it is hardly useful to justify candidates of CB-SMBHs in practice. Moreover, it has been shown that the binary black holes are peeling off gas from the inner edge of the CBD making the accretion rates quite high (the rates still comparable with the accretion rates of the CBD) (Farris et al. 2014;Shi & Krolik 2016;Bowen et al. 2019). For the current case of NGC 1068, the binary black holes are peeling off gas for accretion to radiate gravitational energies.
Searching for systematically periodical shifts of AGNs with double-peaked profiles is expected to find CB-SMBHs (Runnoe et al. 2017;Doan et al. 2020). Similar to classical binary stars, the double BLRs orbiting around the mass center of the binary black holes will lead to opposite shifts of the red and blue peaks with orbital phases. The long term campaigns are expected to carry out in next a few years, but SDSS J0938+0057, SDSS J0950+5128, and SDSS J1619+5011 in this sample showed systematic and monotonic velocity changes consistent with the binary hypothesis (Runnoe et al. 2017).
As a consequence of interaction with CB-SMBHs, as shown by numerical simulations (Ebisuzaki et al. 1991;Milosavljević et al. 2002;Merritt 2006), stars are ejected from galactic center, forming a "mass deficit" in its central part. Actually, it does not appear to have a pronounced and unambiguously detectable effects on the density profile. Even though it works for searching for CB-SMBHs, on the other hand, this method doesn't apply to determine orbital parameters to test properties of gravitational waves, and the same shortcomings for the identifications of AGNs with X-shaped jet morphology.
Reverberation mapping of AGNs is a powerful tool to probe the kinematics and structure of the BLR (Peterson 1993). Usually, BLRs have relatively stable structures justified from the large sample of SDSS quasars (Vanden Berk et al. 2001) though they are variable on timescale of years. In principle, the broad-line gas is governed by the central black hole potential, and therefore it is expected to show signatures of CB-SMBH orbital motion from the 2-dimensional velocity delay map from reverberation mapping of AGNs Songsheng et al. 2020;Kovacěvić et al. 2020). This reasoning also depends on the orbital motion, but the advantage of this scheme is obvious that the campaigns can be performed in seasons and does not require long-term monitoring as the approach of period searching. A campaign called as Monitoring AGNs with Hβ Asymmetry (MAHA) has been conducted with Wyoming Infrared Observatory (WIRO) Telescope targeting on objects with various Hβ profiles (Du et al. 2018). Results are expected to be carried out shortly.
GRAVITY/VLTI provides the highest spatial resolution in NIR. The differential phase curves measured by GRAVITY depends on the geometric structures and kinematics of ionized gas emitting broad lines. For CB-SMBHs, there are several simplest configurations of AM distributions of individual BLRs and binary orbital motion. Details of the phase curves have been explored for signals of CB-SMBHs (Songsheng et al. 2019). We are expecting to jointly observe some Seyfert 2 galaxies through ALMA and GRAV-ITY, which can be selected by checking optical polarized spectra with complicated profiles. We hope discover more Seyfert 2 galaxies with the CRDs like in NGC 1068.
Some AGNs show that rotation curves of masers don't follow the Keplerian law. The radial self-gravity of the maser disc could be important (Hure et al. 2011), however, infrared emissions doesn't support this hypothesis. Moreover, ALMA measurements of molecular gas mass in NGC 1068 is significantly smaller than the black hole. It is a curiosity that the maser disc shares different dynamics from the molecular gas if inspecting details (Greenhill et al. 1997;Impellizzeri et al. 2019). On the other hand, it could be plausible, if speculate, that CB-SMBHs work in these objects. Detailed dynamical modeling of the maser or molecular gas at parsec scales will allow us to determine some of orbital parameters, which are keys to test properties of nano-Hz gravitational waves.

APPENDIX C: CAVITY OF THE CBD
Dynamics of gas at parsec scales revealed by the water megamaser (Greenhill et al. 1996;Gallimore et al. 1996a) and HCN/HCO + (J = 3 − 2) molecular gas (Imanishi et al. 2018;Impellizzeri et al. 2019; provides strong constraints on the accretion disc size of the central black hole and hence spatial distribution of gas in the compact regions. For NGC 1068, the black hole mass is of M • = (0.8 − 1.7) × 10 7 M and the bolometric luminosity is L Bol ≈ 2.0 × 10 45 ergs s −1 from multiwavelength data (Gravity et al. 2020), we thus have the Eddington ratio of λ Edd ≈ 1, implying that NGC 1068 has high accretion rates. The strong Fe emissions in the polarized optical spectra (Miller et al. 1991) shows the key feature of high rates (Du & Wang 2019).
Recent ALMA observations of HCN/HCO + (J = 3 − 2) emission lines (Impellizzeri et al. 2019; show that the extrapolation of the rotation curve of the 0.5 − 1.2 pc HCN molecular gas, which agrees with that of the water maser (Greenhill et al. 1996;Gallimore et al. 1996b), is consistent with the rotation curve extended to 7 pc disc (Impellizzeri et al. 2019). This indicates that the enclosure mass is still dominated by the central black hole at a few parsec scales (Hure et al. 2002). With help of the gas disc mass (M AD ) estimated from Equation (D2), we have upper limit of its outer boundary of R M • ≤ 1.95 × With the gas dynamics, we draw a conclusion that the disc must be truncated at least R M • . Indeed, the HCN molecular gas is M ⊕ gas ∼ 9 × 10 5 M ( 3 pc: corresponding to a surface density Σ ⊕ gas ≈ 3 × 10 4 M pc −2 and a number density of n ⊕ ≈ 2.2 × 10 6 cm −3 ) measured by ALMA (Imanishi et al. 2018;, which is smaller than M • . At this edge of the disc, its temperature is T eff (R M • ) ≈ 120 ( m 10 /M 7 ) 1/4 R/R M • −3/4 K much lower than the dust sublimation temperature, and the disc surface density is Σ AD ≈ 5.7 × 10 7 M pc −2 . There is another critical radius corresponding to the dust sublimation temperature, we have R dust = 0.04 ( m 10 /M 7 ) 1/3 R M • . According to the Kennicutt-Schmidt law (Kennicutt 1998), this high surface density of gas between (R dust , R M • ) with low temperatures are undergoing starbursts with a rate of R ≈ 4.8 M yr −1 , resulting in an infrared luminosity of L * MIR ≈ 1.1 × 10 44 ergs s −1 radiating between 8-100µm. Such a bright mid-and far infrared emission doesn't appear in the spectral energy distributions (the observed MIR luminosity is L obs MIR 1.0 × 10 43 ergs s −1 estimated from Figure 5 in Ref (Gravity et al. 2020)). Considering the constraints of L * MIR < L obs MIR , we have the outer boundary radius R AD L obs MIR /L * MIR R M • 0.1 R M • = 0.019 pc, which is consistent with R dust . This means that the R AD − R M • region should be a cavity of gas (or a low-density region though we are not able to estimate its density currently). On the other hand, GRAVITY (Gravity et al. 2020) measured the radius of the dust sublimation R sub ≈ 0.24 ± 0.03 pc and discovered an extended structure with a central hole with the radius of R sub by imaging the central regions of NGC 1068 at 2µm, agreeing with R M • R sub . All these provide evidence for that the 2µm-imaged hole is a central cavity of gas.
Such a cavity is a natural consequence of tidal interaction of the circumnuclear disc (CND) with the binary black holes (Artymowicz et al. 1994;Escala et al. 2005;MacFadyen & Milosavljevi 2008;Cuadra et al. 2009;Lodato et al. 2009) if they exit subsequently evidenced by the CRD in the galactic center of NGC 1068. According to analytical analysis, the inner radius of the CND, which is denoted the circumbinary disc (CBD), should be around R CBD ≈ R sub = 2A, where A is the separations of the binary black holes. From the 2µmimaged hole, we have A ≈ 0.12 pc.

APPENDIX D: ACCRETION DISCS
Given the black hole mass of M • = 10 7 M in NGC 1068, distance scales are: 0.1 pc ≈ 114 ltd ≈ 10 5 R Sch . The dimensionless accretion rate is defined as m = M/L Edd c −2 = 10 η −1 0.1 λ Edd from L Bol = η Mc 2 , where η 0.1 = η/0.1 is the radiative efficiency, L Edd is the Eddington luminosity, λ Edd = L Bol /L Edd is the Eddington ratio and M is the accretion rates. In the outer part of the Shakura-Sunyaev disc (Shakura & Sunyaev 1973), where gas pressure dominates over radiation pressure and absorption over Thompson scattering, the surface density is given by