© B. Vollmer et al. 2022

1. Introduction

Active galactic nuclei (AGNs) can outshine the optical light of the whole galaxy which surrounds them. This huge luminosity and thus this huge amount of energy is gained through the conversion of gravitational potential energy to electromagnetic radiation. The ingredients for an AGN are the follwing: (i) a supermassive black hole with a mass exceeding about 10⁶ M_⊙; (ii) interstellar gas that surrounds the central black hole; and (iii) gas infall onto the black hole. As the material falls in toward the black hole, angular momentum causes it to form an accretion disk. During infall the disk gas is heated to temperatures between 1 and 4 × 10⁵ K (e.g., Bonning et al. 2007), giving rise to intense optical and UV emission. Moreover, the geometrically thin accretion disk is surrounded by a hot corona giving rise to strong X-ray emission. The typical size of the accretion disk is about 10 lightdays or ∼0.01 pc (e.g., Kokubo 2018). Further out, at distances of a few parsec to a few tens of parsecs, a massive gas disk is present in AGNs (Sani et al. 2012; Hicks et al. 2013; Combes et al. 2019). The estimated disk masses are on the order of 10% of the dynamical mass, exceeding 10⁶ M_⊙ (García-Burillo et al. 2021). Because of the uncertain conversion factor between luminosity and mass, the mass estimates are quite uncertain. Most frequently, these massive gas disks are geometrically thick, that is they have a high gas velocity dispersion (Davies et al. 2007; Hicks et al. 2009; Sani et al. 2012).

The Galactic Center is not an AGN although it hosts a supermassive black hole (M_BH = 4.3 × 10⁶ M_⊙, Gillessen et al. 2009) and a gas disk – Circumnuclear Disk (CND) – of several 10⁴ M_⊙ located between 1 and 5 pc (e.g., Guesten et al. 1987; Christopher et al. 2005; Etxaluze et al. 2011). The reason for the quiescence of the Galactic Center is a negligible mass accretion rate onto the black hole.

With a distance of 8.2 kpc (GRAVITY Collaboration 2019), the gas of the CND can be observed with a resolution of ∼0.1 pc. The CND in the Galactic Center is made of gas clouds with sizes of ∼0.1 pc, masses of a few 10 M_⊙, an area filling factor of Φ_A ∼ 0.1, and a volume filling factor of Φ_V ∼ 0.01 (Jackson et al. 1993; however, Requena-Torres et al. 2012 suggest a much higher volume filling factor of Φ_V ∼ 0.2). Since the densities derived from the HCN lines (Jackson et al. 1993) are close to those of selfgravitating clouds with a thermal sound speed of ∼1 km s⁻¹, Vollmer & Duschl (2001a) and Vollmer et al. (2004) proposed a model of a collisional disk made of stable gas clouds. This model was extended to AGN gas tori in Vollmer et al. (2008). The small linewidths of the selfgravitating gas clouds are in stark contrast to the observed large linewidths in the CND (Montero-Castaño et al. 2009). This means that the clouds are not selfgravitating and thus short-lived. Strong gas turbulence can create the observed large velocity dispersion.

Vollmer & Davies (2013) developed an analytical model for turbulent clumpy gas disks where the energy driving turbulence is supplied by external infall or the gain of potential energy by radial gas accretion within the disk. The gas disk is assumed to be stationary (∂Σ/∂t = 0) and the external mass accretion rate is thought to be close to the mass accretion rate within the disk (the external mass accretion rate feeds the disk at its outer edge). Within the model, the disk is characterized by the disk mass accretion rate Ṁ and the Toomre Q parameter which is used as a measure of the gas content of the disk. It is suggested that the velocity dispersion of the gas disk at a 10-pc scale or torus is increased through adiabatic compression by the infalling gas. The turbulent and collisional model were applied to the CND in the Galactic Center and thick obscuring tori in AGNs. Despite the different gas masses (CND: 10⁴ M_⊙; AGN torus: 10⁶ M_⊙) and velocity dispersions (CND: 20 km s⁻¹; AGN torus: 50 km s⁻¹), both models share disk gas clouds of similar masses and sizes.

With the advent of the Atacama Large Millimeter/submillimeter Array (ALMA), it is now possible to observe the centers of other galaxies with a decent resolution (≲1 pc) in the millimeter regime. NGC 1068 is the archetypical type 2 AGN with a massive black hole of 8 × 10⁶ (Lodato & Bertin 2003), where the massive thick gas disk is seen edge-on and thus it entirely obscurs the central engine.

Imanishi et al. (2016, 2018) observed the center of NGC 1068 in the HCN(3−2) line with a resolution of 0.1″–0.2″ (7–14 pc) and 0.06″ (4 pc), respectively. García-Burillo et al. (2016) observed the same target in the CO(6−5) line with a 0.06″ (4 pc) resolution. These authors detected a massive gas disk of one to a few 10⁵ M_⊙ and a size of ≲10 pc. Subsequent higher resolution observations (García-Burillo et al. 2019; Impellizzeri et al. 2019; Imanishi et al. 2020) with beamsizes on the order of 20 mas (1.4 pc) clearly resolved the massive gas disk and gave access to its detailed kinematics. Compared to its appearance in the HCN(3−2), HCO⁺(3−2), and HCO⁺(4−3) lines, the molecular disk is bigger and more lopsided in the lower-J CO transitions (García-Burillo et al. 2019; Impellizzeri et al. 2019). Moreover, Impellizzeri et al. (2019) revealed a counter-rotating inner disk at radii smaller than ∼1.2 pc (see also Imanishi et al. 2020). This is surprising because counter-rotating gas disks are supposed to be unstable against shear instability (Quach et al. 2015).

Gallimore et al. (2016), García-Burillo et al. (2019), and Impellizzeri et al. (2019) found signs of a molecular outflow in the inner region of the massive gas disk (R ≲ 3 pc). Furthermore, the torus is connected to the 200 pc-size gas ring through a network of gas lanes whose kinematics are accounted for by a 3D outflow geometry (García-Burillo et al. 2019). The latter authors argued that about half of the mass of the massive gas disk is outflowing. Gas compression by the outflow or by radiation pressure may be considered as sources of mechanical heating of the molecular gas at the inner edge of the massive disk.

In this article we continue the comparison of Vollmer & Davies (2013) between the CND in the galactic center and the massive gas disk in the center of NGC 1068. In the following, the term CND is exclusively used for the 1–5 pc gas disk or ring in the Galactic Center. We produced CND-like and NGC 1068-like dynamical models and calculated their CO, HCN, and HCO⁺ emission for multiple transitions. Whereas traditional models derive the gas properties from molecular line ratios, our forward modeling aims at reproducing qualitatively and, to a certain extent, quantitatively the observed molecular brightness temperatures.

The article is structured in the following way: the model is described in Sect. 2, its application to the CND in the Galactic Center in Sect. 3 and to the inner 20 pc of NGC 1068 in Sect. 4. The cosmic ray (CR) ionization fractions, a more quantitative comparison between the model and observations, the role of the continuum emission, and the relevance of the different gas-heating mechanisms are discussed in Sect. 5. Finally, we give our conclusions in Sect. 6.

2. The model

The inner ∼50 pc of galaxies are complex environments involving a dense nuclear disk and star cluster and a central massive black hole. Misalignments between the ionization cones of AGN, and hence the obscuring gas structures, and the disks of the host galaxies are frequently observed (e.g., Fischer et al. 2013). Indeed, the nuclear gas angular momentum content can be significantly different from that of the galactic disk because of large-scale gas fragmentation or secondary stellar bars (Hopkins et al. 2012). Therefore, gas disks or rings with arbitrary inclinations and position angles are expected at scales of ∼10 pc. These rotating gas structures are perturbed by infalling gas from larger distances with a different angular momentum. The simplest scenario is that of a gas cloud falling onto a preexisting gas disk or ring. Even in such a simple configuration there are more than ten free parameters and parameter space is vast. Therefore, we do not aim at reproducing details of the available observations of the Galactic Center and NGC 1068. Our work has a broader scope. We want to find out if such a simple scenario can explain the main characteristics of the observations in terms of gas distribution and kinematics. In a second step, the combination of the large-scale dynamical model with the small-scale analytical model of turbulent clumpy accretion disks allows us to calculate the molecular line emission of the largest turbulent disk gas clouds. It is then possible to investigate if our models can reproduce the characteristics of the available observations in terms of gas density, temperature, and chemistry, clump area filling factor, and brightness temperatures of the molecular lines. In this respect, our models are not expected to exactly reproduce the available observations, but to serve as a guideline to interpret them. With this work we want to make a step forward in the understanding of the large-scale gas dynamics and the physics in terms of turbulence, heating, cooling, and chemistry in the inner ∼30 pc around a central supermassive black hole.

Our modeling of thick gas disks located at distances between 1 and 10 pc from the central black hole comprises five components: (i) a sticky-particles dynamical code for the large-scale gas dynamics, where the particles represent gas clouds, (ii) an analytical model of a turbulent clumpy accretion disk from which the properties (density, size, velocity dispersion) of the largest turbulent gas cloud are derived, (iii) an analytical model for the calculation of the disk cloud temperature via the equilibrium between gas heating and cooling, (iv) a model for the gas chemistry, and (v) the calculation of the molecular line emission from the disk gas clouds. The dynamical model is used to follow the pc-scale evolution of the gas distribution within a distance of ∼30 pc around the central black hole. The dynamical simulations of a collisional gas disk can be regarded as an approximation of a simulation of a turbulent gas disk in terms of large-scale gas density, velocity dispersion, and gas viscosity (Vollmer & Davies 2013). The analytical accretion disk model, which is used to calculate the disk cloud properties, is linked to the dynamical model via the total gas mass and the particle velocity dispersion: the analytical model assumes a Toomre Q parameter and a gas mass of the central gas disk or ring, which is consistent with the gas properties of the timestep of interest of the dynamical model and observations. Furthermore, we made sure that the disk mass accretion rate of the dynamical model is also consistent with that of the analytical model. Given the size, density, and velocity dispersion of a disk gas cloud, its temperature was calculated via the equilibrium between turbulent mechanical and cosmic ray heating and radiative cooling. The molecular abundances were determined for each disk gas cloud independently by the two-phase gas-grain code Nautilus (Hersant et al. 2009). The cloud molecular gas density, surface density, temperature, and velocity dispersion were then used to calculate the CO, HCN, HCO⁺, CS, and CN line emission and datacubes were produced from the simulations. These datacubes are directly compared to the available observations.

2.1. The dynamical model

Following Vollmer & Duschl (2002), the ISM is simulated as a collisional component, namely as discrete particles that possess a mass and a radius and can have partially inelastic collisions. In contrast to smoothed particle hydrodynamics (SPH), which is a quasi-continuous approach where the particles cannot penetrate each other, our approach allows a finite penetration length between particles, which is given by the mass-radius relation of the particles. During the disk evolution, the cloud particles can have partially inelastic collisions, the outcome of which (coalescence, mass exchange, or fragmentation) is simplified following the geometrical prescriptions of Wiegel (1994).

We follow the orbits of these clouds in the three dimensional gravitational potential. The radial distribution of the total enclosed mass M(R) is given by

(1)

where M_BH is the mass of the central black hole, and M₀ describes the mass distribution of the stellar content. We used M_BH = 4 × 10⁶ M_⊙ and M₀ = 1.6 × 10⁶ M_⊙ pc^−5/4 (Genzel et al. 2010) for the CND in the Galactic Center and M_BH = 10⁷ M_⊙ (Greenhill et al. 1996) and M₀ = 0.8 × 10⁶ M_⊙ pc^−5/4 for NGC 1068.

The integration of the ordinary differential equation is done with the Burlisch-Stoer method (Stoer & Bulirsch 2002) using a Richardson extrapolation and Stoermer’s rule. Vollmer & Duschl (2002) chose the error level to match the theoretical collision rates. The global timestep¹ is typically around 50–100 yr. For a velocity of 100 km s⁻¹ this corresponds to ∼0.01 pc. During each cloud-cloud collision the overlapping parts of the clouds are calculated. Let b be impact parameter and r₁ and r₂ the radii of the larger and smaller clouds. If r₁ + r₂ > b > r₁ − r₂ the collision can result into fragmentation (high-speed encounter) or mass exchange. If b < r₁ − r₂ mass exchange or coalescence (low speed encounter) can occur. In this way a cloud mass distribution is naturally produced. The cloud masses and velocities resulting from a cloud-cloud collision are calculated by assuming mass and momentum conservation.

As in Vollmer & Duschl (2002) we allowed for additional energy dissipation by lowering the final gas cloud velocities by a constant fraction ξ: where v_ini and v_end are the velocity vectors of the initial and resulting gas clouds after the collision. The energy dissipation rate of a collisional disk is

(2)

where t_coll = 4 r_cl/(3 Φ_Vv_turb) (Vollmer et al. 2008) is the collision timescale, Φ_V the cloud volume filling factor, and Σ the mean gas surface density of the disk, that of a turbulent gas disk is

(3)

where the turbulent driving length equals the disk height l_driv = H (Vollmer & Davies 2013). Both expressions are equivalent if the dissipated energy fraction per collision is

(4)

For simplicity we used the cloud radii from the analytical model if r_cl ≤ 0.05 pc and r_cl = 0.05 pc elsewhere. In this way our simulated disk is broadly consistent with η = 0.1 and thus ξ = (1 − η) = 0.9 for the NGC 1068-like simulations (Fig. 1).

Fig. 1. Fraction of dissipated energy per cloud-cloud collision as a function of the distance to the central black hole. The dissipation energy fraction is shown for all clouds during the NGC 1068 model evolution presented in Fig. 3.

By construction, the dynamical model does not include a wind or outflow. Such a wind is observed in NGC 1068 in the line emission of ionized gas (the narrow-line region) at scales of 10–100 pc (e.g., May & Steiner 2017), which has an hourglass structure. Most importantly, its kinematics are consistent with an outflow within a hollow-cone structure (Das et al. 2006; Miyauchi & Kishimoto 2020). At smaller scales, the morphology of the point-source subtracted ALMA band 6 continuum emission (Impellizzeri et al. 2019) and the polar dust emission from MIR interferometric observations (López-Gonzaga et al. 2014) are also consistent with an outflow. Gallimore et al. (2016), García-Burillo et al. (2019), and Impellizzeri et al. (2019) interpreted the high radial velocities and large linewidths of the molecular line emission along the minor axis of the molecular torus as a molecular outflow in the inner region of the massive gas disk (R ≲ 3 pc). Furthermore, the torus is connected to the 200 pc-size gas ring through a network of gas lanes whose kinematics are accounted for by a 3D outflow geometry (García-Burillo et al. 2019). Radiation pressure, the radio jet, and the high-velocity ionized outflow certainly help to launch the molecular outflow. A magnetocentrifugal molecular and dusty wind starting at the inner edge of the massive gas disk is also viable (Vollmer et al. 2018). For a discussion of the outflow scenario we refer to Sect. 5.4.

The cloud particle masses of our dynamical simulations range between 1 and 90 M_⊙. The mean mass of the clouds is ∼20 M_⊙. For the simulations corresponding to the CND in the Galactic Center the number of gas clouds was N ∼ 3200, and the total gas mass is M_gas, tot = 7.5 × 10⁴ M_⊙. For the simulations of the gas distribution around the central black hole as in NGC 1068 the number of gas clouds was N ∼ 1.6 × 10⁴ and the total gas mass is M_gas, tot = 3.1 × 10⁵ M_⊙.

For all simulations the preexisting gas ring was located in the x − y plane. A massive approximately spherical gas cloud was located in the same plane with initial conditions that lead to a prograde orbit and apocentric distances of 10, 15, and 20 pc and pericentric distances of 2, 3, and 4 pc for the CND and 3, 4, and 5 pc for NGC 1068. In this way nine models were calculated. The infalling gas cloud was then rotated around the x-axis from β = 20° to 340° in steps of 20° leading to a total number of 162 models for the CND. Because we focused on counter-rotating gas infall in NGC 1068, the infalling cloud was rotated around the x-axis from β = 120° to 240° in steps of 20° leading to a total number of 63 models for NGC 1068. For each model timestep a datacube was produced with the following projections in steps of 10°: CND: 0 ° ≤i ≤ 350°, 0 ° ≤PA ≤ 350°, 0 ° ≤az ≤ 350°; NGC 1068: 40 ° ≤i ≤ 140°, 220 ° ≤PA ≤ 320°, 0 ° ≤az ≤ 350°, where i is the inclination angle, PA the position angle, and az the azimuthal angle within the x − y plane. We had to restrict the number of projections of the NGC 1068 model because of the long computation times.

The model datacubes were directly compared to the IRAM 30m CS(2−1) observations of Güsten et al. presented in Vollmer et al. (2003) for the CND and to the ALMA CO(2−1) observations of García-Burillo et al. (2019) for NGC 1068. The best-fit models were determined by searching for the minimum reduced χ². For the CND, we additionally required the western edge to be closer to the observer (e.g., Liu et al. 2012). Since the environment of the CND is complex, it is not possible to reproduce the gas distribution with a single infalling gas cloud and therefore, the χ² are quite high. Our timestep of interest has χ² = 161. We considered timesteps with up to a 10% higher χ² as acceptable. To better reproduce the gas distribution of the inner 20 pc, we decided to add a second acceptable timestep to our final CND model. Among the acceptable timesteps we selected the timestep with the best-fit kinematics by setting all voxels of the model and observed datacubes to unity for the χ² calculations. The reduced χ² of the combined model (χ² = 148) is smaller than those of all single timesteps. The properties of the two models and the combined model are presented in Table 1.

The evolution of the CND-like is shown in Fig. 2.

Fig. 2. Evolution of the CND-like model. The times since the beginning of the simulations are indicated. Upper half: RA-Dec projection. Lower half: RA-LOS (line-of-sight) projection.

The impact of the first infalling cloud onto the CND occurs at about t = 0.23 Myr, that of the second infalling cloud at about t = 0.65 Myr. The timestep, which has the minimum χ² and which we compare to observations, that is the time of interest, is t = 0.8 Myr. At this timestep the gas mass of the CND within a radius of 4 pc is 1.2 × 10⁴ M_⊙.

The evolution of the NGC 1068-like model is shown in Fig. 3.

Fig. 3. Evolution of the NGC 1068-like model. The times since the beginning of the simulations are indicated. Upper half: RA-Dec projection. Lower half: RA-LOS (line-of-sight) projection.

The impact of the infalling cloud occurs at about t = 0.33 Myr, the time of interest, where the reduced χ² is minimum, is t = 0.47 Myr. At this timestep the gas mass located between 1.5 and 7 pc is 1.4 × 10⁵ M_⊙. The CS(2−1) emission distribution of the combined CND-like model is presented in Fig. 4, the CO(2−1) emission distribution of the NGC 1068-like best-fit model in panel c of Fig. 16.

Fig. 4. CS(2−1) emission distribution. Left panel: IRAM 30 m observations of the CND in the Galactic Center (Güsten et al., priv. comm.). Right panel: best-fit combined model (see Table 1). The levels are 13 to 260 K km s−1 in steps of 13 K km s−1. Both images have the same physical size.

We inspected the 150 model snapshots with the smallest χ² by eye. All CND-like model snapshots with i > 50° where the western edge is the near side of the CND are prograde encounters between the gas ring and the infalling gas cloud (0 ° ≤β ≤ 60° or 320 ° ≤β ≤ 360°). Furthermore, retrograde encounters generally lead to much higher mass accretion rates than prograde encounters. The higher mass accretion rates lead to centrally concentrated gas distributions, which are not observed in the Galactic Center. We thus conclude that the gas clouds interacting with the CND are most probably on prograde orbits with respect to the CND. On the other hand, only a retrograde encounter can explain the observed kinematics in the inner 20 pc of NGC 1068.

2.2. The analytical model of a turbulent clumpy gas disk

We assume a turbulent clumpy gas disk where the energy to drive turbulence is supplied by external infall or the gain of potential energy by radial gas accretion within the disk. The analytical model is fully described in Vollmer & Davies (2013). An overview of the model is given in Appendix A. Within the model, the disk is characterized by the disk mass accretion rate Ṁ and the Toomre Q parameter, which is used as a measure of the gas content of the disk. Both parameters are assumed to be constant within the disk. The size of the largest turbulent gas clouds is determined by the size of a continuous (C-type) shock propagating in dense molecular clouds with a low ionization fraction at a given velocity dispersion. We used the expressions derived by Vollmer & Davies (2013) for the expected volume and area filling factors, mass, density, column density, and velocity dispersion of the disk clouds. The latter is based on scaling relations of intermittent turbulence whose open parameters are estimated for the Circumnuclear Disk in the Galactic Center. The turbulent clumpy gas disks of the CND and NGC 1068 are calibrated in terms of the radial profiles of the gas surface density and velocity dispersion (Fig. A.1) by the dynamical model and observations.

2.3. The cloud temperature

For the determination of the temperature of the turbulent disk gas clouds we follow the model of Vollmer et al. (2017), which is based on the calculations of Neufeld & Kaufman (1993) and Neufeld et al. (1995). Their model for the radiative cooling of molecular gas includes a detailed treatment of the interstellar chemistry that determines the abundances of important coolant molecules, and a detailed treatment of the excitation of the species H₂, CO, H₂O, HCl, O₂, C, O, and their isotopic variants where important. For simplicity, we only take the main cooling agents, CO, H₂, and H₂O, into account.

For the calculation of the thermal balance within molecular clouds one needs to consider processes affecting the gas and the dust in addition to the radiative gas cooling. We assume gas heating via turbulence, X-rays, and cosmic rays:

(5)

Photoelectric heating by UV photons within photodissociation regions is neglected (for a discussion of photodissociation regions in NGC 1068 see Sect. 4.2.4). The turbulent heating is

(6)

(e.g., Mac Low 1999). Following Maloney et al. (1996) the X-ray heating rate at the cloud center is

(7)

where L_X is the X-ray luminosity between 1 and 100 keV, and N_H = n_Hr_cl is half the column density of the cloud. Following Nelson & Langer (1997) we adopted the cosmic ray heating rate

(8)

where ζ_CR is the cosmic ray ionization rate.

2.4. Chemical network

Chemical modeling was carried out using the Nautilus gas-grain code presented in detail in Hersant et al. (2009), Semenov et al. (2010), and Ruaud et al. (2015). This code computes the abundances of chemical species (atoms and molecules) as a function of time by solving the rate equations for a network of reactions. For gas-phase reactions, we used the kida.uva.2014 network (Wakelam et al. 2015²) comprising 489 species and 6992 reactions. For grain surface reactions, we used the desorption, diffusion, activation barrier energies along with a set of grain surface reactions, all from the KIDA database. Both, thermal and nonthermal desorption processes are taken into account, the latter consisting mainly of CR-induced desorption following the formalism presented by Hasegawa & Herbst (1993).

The model parameters are time, density, gas temperature, grain temperature, UV flux, cosmic ray ionization rate, and the elemental abundances of the elements C, O, N, and S (C/H = 1.7 × 10⁻⁴, O/H = 2.4 × 10⁻⁴, N/H = 6.2 × 10⁻⁵, S/H = 10⁻⁶). At the beginning of the calculation all hydrogen is in molecular form.

Grids of models were obtained by varying the disk cloud lifetime (20 log spaced steps between 10³ and 10⁸ yr), the cloud density (20 log spaced steps between 10³ and 10⁹ cm⁻³), and cloud gas temperatures (20 log spaced steps between 10 and 1000 K). For each type of cloud, the CO, HCN, and HCO⁺ abundances were interpolated on the grid given the lifetime, density, and temperature of the cloud.

2.5. The X-ray dominated region

In the vicinity of the central engine, the interstellar medium is exposed to a tremendously strong UV and X-ray radiation field. Both can ionize the gas and dissociate the molecules comprised within the gas. About 50% and 10% of the bolometric luminosity are emitted in the UV and X-rays. Close to the central engine the dust is sublimated if it is heated to temperatures in excess of ∼1500 K. In NGC 1068 this is the case at a distance of 0.25 pc (GRAVITY Collaboration 2020). Inside this radius, the hydrogen can be still in molecular form if its column density is high enough to permit self-shielding. Outside the dust sublimation radius the UV emission is absorbed by dust with an optical depth of τ_UV ∼ N_H/(10²¹ cm⁻²). The regions, which are directly illuminated by the central UV emission, are called photodissociation regions (see, e.g., Tielens & Hollenbach 1985). These regions generally include a hot (T > 100 K) atomic region near the surface with an extinction A_V < 3, a warm (T ∼ 100 K) partially dissociated region at about A_V ∼ 3 − 5, and a cooler (T < 100 K) interior region at A_V ∼ 10 where oxygen is still photodissociated to atomic form. Only the X-ray emission can penetrate into region of higher A_V. In these regions the X-rays dominate the heating of the gas and its chemistry through X-ray ionization and dissociation. Following Maloney et al. (1996), the energy deposition rate per particle at the center of a gas cloud is

(9)

where L_X is the X-ray luminosity and N_H the hydrogen column density. In the low-ionization limit, the X-ray ionization rate is

(10)

For the chemical network we replace the cosmic ray ionization rate by the X-ray ionization rate whenever the latter exceeds the former.

For an X-ray luminosity of L_X ∼ 5 × 10⁴³ ergs s⁻¹ in NGC 1068 (Bauer et al. 2015; Marinucci et al. 2016), a distance of 3 pc, and a column density of N_H = 10²⁴ cm⁻², we obtain ζ_T = 8 × 10⁻¹¹ s⁻¹. This is about a factor 4 × 10⁴ higher than the cosmic ray ionization rate in the Circumnuclear disk in the Galactic Center (Yusef-Zadeh et al. 2013; Harada et al. 2015). The Nautilus model HCN abundances as a function of the gas density and temperature are presented in Fig. 5 for ionization rates between 2 × 10⁻¹⁵ s⁻¹ and 6 × 10⁻¹¹ s⁻¹. It becomes clear that HCN abundances in excess of x_HCN = 10⁻⁸ are absent for ionization rates that are higher than a few 10⁻¹² s⁻¹. The reason for this absence is the destruction of H₂ molecules through photodissociation, which blocks the chemical reactions leading to the formation of the HCN molecule. Thus, the observed high HCN abundances (Imanishi et al. 2020) cannot be reproduced with models including an ionization rate which is expected if the gas is directly illuminated by the X-ray emission of the central engine. The same result is found for the HCO⁺ and CO abundances.

Fig. 5. HCN abundance xHCN = nHCN/nH as a function of density and temperature for different cosmic ray ionization fractions. Panel a: ζCR = 2 × 10−13 s−1, Panel b: 2 × 10−12 s−1, and panel c: 2 × 10−11 s−1.

Since an X-ray ionization rate of ζ_T ≳ 8 × 10⁻¹¹ s⁻¹ is expected, we conclude that the bulk of the central X-ray emission must be absorbed in a gas layer of column density N_H ≳ 10²⁵ cm⁻², a Compton-thick gas layer. The same conclusion was reached by Burtscher et al. (2016) who stated that deviations from the Galactic N_H/A_V can be simply explained by dust-free neutral gas within the broad-line region in some sources.

If the absorbing gas has a radial extent of ∼0.2 pc, these column densities () can be reached in the midplane of the molecular gas distribution if the inner wall of a thick disk or ring is located at R₁ ≲ 2.6 pc for Q = 10, R₁ ≲ 1.5 pc for Q = 30, and R₁ ≲ 0.7 pc for Q = 100. Alternatively, the Compton-thick gas layer can be located inside the dust sublimation radius and therefore be dust-free. The latter possibility is in agreement with the conclusions of GRAVITY Collaboration (2020). These authors stated that large column densities (∼10²⁵ cm⁻²) are found in the X-rays, which may largely originate from dust-free plasma inside the dust sublimation radius at distances smaller than ∼0.25 pc. In the following, we assume that the X-ray emission is entirely absorbed by the dense gas inside or outside the dust sublimation radius and that the X-ray emission does not play a role for the heating and ionization of the molecular gas.

2.6. Molecular line emission

We employ the emission line modeling used by Vollmer et al. (2017). For simplicity we consider only a single collider (H₂). We consider two-level molecular systems in which the level populations are determined by a balance of collisions with H₂, spontaneous decay and line photon absorption, and stimulated emission with τ > 1. The molecular abundances were calculated using the chemical network (see Sect. 2.4). For simplicity, we neglected the hyperfine structure of HCN. A more detailed description of our modeling is given in Appendix B.

2.7. Model datacubes

A brightness temperature and a cloud surface area were assigned to each gas cloud particle of the dynamical simulations according to the assumed gas density, temperature, velocity dispersion, and size obtained from the analytical model (Sects. 2.2–2.6). A first datacube with a voxel size of 0.01 pc in the two spatial axes and a voxel size of 10 or 20 km s⁻¹ according to the available observations was established and the gas clouds were placed into this datacube according to their projected two-dimensional position and their radial velocity. Within this first datacube the vast majority of the gas clouds is spatially resolved. A Gaussian line profile was applied to the clouds according to the analytical model. The emission from gas clouds, which were hidden (spatially and by radial velocity) by other gas clouds whose optical depth exceeds unity, was removed from the datacube. Each of the velocity channels was then convolved to the spatial resolution of the corresponding observations and the datacubes were clipped with the rms of the available molecular line observations (see Table 2 for NGC 1068). Moment maps and position-velocity (pv) diagrams were calculated from the model datacubes.

3. A CND-like model

The underlying dynamical simulation is shown in Fig. 2. We used a cosmic ray ionization rate of ζ_CR = 10⁻¹⁵ s⁻¹ (R/1 pc), which is consistent with the results of Yusef-Zadeh et al. (2013) and Harada et al. (2015). The radial dependence was introduced because the model with a constant cosmic ray ionization rate of ζ_CR = 2 × 10⁻¹⁵ s⁻¹ lead to a too extended HCO⁺ emission distribution. Thus, a centrally concentrated source of cosmic ray particles is assumed, most probably the colliding winds of the massive stars located within the inner half parsec of the Galaxy (Paumard et al. 2006). It was assumed that the relevant dust temperature for the molecular line emission is set by the heating through dust-gas collisions. These dust temperatures correspond to the cold dust component discovered by Etxaluze et al. (2011). In this way we neglected the molecular line emission from the photodissociation regions at the inner edge of the CND which is irradiated by the UV emission of early type stars located within the central cavity (Genzel et al. 2010). The radial distribution of the disk cloud mass, size, density, H₂ column density, and gas and dust temperatures are presented in Fig. 6. At the inner edge of the CND (R ∼ 1 pc) the characteristics of the clouds are M_cl ∼ 20 M_⊙, l_cl ≲ 0.1 pc, n_H2 ∼ 2 × 10⁶ cm⁻³, N_H2 ∼ 4 × 10²³ cm⁻², T_gas ∼ 300 K, and T_dust ∼ 30 K. The volume filling factor of the clouds is Φ_V = 0.035.

Fig. 6. Characteristics of individual model disk gas clouds across the CND at the time of interest. From upper left to lower right: cloud mass, size, density, H2 column density, gas and dust temperatures, optical depth, and brightness temperature.

The resulting radial profiles of the CO(2−1), CO(3−2), CO(6−5), HCN(3−2), HCO⁺(3−2), and HCO⁺(4−3) optical depths and brightness temperatures are shown in the two lower panels of Fig. 6. For R > 1 pc all lines are optically thick (τ_ν > 1). Whereas the optical depths of the CO and HCN lines decrease with increasing radius, the opposite trend is observed for the HCO⁺ lines. The HCO⁺ optical depth is small for R < 1.5 pc. This is due to the low HCO⁺ abundance, which is x_HCO+ ∼ 10⁻⁸ compared to x_HCN ∼ 10⁻⁶ at R = 1 pc. The HCN brightness temperature profile show a narrow peak around R ∼ 0.6 pc with a peak brightness temperature of about 250 K. The characteristics of these inner disk clouds are similar to those of hot cores. Their high HCN brightness temperatures are due to their high HCN abundances, up to x_HCN ∼ 10⁻⁶, as observed in Galactic hot cores (Boonman et al. 2001; Rolffs et al. 2011). The HCN brightness temperatures strongly decrease for radii larger than ∼1 pc. The HCO⁺ brightness temperatures decrease monotonically with increasing radius from ∼100 K in the center to ∼20 K at R = 4 pc. The CO brightness temperature has a peak of ∼280 K at ∼0.7 pc and decreases monotonically with increasing radius. This decrease is much shallower than that of the HCN lines.

As noted by Mills et al. (2013), the HCN line becomes a maser for temperatures higher than T ∼ 100 K and column densities in excess of N_HCN = 10¹⁶ cm⁻². These conditions are fulfilled by the disk clouds at R ≳ 2 pc. The reason for the high column density is the hot core chemistry which produces very high HCN abundances (x_HCN ∼ 10⁻⁷–10⁻⁶). To avoid HCN maser emission, our model clouds should have an about five times lower HCN abundance. This can be achieved by about 20% lower gas temperatures (T_gas ∼ 220 K instead of T_gas ∼ 260 K) or an about two times longer turbulent lifetime of the clouds. Most importantly, the resulting HCN brightness temperatures did not significantly change when these modifications were applied.

3.1. CO(6−5)

The comparison of the CO(6−5) moment 0 map with the observations of Requena-Torres et al. (2012) is presented in Fig. 7. The observed CO(6−5) map shows an inclined ring structure with two prominent lobes. The southern lobe has an about two times higher flux than the northern lobe.

Fig. 7. CO(6−5) emission distribution. Upper panel: model; lower panel: APEX/CHAMP+ observations integrated between −150 and 150 km s−1 (Requena-Torres et al. 2012). The spatial resolution is 9.5″. The flux is given in units of K km s−1, in Tmb.

The model CO(6−5) flux of the southern lobe is 1060 K km s⁻¹ compared to the observed value of 1780 K km s⁻¹, that of the northern lobe (from 60 to 150 km s⁻¹) is 490 K km s⁻¹ compared to the observed value of 790 K km s⁻¹. Thus, the flux ratio between the lobes is well reproduced by the model but the model fluxes are about 60% smaller than the observed fluxes. The inspection of the spectra of the northern and southern lobes (Fig. 8) shows that the observed linewidths are comparable to those of the model. However, the relative emission at low velocities (≲−80 km s⁻¹) is overestimated by the model by about a factor of two. The difference between the observed and model fluxes in the lobes is mainly due to a difference in the CO(6−5) brightness temperatures, which is linked to the disk cloud densities, sizes, and temperatures.

Fig. 8. CO spectra. Upper panel: CO(6−5) model spectra of the southern and northern lobes. Lower panel: CO spectra toward the southern (left) and northern (right) lobe of the CND (from Requena-Torres et al. 2012). All spectra are convolved to the same angular resolution of 22.5″. The flux density is given in units of K, in Tmb.

3.2. HCN(4−3) and CS(7−6)

Whereas the CO abundance in dense gas clouds is never significantly different from the canonical value of x_CO ∼ 10⁻⁴ (e.g., Lacy et al. 1994), the HCN and CS abundances depend on the gas chemistry. The comparison between of the model HCN(4−3) and CS(7−6) maps thus tells us if our model is not only able to reproduce the gas densities, surface densities, velocity dispersions, and temperatures but also the gas chemistry. The comparison of the HCN(4−3) moment 0 map with the observations of Montero-Castaño et al. (2009) is shown in Fig. 9. Nautilus yields HCN abundances of x_HCN ∼ 10⁻⁶ at a radius of 1 pc and x_HCN ∼ 10⁻⁷ at 2 pc. The HCN abundance is thus significantly enhanced compared to typical values in dense gas clouds (x_HCN ∼ 10⁻⁸; e.g., Hirota et al. 1998) due to the high density and temperatures that lead to hot core physics. The HCO⁺ abundance of several 10⁻⁸ is much lower than the HCN abundance in this region.

Fig. 9. HCN(4−3) integrated intensity map. Upper panel: model. The contour levels are in steps of 53 K km s−1 from 27 to 875 K km s−1 The spatial resolution is 4.6″ × 3.0″. Lower panel: observations from Montero-Castaño et al. (2009). The rms is 5.9 K km s−1. The contour levels are in steps of 35 K km s−1 from 18 to 550 K km s−1, but for the highest contour level at 574 K km s−1. Sgr A* is marked with a star.

The model gas ring is about 30% larger and somewhat more inclined than the observed CND. The bulk of the observed HCN(4−3) emission is confined to radii between 1 (25″) and 2 pc (50″) and the southern lobe, which shows an arc structure, is more prominent than the northern lobe. The maximum model HCN(4−3) integrated intensity is about 50% higher than the maximum observed integrated intensity. The observed clumpy structure of the CND is well reproduced by the model. Moreover, the observed arc structure of the southern lobe is much less prominent in the model. The observed faint HCN(4−3) emission east and west of the inner ring is not present in the model.

The comparison of the CS(7−6) moment 0 map with the observations of Montero-Castaño et al. (2009) is shown in Fig. 10. The CS abundance increases from x_CS = 6⁻¹⁰ at a radius of 1 pc to x_CS = 1.5⁻⁹ at 2 pc.

Fig. 10. CS(7−6) integrated intensity map. Upper panel: model. The contour levels are in steps of 22 K km s−1 from 11 to 275 K km s−1. The spatial resolution is 4.5″ × 3.1″. Lower panel: observations from Montero-Castaño et al. (2009). The rms is 2.5 K km s−1. The contour levels are in steps of 15 K km s−1 from 7 to 170 K km s−1, but for the highest contour level at 190 K km s−1. Sgr A* is marked with a star.

The maximum model CS(7−6) integrated intensity is about 35% higher than the maximum observed integrated intensity. The observed north-south asymmetry is not present in the model. The southern and northern extensions beyond offsets of ±50″ in the model moment 0 map are not observed. As for the HCN(4−3) moment 0 map, the observed clumpy structure of the CND is well reproduced by the model. Hsieh et al. (2021), who derived CS column densities from ALMA multi-transition CS observations, found N_CS = 10¹⁵–10¹⁶ cm⁻². Our model CS column densities of N_CS ∼ 1–3 × 10¹⁵ cm⁻² are situated at the lower end but are comparable to the observed distribution. Our assumed sulfur abundance of S/H = 10⁻⁶ is therefore justified.

3.3. HCN(3−2), HCN(4−3), HCO⁺(3−2), and HCO⁺(4−3)

Mills et al. (2013) observed the CND with APEX in the HCN(3−2), HCN(4−3), HCO⁺(3−2), and HCO⁺(4−3) lines (lower four panels of Fig. 11). They detected strong line emission from the southern lobe and rather faint line emission from the northern lobe. The observed mean HCN/HCO⁺ line ratio is about 1.7–1.8. The corresponding model maps are presented in the four upper panels of Fig. 11. The model HCN and HCO⁺ surface brightnesses of the two lobe are about a factor of two higher than the observed surface brightnesses. The model HCN/HCO⁺ line ratio is about 1.3–1.4 there, about 25% lower than the observed value.

Fig. 11. HCN(3−2), HCN(4−3), HCO+(3−2), and HCO+(4−3) moment 0 maps. Panels a–d: model. Contour levels are from 100 to 500 K km s−1, spaced by 50 K km s−1. Panels e–h: observations (Mills et al. 2013). The spatial resolutions are: HCN(3−2): 23.6″, HCN(4−3): 17.7″, HCO+(3−2): 23.5″, HCO+(4−3): 17.6″. Contours are linearly spaced. (e) HCN(3−2), contours from 75 to 265 K km s−1, spaced by 27 K km s−1. (f) HCO+(3−2), contours from 58 to 140 K km s−1, spaced by 17 K km s−1. (g) HCN(4−3), contours from 60 to 262 K km s−1, spaced by 29 K km s−1. (h) HC0+(4−3), contours from 58 to 140 K km s−1, spaced by 17 K km s−1.

Mills et al. (2013) observed four positions in the CND in the H¹³CN(3−2), H¹³CN(4−3), HC¹³O⁺(3−2), and H¹³CO⁺(4−3) lines. Within our model, we set the carbon isotope ratio [¹²C]/[¹³C] = 30, which is close the value of 25 adopted by Mills et al. (2013) to calculate the isotopologue emission. The resulting moment 0 model maps are presented in Fig. 12. The resulting model line ratios together with the observed line ratios are presented in Table 3. The model line ratios are about a factor of two higher than the observed line ratios. We conclude that our model reproduces the observed brightness temperatures and line ratios within a factor two. This represents an additional verification of our model gas chemistry.

Fig. 12. Model H13CN(3−2) (a), H13CN(4−3) (b), H13CO+(3−2) (c), and H13CO+(4−3) (d) moment 0 maps. The contour levels are those of Fig. 11 divided by a factor 4.

3.4. Discussion

Whereas the observed flux ratio between the northern and southern lobes is well reproduced in the CO line, it is significantly overestimated by the model in the HCN, HCO⁺, and CS lines. Thus, the conditions in terms of density, temperature, and molecule abundances of the gas in the northern lobe seem to be different from those of the southern lobe. Since we used a symmetric radial profile of the model brightness temperatures, we could not model such an effect. Our disk cloud densities are in broad agreement with those derived by Mills et al. (2013). These authors derived an about 10 times higher gas density from the HCN emission than from the HCO⁺ emission. We suggest that this difference at least partly stems from the fact that within the ∼1 pc APEX beam the HCN emission mainly comes from the inner edge of the CND, whereas the HCO⁺ emission has a significant contribution from gas at larger radii. The model disk cloud density at a radius of 1 pc is n_H2 ∼ 2 × 10⁶ cm⁻³, which is within the error bars of the density quoted by Mills et al. (2013). The model cloud density at a radius of 1.7 pc is about three times lower, n_H2 ∼ 7 × 10⁵ cm⁻³.

The somewhat smaller model HCN/HCO⁺ line ratio compared to the observations is most probably caused by the assumed cosmic ray ionization rate. The chosen radially dependent cosmic ray ionization rate represents the best compromise between the HCN surface brightness, which increases with increasing ζ_CR, and the HCN/HCO⁺ line ratio, which decreases with increasing ζ_CR. As an alternative we calculated a model with a constant cosmic ray ionization rate of ζ_CR = 2 × 10⁻¹⁵ s⁻¹. This modification did not significantly change the HCN emission but lead to an extended low surface density HCO⁺ envelope emitted by gas located at radii between 3 and 5 pc. Such an envelope is not observed.

We conclude that our forward modeling can reproduce the CO(6−5), HCN(3−2), H¹³CN(3−2), HCN(4−3), H¹³CN(4−3), HCO⁺(3−2), H¹³CO⁺(3−2), HCO⁺(4−3), and H¹³CO⁺(4−3) line emission distributions within a factor of two. Therefore, our model is able to reproduce in a satisfactory way the gas cloud sizes, densities, temperatures, velocity dispersions, and chemistry in the central 10 pc around the supermassive black hole of the Galaxy. Moreover, the assumed carbon isotope ratio of [¹²C]/[¹³C] = 30 is justified. The available molecular line observations of the CND can be reasonably reproduced by two massive infalling gas clouds on prograde orbits, which interact with a preexisting molecular gas disk.

3.5. The intercloud gas

The dense clouds within the CND might not be the only sources of molecular line emission. According to our analytical model (Sect. 2.2), the gas is of turbulent and clumpy nature. This means that turbulent eddies at scales between the driving length and the cloud size with different densities and temperature coexist. A consistent treatment of the scales as attempted by Vollmer et al. (2017) is beyond the scope of this article. For simplicity, we divided the ISM located in the CND into clouds and intercloud gas. The available interferometric molecular line observations of the CND indicate that the clumps are prominent and well-resolved. This implies that neither the brightness temperature of the intercloud gas dominates the overall emission nor the intercloud gas is optically thick. In the latter case, the intercloud gas, which has a high area filling factor, is expected to hide a significant fraction of the clouds via self-absorption in the molecular lines. The condition that the optical thickness of the intercloud medium should be below unity has important implications for its properties.

To calculate the mass fraction between of the dense clouds, we use the density probability distribution function of Padoan et al. (1997) for overdensities x:

(11)

where the standard deviation, σ, is given by

(12)

and ℳ = v_turb/c_s is the Mach number with the sound speed c_s. The mass fraction of gas with overdensities exceeding x is then

(13)

For ℳ = 20 and x = 30 we find . Thus, about 70% of the gas mass is expected to have lower gas densities than those of the dense clouds.

For the intercloud gas we assumed a size of l_driv = H and area and volume filling factors of unity. It turned out that the CO, HCN and HCO⁺ optical depths exceed unity for the turbulent heating presented in Eq. (6). The molecular lines emitted by the intercloud gas only become optically thin if the heating efficiency is decreased by a factor of 50. We therefore used

(14)

for the intercloud medium. The properties of these largest turbulent eddies are shown in Fig. 13. Compared to the gas clouds (Fig. 6) these eddies are larger, have an about 30 times lower density, an about 10 times lower column density, and an about two times lower gas temperature due to the low turbulent heating efficiency. Most importantly, the optical depth of these eddies is below unity for distances greater than ∼1 pc from the central black hole. Therefore, it cannot be detected in the CO, HCN, and HCO⁺ lines.

Fig. 13. Characteristics of individual model intercloud gas across the CND at the time of interest. From upper left to lower right: cloud mass, size, density, H2 column density, gas and dust temperatures, optical depth, and brightness temperature.

ISM turbulence is thought to be intermittent (see, e.g., Lazarian 2006): turbulence self-similarity is not exactly true even along the inertial range, that is at length scales greater than the dissipation length. Instead the turbulent fluctuations become increasingly sparse in time and space at smaller scales. Pan & Padoan (2009) found that in compressible 3D hydrodynamic simulations the turbulent dissipation is characterized by strongly intermittent fluctuations. A significant fraction of the kinetic energy is viscously dissipated in the finest, most intermittent structures. Dense filaments are sites of strong dissipation, but filaments of high dissipation rate are also found at the interface of low- and high-density regions. As a result, the dissipation rate and the gas density are practically uncorrelated in their simulations (Fig. 3 of Pan & Padoan 2009). Due to the intermittent fluctuations in the turbulent heating rate, a significant mass fraction of a molecular cloud is not heated by turbulence and the mass-averaged cloud temperature is decreased by a factor of two to three.

In the presence of magnetic fields ambipolar diffusion can represent an important contribution to the local intermittent turbulent heating if the ratio between the ambipolar diffusion and the turbulent crossing times is between about one and hundred (Li et al. 2012; Momferratos et al. 2014). The latter condition is approximately fulfilled only by the clouds. Stone et al. (1998) found that the dissipation time of MHD turbulence is on the order of the flow crossing time or smaller, even in the presence of strong magnetic fields. Weak magnetic fields are amplified and tangled by the turbulence, while strong fields remain well ordered. Moreover, these authors found that the density contrasts are larger for strong fields at fixed turbulent Mach number. This might indicate a lower turbulent dissipation rate at larger scales.

If we want the intercloud medium to be transparent for the molecular line emission, the area filling factor of optically thick portions has to be small. The intermittent heating of low-density gas certainly helps to decrease the area filling factor of optically thick gas at these densities. However, it might only be the fact that low-density, high-dissipation regions are found close to high density regions which ensures the low area filling factor of low-density, optically thick gas. Moreover, we can only speculate that a strong uniform magnetic field, as it is observed in the CND by Hsieh et al. (2018), helps to suppress the turbulent heating rate within the inertial range l_cl ≲ l ≲ H.

4. An NGC 1068-like model

Since we are able to reproduce the observed line emission distributions in terms of area filling factor, linewidth, and brightness temperature of the Circumnuclear Disk in the Galactic Center within a factor two, we are confident that the same model can be applied to the distribution of the molecular gas in the central 10 pc of NGC 1068, which is about 1800 times farther away than the Galactic Center. With a distance of 14.4 Mpc, 1 pc corresponds to 0.014″.

The underlying dynamical simulation is shown in Fig. 3. During the evolution of the simulation all gas cloud particles whose distances to the central black hole are smaller than 0.25 pc, which is close to the dust sublimation radius (GRAVITY Collaboration 2020), are removed from the simulation and counted as accreted. The so defined mass accretion rate is presented in Fig. 14. In addition, we show the mass accretion rate of a simulation where the infalling cloud has been removed, that is a simulation of an unperturbed accreting disk.

Fig. 14. Model mass accretion rate measured at R = 0.25 pc. Solid line: preferred model with ξ = 0.9, dashed line: preferred model without the massive infalling gas cloud, vertical line: time of interest.

In the case of an isolated gas ring the inner edge of the ring approaches the center and reaches a radius of 0.25 pc at t = 0.6 Myr. The mass accretion rate then increases to Ṁ = 0.07 M_⊙ yr⁻¹ at 80 Myr and stays constant to 95 Myr. In the case where the infalling cloud hits the gas disk, the inner edge of the ring reaches a radius of 0.25 pc at t = 0.2 Myr. The mass accretion rate rapidly increases to Ṁ = 0.5 M_⊙ yr⁻¹ within 0.2 Myr. It then decreases more slowly to Ṁ = 0.3 M_⊙ yr⁻¹ at t = 0.6 Myr and Ṁ = 0.2 M_⊙ yr⁻¹ at t = 0.9 Myr. At the time of interest (t = 0.47 Myr corresponding to the minimum χ²; see Sect. 2) the mass accretion rate is close to its maximum (Ṁ ∼ 0.5 M_⊙ yr⁻¹).

For the central region of NGC 1068 we assume γ = 300 cms, v_A = 1.5 km s⁻¹, and ζ_CR = 2 × 10⁻¹³ s⁻¹ (Table A.1). This ionization rate leads to molecular line emission distributions from the dense disk clouds, which are comparable to the available observations. We refer to Sect. 5.2 for a discussion on the different ionization fractions. In Sect. 2.5 we argue that the X-ray emission is entirely absorbed by dense gas located inside the molecular emitting region at (R ≲ 0.3 pc). Therefore, X-ray heating and ionization can be ignored. The molecular gas is heated by cosmic rays and turbulent mechanical energy injection. In all presented models the mechanical energy injection is by far the dominant heating mechanism. The effects of optical and UV emission on the gas are treated separately in Sect. 4.2.4. For all comparisons with observations the centroid coordinates of the central continuum source S1 are 02^h42^m40.709444″, −00 ° 00′47.9446″, J2000.0 (Gallimore et al. 2004).

The face-on view of the cloud distribution at the time of interest is shown in Fig. 3. A very dense gas disk with a radius of ∼1.5 pc might correspond to the maser disk observed in NGC 1068 (Greenhill et al. 1996; Gallimore et al. 1996). A rather symmetric distribution of clouds with a lower cloud density is observed at radii between 1.5 and 5 pc. As in Sect. 2.7 brightness temperature and a cloud surface area were assigned to each gas cloud particle of the dynamical simulations according to the assumed gas density, temperature, velocity dispersion, and size obtained from the analytical model. Strictly speaking, our analytical model, which assumes a symmetric disk with fully developed turbulence, is only applicable within this region. To take into account an intercloud medium (see Sect. 4.2.3) we assigned half of the total gas mass to the dense clouds and the other half to the intercloud medium.

4.1. Continuum emission and line absorption

We assume that the millimeter continuum emission is free-free emission produced by thermal electrons (Gallimore et al. 2004) whenever a gas cloud is directly illuminated by the AGN, and thermal dust emission otherwise. Within the dynamical models all gas clouds are optically thick with respect to the UV and optical emission of the AGN. We therefore identified the gas clouds in our simulation, which are not hidden by other clouds. In Vollmer et al. (2018) we suggested that there is a transition between the thick outer gas disk and the thin inner maser disk at a radius of ∼1.5 pc where a magnetocentrifugal wind is launched. We expect that the gas clouds become disrupted within this region leading to larger volume and area filling factors. To mimic this effect, we increased the cloud sizes for radii smaller than 2 pc to 0.01 pc for the calculation of the continuum emission. This also includes the location of the observed molecular maser disk at R ≲ 1.5 pc.

It is assumed that the ionization front has a constant surface density where the cloud becomes optically thick to the optical emission of the AGN. The free-free emission is proportional to the emission measure . We adopted the relation between the emission measure and the gas density n found by Vollmer & Duschl (2001a,b; panels a and c of their Fig. 5): EM ∝ n^4/3. A 3D datacube with the three position axes with the emission measure of all model clouds was created and a 2D map was created which was then convolved to the beamsize of the Impellizzeri et al. (2019) observations and normalized to their maximum flux density (see, e.g., lower left panel of Fig. 19). The 3D datacube was used to calculate molecular line absorption from clouds with optical depths higher than unity via

(15)

where v is the radial velocity. The excitation temperatures T_ex comes from our molecular line emission calculations (Sect. 4.2) and the background temperature T_bg from the normalized 3D datacube.

4.2. Molecular line emission

As for the CND, we first present the results for the distribution of the dense disk clouds (Sect. 4.2.1). The model moment 0, moment 1, pv diagrams, and spectra along the major axis are compared to the HCN(3−2), HCO⁺(3−2), HCO⁺(4−3), CO(2−1), CO(3−2), and CO(6−5) observations. We then discuss the influence of a putative intercloud medium (Sect. 4.2.3). In addition, we separately evaluated the influence of the photodissociation regions close to the central engine (Sect. 4.2.4).

HCN has a large dipole moment and therefore does not trace dense gas if there is another excitation mechanism that is faster than the H₂ collisions and independent of gas density. One such excitation path is through a vibrationally excited state, to which molecules can be pumped by infrared radiation (Carroll & Goldsmith 1981). The first vibrationally excited state of HCN is its bending state (v₂ = 1) 1024 K above the ground with an emitting wavelength of λ = 14 μm (Sakamoto et al. 2010). Following Sakamoto et al. (2010), we define an equivalent gas density

(16)

where T₀ = 1024 K corresponds to the energy gap between the two vibrational levels v = 0 and 1, A_vib = 3.7 s⁻¹ is the Einstein coefficient for the vibrational transition, and γ_{J, J − 1} is the collisional rate coefficient. Following Vollmer et al. (2017), HCN IR-pumping is implemented in the model by replacing the cloud density n_cl by n_equiv if n_equiv > n_cl in the HCN emission calculations.

4.2.1. The dense clouds

We used a cosmic ray ionization rate of ζ_CR = 2 × 10⁻¹³ s⁻¹, which lead to HCO⁺ brightness temperatures closest to the available observations. As stated in Sect. 2.5, we assumed that the X-ray emission is absorbed within the Compton-thick, fully ionized, and dust-free gas within the dust sublimation radius. The radial distribution of the disk cloud mass, size, density, H₂ column density, and gas and dust temperatures of the Q = 30 model are presented in Fig. 15. At a radius of R ∼ 2 pc the characteristics of the clouds are M_cl ∼ 1 M_⊙, l_cl ∼ 0.02 pc, n_H2 ∼ 10⁷ cm⁻³, N_H2 ∼ 3 × 10²³ cm⁻², and T_gas ∼ 600 K. We adopted the radial profile of the dust temperature of the disk midplane from the 3D radiative transfer calculations of Vollmer et al. (2018) for all models:

(17)

Fig. 15. NGC 1068 Q = 30 analytical model of a turbulent clumpy gas disk (Sect. 2.2). Characteristics of individual model disk gas clouds across the disk at the time of interest. From upper left to lower right: cloud mass, size, density, H2 column density, gas and dust temperatures, optical depth, and brightness temperature. IR pumping is included for the HCN and HCO+ lines.

The resulting radial profiles of the CO(2−1), CO(3−2), CO(6−5), HCN(3−2), HCO⁺(3−2), and HCO⁺(4−3) optical depths and brightness temperatures are shown in the two lower panels of Fig. 15. The optical depths of the CO lines decrease monotonically with radius. The CO(2−1) is optically thin for R > 1 pc. The HCN(3−2) line is optically thick everywhere. As for the CND in the Galactic Center, the small HCO⁺ optical depth for R ≲ 4 pc is due to the low HCO⁺ abundance: at R = 2 pc the model yields x_HCO+ ∼ 4 × 10⁻⁸ whereas x_HCN ∼ 2 × 10⁻⁶ due to the hot core physics.

All brightness temperature radial profiles show a central depression. Whereas the maximum of the HCN and CO lines is T_b ∼ 500–600 K, the maximum brightness temperature of the HCO⁺ lines is about ∼200 K. The maxima are reached at a radius of R = 1–3 pc.

4.2.2. A quantitative comparison between the model and observations

We anticipate our conclusion that the model can reproduce the available molecular line observations within a factor of two (Table 4).

The observed emission of all lines is reproduced by the model within 40%, except that of the CO(6−5) line, which is overestimated by about a factor of two by the model. The reason for this might be an overestimation of the gas densities or temperatures. We suggest that the small integrated HCO⁺(4−3) intensity indicates the presence of an intercloud medium (Sect. 4.2.3). In the following, all observed emission distributions were rotated counter-clockwise by 20°.

4.2.3. CO(2−1) and CO(3−2)

García-Burillo et al. (2019) observed NGC 1068 in the CO(2−1) and CO(3−2) line with a resolution of 40 mas (2.8 pc; Fig. 16). The CO(2−1) and CO(3−2) emission distributions are extended (0.12″ or 8.4 pc to the east and 0.1″ or 7.0 pc to the west). The pv diagram shows a velocity gradient in the outer disk with positive velocities to the east and negative velocities to the west. This gradient is more prominent in the CO(2−1) line than in the CO(3−2) line. A velocity gradient in the opposite direction is seen in the inner disk. This velocity gradient is more prominent in the CO(3−2) line than in the CO(2−1) line. The identification of this apparent counter-rotation signature in the major axis pv diagram is explained by García-Burillo et al. (2019) as due to the presence of outflowing gas inside the torus superposed to rotation. The highest brightness temperatures of both lines are found at velocities close to the systemic velocity.

Fig. 16. NGC 1068. Panel (a) observed CO(2−1) moment 0 map (I0 = 7923 K km s−1), (b) pv diagram (contours are Tb0 = 37.8 K) from García-Burillo et al. (2019), (c) dense disk cloud model CO(2−1) moment 0 maps, (d) pv diagrams along the major axis, (e) observed CO(3−2) moment 0 map (I0 = 7923 K km s−1), (f) pv diagram (Tb0 = 37.8 K) from García-Burillo et al. (2019), and (g) dense disk cloud model CO(3−2) moment 0 maps, (h) pv diagrams along the major axis. The model levels are the same as those of the observed maps.

The corresponding model CO(2−1) and CO(3−2) moment 0 maps and pv diagrams are presented in Fig. 16. The extent of the model CO(2−1) and CO(3−2) emission distributions is well comparable to those of the observations. The linewidths of the observed CO(2−1) emission of the pv diagram are somewhat smaller but comparable to those of the model. Especially the observed western extension at negative velocities is present in the model CO(2−1) diagram. As observed, the model CO(3−2) emission distribution is more centrally concentrated than the model CO(2−1) emission distribution. The observed CO(3−2) linewidths are about 30% smaller than the model linewidths. Moreover, the western extension at negative velocities present in the model CO(3−2) pv diagram is absent in the observed CO(3−2) data.

4.2.4. CO(6−5)

Gallimore et al. (2016) observed the center of NGC 1068 in the CO(6−5) line with a resolution of 80 mas (5.6 pc). The moment 0 map and the pv diagram are presented in the upper panels of Fig. 17. The pv diagram revealed CO(6−5) emission at very high velocities close to the central black hole (v > 400 km s⁻¹). The highest CO(6−5) brightness temperatures are found in the center at the systemic velocity. The corresponding model CO(6−5) moment 0 map and pv diagram are shown in Fig. 17. The extent of the observed CO(6−5) distribution and the linewidth of the high brightness temperature emission (pv diagram) are reproduced by the model. As for the CO(3−2) emission, the western extension at negative velocities present in the model CO(6−5) pv diagram is absent in the observed CO(6−5) diagram. The observed high-velocity emission is absent in the model pv diagram.

Fig. 17. NGC 1068. Panel (a) observed CO(6−5) moment 0 map (Gallimore et al. 2016) (I0 = 6657 K km s−1), (b) CO(6−5) pv diagram along the major axis (Tb0 = 38.0 K), (c) dense disk cloud model CO(6−5) moment 0 maps, and (d) pv diagrams along the major axis. The model levels are the same as those of the observed maps.

4.2.5. HCN(3−2)

Imanishi et al. (2018, 2020) and Impellizzeri et al. (2019) observed NGC 1068 in the HCN(3−2) line with ALMA. Whereas the Imanishi et al. (2018) data have a spatial resolution of 40–70 mas (2.8–4.9 pc), the data of Impellizzeri et al. (2019) and Imanishi et al. (2020) have a resolution of 20 mas (∼1.4 pc). The comparison between the model and the high resolution HCN(3−2) datacube is presented in Fig. 18 together with the observations by Impellizzeri et al. (2019). The main characteristics of the observations are: (i) HCN(3−2) emission is detected in absorption against the continuum in the central resolution element, (ii) the bulk of the emission comes from a ≲3 pc region, with (iii) the western side being significantly brighter than the eastern side, (iv) emission at radii ≳3 pc is blueshifted in the west and redshifted in the east. On the other hand, emission at radii ≲2 pc is redshifted in the west and blueshifted in the east and coincides with the maser disk. The dense molecular gas in the gas disk appears counter-rotating between the innermost and outer parts, and (iv) a redshifted high-velocity (170–370 km s⁻¹) emission component is present at the innermost western side of the gas disk.

Fig. 18. NGC 1068. Left half: high resolution HCN(3−2) observations (Impellizzeri et al. 2019). Panel (a) moment 0 map (I0 = 17300 K km s−1), (b) pv diagram along the major axis, (c) pv diagram along the minor axis (Tb0 = 70 K), (d) moment 1 map (in km s−1), (e) continuum map (contour levels are (0.019, 0.037, 0.075, 0.15, 0.3, 0.6, 1, 2, 4, 8)×43 K). Right half: high resolution HCN(3−2) emission of the dense disk cloud model with IR pumping. Panel (f) moment 0 map, (g) pv diagram along the major axis, (h) pv diagram along the minor axis, (i) unsmoothed pv diagram along the major axis, (j) moment 1 map (in km s−1), (k) continuum map. The model levels are the same as those of the observed maps. The additional white contour is at 80 K.

The observed east-west asymmetry of the moment 0 map is present but less prominent in the model. Qualitatively, points (i), (ii) and (iv) are reproduced by the model. The observed western high-velocity emission is partly seen in absorption in the model. In particular, the observed counter-rotation between the inner and outer gas disk in the model pv diagram and the velocity field (moment 1 map) is present in the model. The model and observed linewidths are comparable.

On the other hand, the model has the following shortcomings: The model moment 0 map does not show a central absorption. Whereas the maxima of the emission lines are slightly negative and close to the systemic velocity in the observations, they are located at ∼ ± 70 km s⁻¹ in the model disk. The model pv diagram along the minor axis shows HCN in emission on the eastern side, whereas it is observed on the western side.

4.2.6. HCO⁺(3−2)

Imanishi et al. (2018, 2020) observed NGC 1068 in the HCO⁺(3−2) line with ALMA at a resolution of 50–70 mas (3.5–4.9 pc) and 20 mas (1.4 pc), respectively. The comparison between the model high resolution HCO⁺(3−2) datacube is presented in Fig. 19 together with the observations by Imanishi et al. (2020). The observed counter-rotation between the inner and outer parts of the gas disk is also present in the model HCO⁺(3−2) pv diagram and velocity field (moment 1 map). In contrast to the observations, the model HCO⁺(3−2) pv diagram shows a significantly deeper absorption line in the central resolution element than the model HCN(3−2) pv diagram. As for the HCN(3−2) line, the HCO⁺ emission of the model pv diagram along the minor axis is located on the opposite side of the observed HCO⁺ emission.

Fig. 19. NGC 1068. Left half: high resolution HCO+(3−2) observations (Imanishi et al. 2020). Panel (a) moment 0 map (I0 = 18420 K km s−1), (b) moment 1 map (in km s−1), (c) pv diagram along the major axis, and (d) pv diagram along the minor axis (Tb0 = 75 K). Right half: dense disk cloud model. Panel (e) moment 0 map, (f) moment 1 map (in km s−1), (g) pv diagram along the major axis, and (h) pv diagram along the minor axis. The model levels are the same as those of the observed maps.

4.2.7. HCO⁺(4−3)

García-Burillo et al. (2019) observed NGC 1068 in the HCO⁺(4−3) line with a resolution of 40 mas or 2.8 pc (left side of Fig. 20). The emission distribution is concentrated within the inner ±50 mas or 3.5 pc. The corresponding model HCO⁺(4−3) moment 0 map and pv diagram are presented in Fig. 20. The observed HCO⁺(4−3) emission is about three times stronger than the model emission. Contrary to the observations, most of the model HCO⁺(4−3) line is seen in absorption at negative velocities. We suggest that the small integrated HCO⁺(4−3) intensity indicates the presence of an intercloud medium (Sect. 4.2.3).

Fig. 20. NGC 1068. Panel (a) observed HCO+(4−3) moment 0 map (I0 = 7217 K km s−1), (b) pv diagram (Tb0 = 47.5 K) from García-Burillo et al. (2019), (c) dense disk cloud model HCO+(4−3) moment 0 maps, and (d) pv diagrams along the major axis. The model levels are the same as those of the observed maps.

4.2.8. HCN, HCO⁺, H¹³CN, and H¹³CO⁺ spectra along the major axis

To make a more quantitative comparison between the model and the observations, we present the observed HCN(3−2) spectra (Impellizzeri et al. 2019) and HCO⁺(3−2), H¹³CN(3−2), and H¹⁴CO⁺(3−2) spectra (Imanishi et al. 2020) from resolution elements along the major axis in Appendix C. The model resembles the observed HCN(3−2) and HCO⁺(3−2) spectra only qualitatively. Given the tiny amount of H¹³CN(3−2) emission and the absence of H¹³CO⁺(3−2) emission in the model spectra, our assumed isotope ratio ¹²C/¹³C = 30 is justified.

4.2.9. Infrared pumping

The results for the HCN(3−2) emission in the absence of IR pumping (see Sect. 4.2) of the Q = 30 model are presented in Fig. 21.

Fig. 21. IR pumping. Left half: high resolution HCN(3−2) emission of the dense disk cloud model. Panel (a) moment 0 map, (b) moment 1 map, (c) pv diagram along the major axis, and (d) pv diagram along the minor axis. The model levels are the same as those of the observed maps (Fig. 18). Right half: model high-resolution HCN(3−2) spectra along the major axis.

The influence of IR-pumping is best appreciated in the pv diagrams and spectra: the brightness temperatures decrease by about a factor of two and thus the area of regions with detected HCN(3−2) emission decreases. The total integrated intensity is 1.2 Jy km s⁻¹ instead of 2.7 Jy km s⁻¹ with IR pumping. We conclude that IR-pumping increases the HCN(3−2) emission by about a factor of two. This is consistent with the small amount of HCN-VIB J = 3-2 emission detected by Imanishi et al. (2020).

4.2.10. The intercloud gas

As for the CND (Sect. 3.2), the contribution of the intercloud gas to the molecular line emission has to be evaluated. Assuming a lognormal probability density function for the gas density (Eq. (13)) with ℳ = 35 and x = 50 we find . Thus, as in the CND about 70% of the gas mass is expected to have lower gas densities than that of the dense clouds. Since the mean gas density is proportional to Q⁻¹, the mean gas density of the thick gas disk in NGC 1068 is expected to be about 10 times higher than in the CND. The probability that the intercloud gas in the thick disk of NGC 1068 is optically thick leading to important self-absorption in the molecular lines is thus higher than in the CND.

The NIR interferometric observations of the GRAVITY Collaboration (2019) revealed an extinction of A_V ∼ 100 in front of the elongated central emission corresponding to a gas column density of ∼10²³ cm⁻². The same column density was measured toward the position of the AGN from the CO(2−1) and CO(3−2) intensities by García-Burillo et al. (2019). This gas has to have a high area filling factor because a partial coverage of the extended NIR emission by optically thick clouds would not lead to the observed reddening. It is remarkable that the model of continuous thick stratified gas disk of Vollmer et al. (2018) reproduces the observationally derived column density in front of the central source. Based on the available observations it cannot be excluded that the intercloud gas dominates the molecular line emission. In order to obtain an optically thin intercloud gas, the Toomre Q parameter has to be set to Q = 100 and the dissipation rate at the driving length must be decreased by at least a factor of 1/200 compared to Eq. (6) (left panel of Fig. 22). This is four times less than the factor required for the CND. As for the CND, we can only speculate that a strong uniform magnetic field, as it is observed by Lopez-Rodriguez et al. (2020), helps to suppress the turbulent heating rate within the inertial range l_cl ≲ l ≲ H. Indeed, these authors determined the ratio between the turbulent and the uniform magnetic field strengths to be 0.06.

Fig. 22. NGC 1068 intercloud gas model. Characteristics of individual model intercloud gas across the disk at the time of interest. Q = 30, ζCR = 2 × 10−13 s−1, and full turbulent heating. Gas cloud characteristics. From upper left to lower right: cloud mass, size, density, H2 column density, gas and dust temperatures, optical depth, and brightness temperature.

On the other hand, for an optically thick intercloud gas with brightness temperatures high enough to be detected by the available observation setups, we had to assume a Q = 30 disk with full heating (Eq. (6)) and diffuse clouds of size H/2 and density of three times the mean density. As for the dense clouds we assumed a CR ionization rate of ζ_CR = 2 × 10⁻¹³ s⁻¹. The resulting properties of the diffuse clouds are presented in the right panel of Fig. 22. For radii ≲2 pc all lines are optically thick. The maximum of the molecular brightness temperatures is located at a radius of ∼0.5 pc.

The resulting high-resolution model HCN(3−2) and HCO⁺(3−2) observations are presented in Fig. 23. The HCN(3−2) emission of the intercloud gas is comparable to that of the dense disk clouds (Fig. 18). As expected from the surface brightness profiles, the line emission is concentrated within the inner ∼3 pc. The extent of the model HCN intercloud emission is thus smaller than that of the observations (lower panels of Fig. 18). The model HCN linewidth is smaller but comparable to that of the observations. The highest observed velocities (v > 200 km s⁻¹) are not reproduced by the model. The intercloud model does not show the HCN(3−2) line in absorption.

Fig. 23. NGC 1068 intercloud gas model. Left half: HCN(3−2); right half: HCO+(3−2). Panels (a) and (e) model moment 0 maps, (b) and (f) model moment 1 maps, (c) and (g) pv diagrams along the major axis, (d) and (h) pv diagrams along the minor axis. The model levels are the same as those of the observed maps (Figs. 18 and 19), except for the white contour at 80 K.

As the model HCN(3−2) emission, the model HCO⁺(3−2) emission is more centrally concentrated. Otherwise, it shows similar spatial and kinematic structures with somewhat higher brightness temperatures.

The model HCO⁺(4−3), CO(2−1), CO(3−2), and CO(6−5) are presented in Fig. 24. The corresponding observations are shown in Figs. 20, 16, and 17. The line emission is concentrated within a radius of ∼3 pc for all lines. In line with the observations, the HCO⁺(4−3), CO(2−1), CO(3−2), and CO(6−5) lines are not seen in absorption.

Fig. 24. NGC 1068 intercloud gas model. Panels (a) and (b) HCO+(4−3) moment 0 maps and pv diagrams, (c) and (d) CO(2−1) moment 0 maps and pv diagrams, (e) and (f) CO(3−2) moment 0 maps and pv diagrams, (g) and (h) CO(6−5) moment 0 maps and pv diagrams. The model levels are the same as those of the observed maps (Figs. 16, 17, and 20).

We conclude that the contribution of the intercloud gas to the total emission is expected to be highest at radii ≲3 pc. The observed absence of line absorption in the HCO⁺(4−3) lines might be a hint to the existence of an intercloud medium.

4.2.11. The photodissociation region

In Sect. 2.5 we argue that the X-ray emission is entirely absorbed by dense gas located inside the region where the gas emits molecular line emission at R ≲ 0.3 pc. This gas is most-probably located within the dust sublimation radius and is therefore dust-free and optically thin for the UV and optical emission. All gas surfaces, which are directly illuminated by UV and optical emission from the AGN, give rise to photodissociation regions (PDR). The models presented in Sects. 4.2.1 and 4.2.3 do not include the effects of optical and UV emission on the gas. To determine the influence of PDRs on the observed molecular emission distributions, we treat the PDR emission separately and assume in the following that the gas is only heated by the UV and optical emission, meaning that there is no turbulent mechanical heating. These PDR strongly emit in molecular lines (e.g., Meijerink & Spaans 2005; Meijerink et al. 2007). We assumed a bolometric luminosity of L ∼ 3 × 10⁴⁴ erg s⁻¹ for NGC 1068 (Vollmer et al. 2018). The HCN, HCO⁺, and CO line emission were calculated with the Meudon PDR code (Le Petit et al. 2006; Le Bourlot et al. 2012). The code considers a stationary plane-parallel slab of gas and dust illuminated by a radiation field. It solves at each point in the cloud, the radiative transfer in the UV taking into account the absorption in the continuum by dust and in discrete transitions of H and H2. The model computes the thermal balance taking into account heating processes such as the photoelectric effect on dust, chemistry, cosmic rays, etc. and cooling resulting from infrared and millimeter emission of the abundant species. Chemistry is solved for any number of species and reactions. Once abundances of atoms and molecules and level excitation of the most important species have been computed at each position in the cloud, line intensities are deduced by a post-processor code. We assumed A_V = 30 for the plane-parallel slab and a cosmic ray ionization rate of ζ_CR = 10⁻¹⁴ s⁻¹. Such a high rate was found by Le Petit et al. (2016) who compared their PDR models with emission from dense clouds in the central molecular zone (CMZ) of the Galaxy.

The radial profiles of the molecular line emission for the dense clouds in Q = 30 gas disk are presented in Fig. 25. The CO and HCO⁺(4−3) line emission from the PDR is stronger than, the PDR HCO⁺(3−2) line emission is comparable to that due to mechanical and CR heating for all observed transitions. On the other hand, the HCN(3−2) line emission from the PDR is much weaker than that due to mechanical and CR heating heating. The reason for the low HCN(3−2) line emission from the PDR is caused by the temperature gradient of the PDR: the temperature of the emitting region with an optical depth of τ_HCN ∼ 1 is much lower in the PDR model than in the model with turbulent mechanical heating where all gas is uniformly heated. The high gas temperature at high densities leads to a strong increase of the HCN abundance via a hot core chemistry.

Fig. 25. Radial profiles of the molecular line emission from the dense clouds (thin lines) and the associated PDR regions (thick lines). The CR ionization rate is ζCR = 10−14 s−1.

The 3D configuration of the clouds is a thick ring with an inner wall which is directly illuminated (Vollmer et al. 2018). Since the observed brightness temperature of the directly illuminated clouds is proportional to the cloud brightness temperature and their area filling factor, we expect the contribution of the PDR to the total emission to be dominant close to the edge of the thick ring at a distance of about 1.5 pc from the central black hole (Vollmer et al. 2018). The mean free path of the clouds at the inner edge is about 1 pc corresponding to the thickness of the PDR-dominated region.

As before, we assume that the gas in the ring is turbulent and clumpy. The dense clouds are thus embedded in an intercloud medium. If the intercloud gas is present, it is more likely that it is directly illuminated than the embedded clouds. We therefore calculated the PDR molecular line emission for the intercloud medium discussed in Sect. 4.2.3 (Fig. 26). As for the dense clouds, the PDR CO line emission is stronger, the PDR HCN emission is weaker than that of the mechanically and CR heated intercloud medium. For R < 5 pc the emission of both PDR HCO⁺ lines is weaker than, for R ≥ 5 pc it is comparable to and weaker than that of the intercloud medium.

Fig. 26. Radial profiles of the molecular line emission from the intercloud gas (thin lines) and the associated PDR regions (thick lines). The CR ionization rate is ζCR = 10−14 s−1.

We conclude that the central AGN most probably directly illuminates the inner edge of the thick gas ring. Based on our model calculations we do not expect that the emission from the PDR regions dominates the total molecular line emission at the inner edge of the ring but it can increase the emission by about a factor of two.

4.2.12. HCN(4−3), CN(3−2), and CS(7−6)

Unpublished ALMA Cyle 2 Band 7 HCN(4−3), CN(3−2), and CS(7−6) observations are compared to our models. Details of the observations are given in Appendix D.

4.2.13. Model and observed central spectra

Since the HCN(4−3), CN(3−2), and CS(7−6) emission distributions are barely resolved we only focus on the spectrum of the central pixel (S1). The observed and model integrated spectra for the dense cloud and diffuse gas models are presented in Fig. 27.

Fig. 27. ALMA integrated spectra of the central source S1 of NGC 1068 (black line). Upper panel: HCN(4−3); middle panel: CN(3−2); lower panel: CS(7−6). The gray line corresponds to the second component of the CN(3−2) triplet. Colored lines: dashed: S1 spectrum from the dense disk cloud model; dotted: S1 spectrum from the diffuse gas model.

The height and the width of the ALMA HCN(4−3) spectrum is well reproduced by the dense disk cloud and intercloud models with IR pumping, as it is expected from Sect. 4.2.1. The integrated intensities of the models without IR pumping are about two times for the dense cloud and four times for the diffuse gas smaller than the corresponding models with IR pumping.

For the CN molecule, we focussed our attention to the CN(N = 3–2, J = 7/2–5/2) transition at 340.248 GHz, which is strongest line of the CN(3−2) triplet. The models reproduce the observed linewidth. The model CN emission of the dense disk clouds is about nine times smaller than the observed CN emission. The CN emission of the diffuse gas model is still smaller. CN has a high critical density (about half of the HCN critical density) and is primarily formed from photodissociation of HCN and neutral-neutral reactions with N, C₂, CH₂, and CH (Aalto et al. 2002; Boger & Sternberg 2005; Chapillon et al. 2012). Intermediate stages in the reaction pathways involve neutral and ionized carbon (Boger & Sternberg 2005). CN is thus thought to preferentially form in regions illuminated by intense radiation fields (PDRs), whereas it is destroyed in dense gas via CN+H₂ → HCN+H. It is thus expected that the PDRs (Sect. 4.2.4) should thus significantly contribute to the CN emission.

The CS(7−6) emission of the dense disk cloud model is about a factor of three smaller than that of the intercloud model is well comparable to the observed CS(7−6) emission. Given the uncertainties of the model, we conclude that the model is able to reproduce the observed CS(7−6) spectrum with the same sulfur abundance as for the CND in the Galactic Center (S/H = 10⁻⁶). This abundance is intermediate between that observed in cold dense gas (S/H ∼ 10⁻⁷; e.g., Quan et al. 2008) and the solar sulfur abundance of S/H = 1.5 × 10⁻⁵ (Jenkins 2009) used in the hot core models of Vidal & Wakelam (2018). It is consistent with the sulfur abundance found in the hot core in the star-forming region Orion KL (Esplugues et al. 2014). We can only speculate that the thermal evaporation of icy grain mantles at dust temperatures > 100 K might be responsible for the enhanced sulfur abundance with respect to cold cores. In addition, sputtering induced by collisions may be efficient in this turbulent environment to erode the grain mantles and release the sulfuretted species to the gas phase as proposed by Fuente et al. (2016) for the Galactic dark cloud Barnard 1.

5. Discussion

The comparison between the moment 0 and moment 1 maps and the pv-diagrams along the major and minor axes of the different models and the HCN(3−2), HCO⁺(3−2), HCO⁺(4−3), CO(2−1), CO(3−2), and CO(6−5) line observations of NGC 1068 showed that the model reproduces the available observations within a factor of about two. In the following, we discuss the stability of counter-rotating disks, the CR ionization fractions, elucidate the role of the continuum emission, and discuss the role of outflows or winds. Furthermore, the implications of the proposed models are discussed in terms of a twofold outflow and AGN variability.

5.1. On the stability of counter-rotating disks

Misalignments between the ionization cones of AGN, and hence the obscuring tori, and the disks of the host galaxies are frequently observed (e.g., Fischer et al. 2013). Hopkins et al. (2012) stated that these misalignments occur for at least two reasons: first, large-scale gas fragmentation can occur (part of a spiral arm or other instability fragments and sinks to the center), which can dramatically change the nuclear gas angular momentum content. Secondly, even in perfectly smooth flows, secondary bars in the presence of dissipative gas processes will tend to decouple their angular momentum from the primary bar. Inflow and dissipation lead to a decoupling the inner mode angular momentum and orbit plane from that of the outer mode. Angular momentum exchange in the gas in the central regions of galaxies can be strongly dominated by supersonic gas shocks surrounding strong torquing regions in the stellar nuclear disk with lopsided/eccentric (m = 1) modes (Hopkins & Quataert 2011; Bacon et al. 2001; Jacobs & Sellwood 2001; Salow & Statler 2001; Sambhus & Sridhar 2002). These modes can resonantly exchange angular momentum with the pattern at larger radii in the manner of nested bars, leading (in plane) to possible reversals and counter-rotation of the pattern, which in turn reverses the sense of torques on the gas. An infall of counter-rotating gas is thus expected, but should be rare.

Quach et al. (2015) and Dyda et al. (2015) studied counter-rotating gas disks analytically and through numerical simulations. They found that strong supersonic Kelvin-Helmholtz instabilities develop at the interface between the two disks. The growth rates are on the order of or larger than the angular rotation rate at the interface. Dyda et al. (2015) showed that the onset of the instability leads to the creation of a gap between the two disks (their Fig.11). The size of the gap is about the vertical size of the outer disk. The gap is created by the mixing of the two counter-rotating components, strong heating of the mixed gas, vertical expansion of the gas, and subsequent annihilation of the angular momenta of the two components. As a result, the heated gas free-falls toward the center over the surface of the inner disk. Thereafter, the gap begins to fill as matter from the outer disk moves inwards due to its viscosity. In addition, gas from the inner disk will also fill the gap from the other side due to viscous diffusion. After a time of t ∼ h/v_r the instability develops again. Our Q = 30 disk model of NGC 1068 has a height of about 0.5 pc at R = 1.5 pc and a radial velocity of v_r = 12 km s⁻¹. The width of the gap is thus about 0.5–1 pc and the time to fill it is t = 4 − 6 × 10⁴ yr corresponding to 7–11 dynamical timescales. This is quite short, but certainly long enough to be observable. By construction, in our simple sticky particle simulations instabilities cannot develop. The expected gap between the counter-rotating disks is thus not present in the simulations. However, the mass accretion rate of the outer thick disk giving rise to v_r ∼ 10–15 km s⁻¹ is present. It is this radial accretion which sets the lifetime of the system.

5.2. The cosmic ray ionization and gas heating

The CR ionization fraction is mainly constrained by the HCO⁺ emission. For the models where only dense disk clouds are present in the thick gas disk (Sect. 4.2.1) a CR ionization fraction of ζ_CR ∼ 10⁻¹³ s⁻¹ yields the best results. For a significant contribution to the molecular line emission from the intercloud medium, the same CR ionization rate of is required (Sect. 4.2.3). The CR ionization rate within the thick gas disk of NGC 1068 is thus 100 times higher than that in the CND (ζ_CR ∼ 10⁻¹⁵ s⁻¹; Sect. 3). Within the inner 8 pc of NGC 1068 the gas heating rate due to cosmic rays is less then 10% of the turbulent heating rate for the dense gas clouds and less than 20% for the intercloud gas. We thus conclude that within the framework of the model turbulent mechanical heating dominates the heating in the thick gas disk of NGC 1068. In the present model turbulence is maintained by the gain of gravitational potential of the accreting gas.

5.3. The role of the continuum emission

Our modeling allows us to investigate the role of the continuum emission in a simple way: we can remove the continuum emission, in other words we can switch off the central AGN. The resulting moment 0 and moment 1 maps together with the pv diagrams are presented in Fig. 28 for the Q = 30 model with and without IR-pumping.

Fig. 28. NGC 1068 dense disk cloud model without absorption. Left panel: Q = 30 model high resolution HCN(3−2) emission of the dense disk cloud model. Right panel: same as left panel but with IR-pumping included in the model. Panels (a) and (e) moment 0 maps, (b) and (f) moment 1 maps, (c) and (g) pv diagrams along the major axis, (d) and (h) pv diagrams along the minor axis. The model levels are the same as those of the observed maps (Fig. 18). The additional white contours levels are 80 and 90 K.

By comparing Fig. 28 to Fig. 18 it becomes clear that only the central resolution element is affected by absorption, as expected. HCN(3−2) emission is detected within a velocity range between −160 km s⁻¹ and 180 km s⁻¹. Whereas the velocity range at negative velocities is compatible with the observed velocity range, HCN(3−2) emission is observed up to a velocity of ∼300 km s⁻¹ (panel b of Fig. 18). HCN(3−2) emission at these high velocities is not present in the model.

We interpret this finding in the following way: we adopt the framework of the Vollmer et al. (2018) model, where the continuum emission is produced by the thin maser disk and the inner edge of the thick gas disk. Whenever a gas cloud of the thick disk is located in front of the continuum emission from the point of view of the observer, HCN(3−2) is seen in absorption. We suggest that the observed HCN(3−2) emission at velocities > 200 km s⁻¹ is emitted by the thin maser disk located within a distance of ∼1.5 pc around the central black hole. For a black hole mass of 10⁷ M_⊙ the Keplerian velocity is about 210 km s⁻¹ at 1 pc and 190 km s⁻¹ at 1.25 pc. This corresponds to a channel width of 20 km s⁻¹. The size of the emitting region is thus about 0.2 pc whereas the observational beam is 1.4 pc. The area filling factor is thus about (0.25/1.4)² = 0.03. The observed flux density is about 0.7 mJy beam⁻¹ = 30 K. Thus, the brightness temperature of the emitting gas is about 30/0.03 K = 1000 K. We think that such a high brightness temperature can only be reached if we observe at least some HCN maser emission from this region. The regions of H₂O and HCN maser emission are thus coincident. In addition, Gallimore et al. (2004; their Fig. 7) showed that a part of the maser spots or maser disk extend beyond the VLBA 5 GHz continuum emission distribution in the northwest of the central black hole. The velocities of these maser spots are indeed higher than 200 km s⁻¹. A significant part of the HCN emission might thus not be absorbed by the central continuum emission.

5.4. An outer and inner outflow

A wind or outflow is observed in NGC 1068 at different scales and in different gas phases: The line emission of ionized gas at scales of 10–100 pc has an hourglass structure and its kinematics are consistent with an outflow within a hollow-cone structure (May & Steiner 2017; Das et al. 2006; Miyauchi & Kishimoto 2020). At smaller scales, the morphology of the point-source subtracted ALMA band 6 continuum emission (Impellizzeri et al. 2019) and the polar dust emission from MIR interferometric observations (López-Gonzaga et al. 2014) are also consistent with an outflow. Gallimore et al. (2016), García-Burillo et al. (2019), and Impellizzeri et al. (2019) interpreted the high radial velocities and large linewidths of the molecular line emission along the minor axis of the molecular torus as a molecular outflow in the inner region of the massive gas disk (R ≲ 3 pc). Furthermore, the torus is connected to the 200 pc-size gas ring through a network of gas lanes whose kinematics are accounted for by a 3D outflow geometry (García-Burillo et al. 2019). We call the latter molecular and dusty wind observed at scales between a few parsecs and ∼100 pc the outer wind.

Such an outflow, which is definitely present in NGC 1068, is not included the present model. The results presented in Sect. 4.2 suggest that a good part of the kinematical features of the molecular line emission, which were interpreted as an outflow, might be due to the perturbed kinematics caused by the gas-gas collision of two counter-rotating entities. On the other hand, the present model cannot account for all observed emission: the large linewidth and velocity-asymmetry of the HCN(3−2) and HCO⁺(3−2) lines observed in emission in the pv diagram along the minor axis of the molecular torus at R ≲ 2 pc (Figs. 18 and 19) and the absorption wing at negative velocities in the HCN(3−2) spectrum of the central position (Fig. C.1). At least these two features are most probably caused by the molecular outflow. It is not excluded that the signatures of the gas-gas collision and the outflow are superimposed in the inner R ≲ 3 pc of NGC 1068.

The mass inflow rate from the dynamical model at 0.25 pc (Sect. 4), although rather variable, is around 0.5 M_⊙ yr⁻¹ at the selected timestep. This matches the Vollmer et al. (2018) model quite well for the inflow rate of 0.7 M_⊙ yr⁻¹ in the thin disk. Most of the inflowing gas must be blown out rather than accreted onto the central black hole. Within the scenario of Vollmer et al. (2018) a magnetocentrifugal (outer) wind aided by radiation pressure is launched from the thick gas disk at a radius of R ∼ 1.5 pc. This wind is naturally dusty and is assumed to be responsible for the observed polar infrared emission at the ten parsec scale (López-Gonzaga et al. 2014). The model of Vollmer et al. (2018) yielded a mass outflow rate of Ṁ_wind ∼ 0.9 M_⊙ yr⁻¹ of the magnetocentrifugal wind compared to a mass accretion rate of Ṁ_thick disk ∼ 1.6 M_⊙ yr⁻¹ of the thick gas disk. The thin gas disk has still a high mass accretion rate³Ṁ_thin disk ∼ 0.7 M_⊙ yr⁻¹. This has to be compared to the estimate of the mass accretion rate onto the central black hole of Ṁ_BH ∼ L_bol/(0.1 c²) ∼ 0.05 M_⊙ yr⁻¹ given the bolometric luminosity of L ∼ 3 × 10⁴⁴ erg s⁻¹ (Vollmer et al. 2018). It becomes clear that the radial gas inflow within the thin maser disk cannot and is not maintained at radii ≲0.1 pc where the optically emitting hot accretion disk is located. A natural way to explain the ratio Ṁ_BH/Ṁ_thin disk ∼ 0.1 is to assume a second mass outflow inside the dust sublimation radius of 0.25 pc (GRAVITY Collaboration 2020). This massive second outflow with a mass outflow rate of Ṁ_outflow ∼ 0.6 M_⊙ yr⁻¹ is naturally consistent with the suggestion that the Compton-thick gas layer is located within the dust sublimation radius (Sect. 2.5).

5.5. Sub-Keplerian rotation

It is quite surprising that the observed radial velocities of the HCN, HCO⁺, and CO lines do not reach the values expected from the rotation curves at radii R ≳ 2 pc, especially at positive velocities (García-Burillo et al. 2019; Imanishi et al. 2020, see also panels b and f of Fig. 16 and panel b of Fig. 17). There are three possible explanations for this finding: (i) a relatively low inclination angle of the outer gas disk, (ii) noncircular motions of or within the outer disk, (iii) a significant radial pressure support, and (iv) different emission conditions in the regions of high radial velocities. (i) A relatively low inclination angle of the outer disk is excluded by the observed moment 0 maps. García-Burillo et al. (2019) give an inclination angle i ≥ 60°, which is consistent with our results (Table 1). (iii) According to Burkert et al. (2010) the rotation velocity in the presence of a radial gas pressure gradient is

(18)

where v₀ is the zero-pressure velocity curve and the thermal pressure gradient neglected. For a model with constant Toomre parameter Q and vertical hydrostatic pressure equilibrium the velocity dispersion is constant and the density is given by ρ = Ω²/(π G Q). In the case of a constant rotation curve the pressure gradient is maximum:

(19)

With v₀ ∼ 100 km s⁻¹ and radial velocity dispersion of v_turb ∼ 50 km ⁻¹ assuming isotropic turbulence, one obtains v_rot = 70 km s⁻¹ or a velocity difference of Δv = 30 km s⁻¹. Within the radius of influence of the black hole, where v₀ corresponds to the Keplerian velocity, the velocity difference is smaller. This velocity difference can in principle explain the observed velocity difference. (iv) An intriguing possibility is that the emission of the high radial velocity gas is much weaker than that of the low radial velocity gas. This can be the case if the molecular emission is dominated by the PDR even in the outer parts of the disk because the high radial velocity gas is located in the plane of the sky and we therefore observe these illuminated surfaces almost edge-on. This orientation leads to a suppressed molecular line emission. This scenario implies that the parts of the outer gas disk (R > 2 pc), for which molecular line emission is detected, are directly illuminated by the AGN. In this case it is expected that deeper observations reveal gas of higher radial velocities.

We conclude that explanation (ii)–(iv) can explain the difference between the observed radial velocities and those expected from the mass profile. Given that the distribution of the molecular line emission outside R = 2 pc is asymmetric, we think that non-circular motions together with different emission conditions in the regions of high radial velocities might be the preferred explanations.

5.6. Plausibility of the proposed scenario

As stated in Sect. 1, our models are not intended to reproduce the NGC 1068 observations in detail, but to serve as a guideline to interpret the available observations. In the following we investigate if the proposed accretion scenario where gas clouds fall onto a preexisting gas disk or ring is plausible. It is certainly the simplest model configuration that can qualitatively, and to a certain degree quantitatively, explain the available molecular line observations. The mass accretion rates of the dynamical model are quite high (Ṁ ∼ 0.1–0.5 M_⊙ yr⁻¹; Fig. 14). The inspection of the simulation without infalling gas cloud showed that the outer edge of the ring moves from a radius of 7 pc to 3 pc within ∼0.5 Myr. This means that the gas cloud has to hit the preexisting gas disk or ring at the right moment. If the preexisting gas disk or ring had a larger size, the time window would be larger. Nevertheless, there has to be a certain timing between the accreting inner gas disk or ring and the infalling gas cloud.

We have a scenario in mind where there is permanent chaotic infall of gas clouds of different angular momentum into the central 20 pc. If there is no preexisting gas disk or ring, such structure will form through a dispersion ring (Sanders 1998). If a second gas clouds falls into the central 10 pc within the accretion time of the dispersion ring (∼1 Myr), we would find a NGC 1068-like scenario. The accretion time of the CND in the Galactic Center is a few Myr (Vollmer & Davies 2013). The time window for the infall of an external gas cloud can thus be about 4 times larger than that of NGC 1068. It is quite natural that the infall frequency of external gas clouds is significantly higher in AGNs compared to quiescent galactic centers.

The center of the Circinus galaxy also harbors an AGN with a maser disk around it (Greenhill et al. 2003). The ALMA CO, HCN, and HCO⁺ observations of Izumi et al. (2018) and Kawamuro et al. (2019) showed that the gas located within the inner 20 pc rotates in the same sense as the maser disk. The expected rotation velocity of the obscuring thick gas disk is clearly detected (Izumi et al. 2018). This disk is fed by a gaseous bar structure of a size of ∼100 pc. Thus, contrary to NGC 1068, there seems to be continuous infall over more than a dynamical timescale of t_dyn ∼ 100 pc/100 km s⁻¹ ∼ 1 Myr into the central 20 pc of the Circinus galaxy.

5.7. AGN variability

Infalling gas clouds thus seem to feed the gas disk found at distances of ∼10 pc around the central black hole. Due to these infall the mass accretion rate into the central ∼0.1 pc will be highly variable. As seen before, a mass accretion rate of Ṁ_BH ≲ 0.1 M_⊙ yr⁻¹ into the central 0.1 pc is sufficient to make a galactic nucleus active. Without the infall of an external gas cloud, this mass accretion rate can be maintained over ∼10 Myr. If a mass accretion higher than ≳0.3 M_⊙ yr⁻¹ is needed to make the galaxy nucleus active, the timescale becomes ∼0.2 Myr (see Fig. 14). The latter timescale is consistent with the flickering timescale found by Schawinski et al. (2015). However, we caution the reader not to overinterpret these findings because the interaction between the outer and inner disks and the outer and inner outflows, which set the accretion rate onto the central black hole giving rise to the AGN luminosity, is not clear at all.

6. Summary and conclusions

Recent ALMA high-resolution CO, HCO⁺, and HCN line observations of the nucleus of NGC 1068 found compact emission at a scale of ∼10 pc (Imanishi et al. 2018, 2020; García-Burillo et al. 2019; Impellizzeri et al. 2019). The emission distribution is elongated and shows signs of rotation. The interferometric observations of the highest angular resolutions revealed a counterrotating gas disk at distances ≲1.5 pc from the central black hole. The observed linewidths are large indicating that the gas disk is geometrically thick.

The gas disk in the nucleus of NGC 1068 has several properties in common with the Circumnuclear Disk (CND) in the Galactic Center: both are thick rotating gas rings, which are kinematically perturbed. The gas in the CND is resolved into clouds with sizes of ∼0.1 pc (Montero-Castaño et al. 2009) and densities between 10⁵ and 10⁶ cm⁻³ (Harada et al. 2015).

A simple dynamical model consisting of gas clouds that can have partially inelastic collisions is used to reproduce the main kinematical properties of the available observations (Sect. 2.1). The simplest configuration to explain the kinematics of the CND and the thick gas disk in NGC 1068 is gas clouds, which fall onto a preexisting gas ring. We made 162 simulations for the CND and 63 simulations for NGC 1068 and compared all meaningful projections to existing molecular line observations. The best-fit cloud-ring collision is prograde for the CND and retrograde for NGC 1068.

We applied the analytical model of a turbulent clumpy gas disk of Vollmer & Davies (2013) calibrated by the dynamical model and observations to the CND and the thick gas disk with the inner 20 pc of NGC 1068 (Sect. 2.2). Within the model framework the energy driving turbulence is supplied by external infall or the gain of potential energy by radial gas accretion within the disk. The external mass accretion rate is assumed to be close to the mass accretion rate within the disk (the external mass accretion rate feeds the disk at its outer edge). Within the model, the disk is characterized by the disk mass accretion rate Ṁ and the Toomre Q parameter, which is used as a measure of the gas content of the disk. Both quantities are assumed to be constant within the disk. The analytical model yields the average disk properties and those of the largest turbulent disk clouds whose size is given by the thickness of C-shocks within the partially ionized gas permeated by a relatively strong magnetic field.

The gas temperature is calculated through the local equilibrium between mechanical and cosmic ray heating on the one hand, and molecular line cooling (H₂O, H₂, CO) in the other hand (Sect. 2.3). The time-dependent molecule abundances are derived by the Nautilus gas-grain code (Sect. 2.4). Radial profiles of multi-transition CO, HCO⁺, and HCN brightness temperature were calculated.

To each disk gas cloud of the dynamical model the brightness temperature from the analytical model was assigned according to its radius. In this way datacubes of the molecular lines were produced, which are directly compared to the available observations via the moment 0 maps, pv diagrams, and spectra along the major axis. For simplicity we divided the turbulent clumpy gas into an intercloud gas and dense clouds. Furthermore, we evaluated the role of infrared pumping for the HCN emission and the role of X-ray dominated and photodissociation regions. No wind our outflow component is included in the model.

Our models are not expected to exactly reproduce the available observations, but to serve as a guideline to interpret these observations. Whereas traditional models derive the gas properties from molecular line ratios, our forward modeling tries to reproduce qualitatively and, to a certain extent, quantitatively the observed molecular brightness temperatures.

Based on our modeling effort we found the following results for the CND in the Galactic Center:

1. in the Galactic Center we suggest collisions of two gas clouds with a preexisting gas ring. This is certainly not a unique solution. The CND is observed ∼0.57 Myr and ∼0.15 Myr after the impacts (Sect. 3).

2. The CND-like model yields Q = 300 and has a gas mass of M_gas ∼ 10⁴ M_⊙ between 1 pc ≲ R ≲ 4 pc. The mass accretion rate of such a gas ring is a few 10⁻³ M_⊙ yr⁻¹.

3. Our model containing only dense disk clouds reproduces within a factor of two the available CO(6−5) (Requena-Torres et al. 2012), single dish HCN(3−2), HCN(4−3), HCO⁺(3−2), HCO⁺(4−3), H¹³CN(3−2), H¹³CN(4−3), H¹³CO⁺(3−2), H¹³CO⁺(4−3) (Mills et al. 2013), and interferometric HCN(4−3) and CS(7−6) (Montero-Castaño et al. 2009) observations of the CND (Sect. 3). The cosmic ray ionization rate within the CND is ζ_CR = 2 × 10⁻¹⁵ s⁻¹.

4. The HCN abundance is significantly enhanced in the inner 2 pc due to hot core physics.

5. To avoid self-absorption of the intercloud medium in the CND, the turbulent heating at the largest scales, which are comparable to the disk height, has to be decreased by a factor of ≳50 (Sect. 3.2). This might be achieved by intermittent turbulent heating and potentially by a strong uniform magnetic field.

The following results were found for NGC 1068:

1. The infalling cloud collides with the gas ring on a retrograde orbit. The thick gas disk (Q = 30) is observed ∼0.14 Myr after the impact (Sect. 4).

2. The model gas disk at 1.5 pc ≲ R ≲ 7 pc has a mass of M_gas ∼ 1.4 × 10⁵ M_⊙. The mass accretion rate of such a thick gas disk is ∼0.5 M_⊙ yr⁻¹.

3. Our model reproduces within a factor of ∼2 the available CO(2−1), CO(3−2), CO(6−5), HCN(3−2), HCN(4−3), HCO⁺(3−2), and HCO⁺(4−3), and CS(7−6) observations of the gas distribution in the inner 20 pc of NGC 1068 (Sect. 4). The cosmic ray ionization rate within the dense disk clouds in the thick gas disk is ζ_CR = 2 × 10⁻¹³ s⁻¹.

4. As for the CND in the Galactic Center, the HCN abundance is greatly enhanced due to hot core physics. The CS(7−6) observations of NGC 1068 can be reproduced with a sulfur abundance of S/H = 10⁻⁶.

5. Based on the available observations it cannot be excluded that the intercloud gas dominates the molecular line emission in NGC 1068 at radii ≲2 pc (Sect. 4.2.3). In order to obtain an optically thin intercloud gas, the Toomre Q parameter has to be set to Q = 100 and the dissipation rate at the driving length must be decreased by at least a factor of 1/200 compared to Eq. (6). This is four times less than the factor required for the CND. On the other hand, for an optically thick intercloud gas with brightness temperatures high enough to be detected by the available observation setups, we had to assume a Q = 30 disk with full heating (Eq. (6)) and diffuse clouds of size H/2 and density of three times the mean density.

6. Within the framework of the model turbulent mechanical energy input is the dominant gas-heating mechanism within the thick gas disk (Sect. 5.2). Turbulence is maintained by the gain of potential energy of the accreting gas.

7. A good part of the kinematical features of the molecular line emission in NGC 1068, which were interpreted as an outflow, might be due to the perturbed kinematics caused by the gas-gas collision of two counter-rotating entities. It is not excluded that the signatures of the gas-gas collision and the outflow are superimposed in the inner R ≲ 3 pc of NGC 1068.

8. Since HCN abundances in excess of x_HCN = 10⁻⁸ are absent for ionization rates that are higher than a few 10⁻¹² s⁻¹, the observed high HCN abundances in NGC 1068 cannot be reproduced with models including an ionization rate which is expected if the gas is directly illuminated by the X-ray emission of the central engine (Sect. 2.5). We suggest that the bulk of the central X-ray emission is absorbed in a gas layer of Compton-thick gas inside the dust sublimation radius, maybe by an outflow of dense fully ionized gas. In this case, NGC 1068 would harbor two outflows, an inner ionized and an outer dusty outflow (see Vollmer et al. 2018).

9. Infrared pumping increases the HCN(3−2) line emission of the whole gas disk or ring by up to a factor of two (Sect. 4.2.2).

10. It is expected that the central engine directly illuminates inner edge of the thick gas disk. Based on our model calculations we concluded that the molecular line emission can be enhanced by the PDR emission by at most a factor of two at radii between 1 and 2 pc (Sect. 4.2.4). We do not think that the PDR emission dominates the total emission at these radii.

The role of the intercloud gas in the thick disk and its turbulent heating has to be further investigated. The determination of the optical depths and brightness temperatures at all scales of the turbulent clumpy ISM as attempted by Vollmer et al. (2017) is certainly desirable. However, such a treatment requires a full 3D radiative transfer model for the molecular line emission. In the present model the radiative transfer is only calculated within the gas clouds and optically thick clouds can hide the emission of other clouds located behind them. We conclude that the scenario of infalling gas clouds onto preexisting gas rings in galactic centers is viable and consistent with available observation of the CND in the Galactic Center and the dense gas distribution within the inner 20 pc of NGC 1068. We suggest that the observed AGN flickering (Schawinski et al. 2015) might be linked to the proposed scenario of permanent chaotic infall of gas clouds of different angular momentum into the central 20 pc of galactic centers.

References

Appendix A: The analytical model

The analytical model is fully described in Vollmer & Davies (2013). Here we only give an overview. We assume a turbulent clumpy gas disks where the energy to drive turbulence is supplied by external infall or the gain of potential energy by radial gas accretion within the disk. The size of the largest turbulent gas clouds is determined by the size of a continuous (C-type) shock propagating in dense molecular clouds with a low ionization fraction at a given velocity dispersion. We use the expressions derived by Vollmer & Davies (2013) for the expected volume and area filling factors, mass, density, column density, and velocity dispersion of the clouds. The latter is based on scaling relations of intermittent turbulence whose open parameters are estimated for the Circumnuclear Disk in the Galactic Center. The meaning of all model parameters is given in Table A.1.

The averaged overall disk is treated as an accretion disk in hydrostatic equilibrium (e.g., Pringle 1981) with a given mass accretion rate Ṁ and Toomre parameter

(A.1)

where v_turb and v_rot are the turbulent and rotation velocities, and M_dyn and M_gas are the total enclosed mass and the gas mass. The Toomre Q parameter is used as a measure of the gas content of the disk with Q = 1 for the maximum disk gas mass. If the disk is in hydrostatic equilibrium the midplane density is given by ρ = Ω²/π/G/Q, where Ω = v_rot/R is the angular velocity.

Thermal instabilities triggered by small-scale gas compression lead to the condensation of dense cold clouds (Parravano 1987). Within the strong shock approximation the gas compression rate in the shock is given by

(A.2)

where the is the neutral Alfvén speed in the diffuse medium with magnetic field B (see, e.g., Chen & Ostriker 2012). The neutral Alfvén speed of the ISM clouds of different densities is almost constant v_A, 0 ∼ 1 km s⁻¹ (Crutcher 1999). Assuming recombination-ionization equilibrium, the shock thickness is given by Eq. 43 of Chen & Ostriker (2012):

(A.3)

where α = 3.7 × 10¹³ cm³s⁻¹g⁻¹ (Draine 1983) is the collision coefficient, x_i the ionization fraction, and (μ_i/μ_n) = (30/2.3) = 13 the fraction between the mean ion and neutral molecular weights. We assume a cloud radius r_cl = l.

The ionization fraction is given by

(A.4)

where γ = 600-2000 cms (Williams et al. 1998; McKee et al. 2010) and n_H = ρ/(2.3 × m_p). We adopted γ ∼ 600 cms. If the neutral Alfvén speed is higher, γ and ζ_CR have to be adapted such that (v_A,0/1 km s⁻¹) = , where ζ_CR, 0 is the adopted cosmic ray ionization rate.

The disk mean density is linked to the cloud density via the volume filling factor ρ = ρ_clΦ_V. An approximation for the volume filling factor is

(A.5)

In the case of compressible turbulence, the mean volume rate of energy transfer becomes (e.g., Fleck 1996). If the density scales with ρ_l ∝ l^−3α, the turbulent velocity dispersion obeys the relation v_l ∝ l^1/3 + α (Fleck 1996; Kritsuk et al. 2007). Extending the β-model to compressible turbulence leads to

(A.6)

With ρ_l ∝ l^−3α, the turbulent velocity dispersion obeys the relation . For D = 2 one obtains v_l ∝ l^α and the energy flux is conserved. In this case the scaling relation (Eq. A.6) becomes

(A.7)

With the area filling factor Φ_A = Φ_VH/r_cl = ρ/ρ_clH/r_cl one obtains

(A.8)

For each distance from the central black hole the size, density, and velocity dispersion of the turbulent gas clouds are determined using Eqs. A.3, A.5, and A.8, respectively. Each model is characterized by two input parameters: the Toomre Q parameter or the gas surface density and the turbulent velocity dispersion of the gas disk v_turb. These parameters where chosen such that the analytical models are consistent with the dynamical model at the time of interest. The comparison between the gas surface densities and gas velocity dispersions of the CND and NGC 10168 models are presented in Fig. A.1. Since the NGC 1068 model has an asymmetric mass distribution, we show surface density cuts along the line of sight; radial profiles are shown otherwise. The gas surface densities of the analytical models are on the higher end of but well comparable to the model surface densities. We decided to set the gas velocity dispersions of the analytical models to those derived from observations (Guesten et al. 1987; Hicks et al. 2009), which are about 25 % higher than those of the dynamical model.

Fig. A.1. Calibration of the analytical model (dashed lines) by the dynamical model (solid lines) at the time of interest. Panel (a) radial profile of the gas surface density of the CND-like model. Panel (b) radial profile of the gas velocity dispersion of the CND-like model. Panel (c) gas surface density cut along the line of sight of the NGC 1068-like model. Panel (d) radial profile of the gas velocity dispersion of the NGC 1068-like model.

Appendix B: Molecular line emission

We employ the emission line modeling used by Vollmer et al. (2017). A molecular line source is usually observed by chopping the telescope’s beam between on- and off-source positions and measuring the difference in antenna temperatures. In general, the difference in brightness temperatures is

(B.1)

where τ is the optical depth of the line, ν the frequency of the observations, h and k the Planck and Boltzmann constants, and T_ex and T_bg the excitation and background brightness temperatures, respectively.

Considering only a single collider (H₂) for simplicity, the excitation temperature is

(B.2)

where T_* = hν_ul/k, n is the gas density, nq_ul the collisional de-excitation rate, and A_ul the Einstein coefficients of the transition ul. The background brightness temperature T_bg is the sum of the effective emission temperatures of the galaxy’s dust T_{eff dust} and the cosmic background at the galaxy redshift T_CMB (see Eq. 17 of da Cunha et al. 2013).

Following Vollmer et al. (2017), we consider two-level molecular systems in which the level populations are determined by a balance of collisions with H₂, spontaneous decay and line photon absorption, and stimulated emission with τ > 1. The molecular abundances were calculated using the chemical network (see Sect. 2.4). For simplicity, we neglected the hyperfine structure of HCN.

The rotation constants, Einstein coefficients, and collision rates were taken from the Leiden Atomic and Molecular Database (LAMDA; Schöier et al. 2005. The CO collision rates were provided by Yang et al. (2010). The HCN collision rates were taken from the He–HCN rate coefficients calculated by Dumouchel et al. (2010), scaled by a factor of 1.36 to go to HCN–H₂ (see Green & Thaddeus 1976). The HCO⁺ collision rates were taken from Flower (1999).

The comparison of our results with those of RADEX (van der Tak et al. 2007) showed that our calculations tend to overestimate the CO optical depths by up to a factor of two and those of HCN and HCO⁺ by a factor of two to three. This is caused by the approximate treatment of the line emission, which is more accurate for the CO optical depth than for the HCN and HCO⁺ optical depths. This approximation was used because it gave us a high flexibility for the model calculations. Whereas the HCN and HCO⁺ brightness temperatures are underestimated by up to a factor of two, the CO brightness temperatures agree with those of RADEX with typical differences of about 10 to 30 %. Given the overall uncertainties of the analytical and chemical models, we think that these uncertainties are acceptable. We estimate the overall uncertainty of the brightness temperatures caused by the uncertainties of the analytical and chemical models and our approximate treatment of the molecular line emission to be about a factor of two. To be consistent with RADEX, we divided all model optical optical depths by a factor of two. We used RADEX for the calculations of the CS and CN optical depths and brightness temperatures.

Appendix C: HCN, HCO⁺, H¹³CN, and H¹³CO⁺ spectra along the major axis

We present the observed HCN(3−2) spectra from Impellizzeri et al. (2019) from resolution elements along the major axis in Fig. C.1. It is remarkable that the observed spectra outside the central resolution element have linewidths in excess of 150 km s⁻¹ and the velocities of their maxima are close to the systemic velocity, that is rotation is barely visible in these spectra. Moreover, the central resolution element shows the HCN(3−2) line in absorption over a velocity range between −400 and 150 km s⁻¹. The line is seen in emission around a radial velocity of ∼250 km s⁻¹.

Fig. C.1. NGC 1068. Left panel: HCN(3−2) high resolution spectra along the major axis (from Impellizzeri et al. 2019 data). Right panel: high-resolution HCN(3−2) spectra along the major axis of the dense disk cloud model with IR pumping.

The corresponding Q = 30 model HCN(3−2) spectra are presented in Fig. C.1. Whereas the emission peaks and the covered velocity range of the models are higher but comparable to the observed emission peaks, the shapes of the model lines are different from the observed line shapes: the model spectra have their maximum at ∼80 km s⁻¹ and show a secondary maximum at ∼ − 70 km s⁻¹, whereas the observed spectra have their maximum much closer to the systemic velocity (|Δv|≲50 km s⁻¹). The model HCN(3−2) line is seen in absorption between −100 and 0 km s⁻¹ in the central spectrum. The model absorption is thus much narrower and less deep than the observed absorption. Moreover, the model emission at positive velocities has no observed counterpart and there is no model emission at a velocity of 250 km s⁻¹ in the central spectrum. The signatures of rotation are thus more visible in the model than in the observations.

The high-resolution HCO⁺(3−2) spectra corresponding to Fig. C.1 are presented in Fig. C.2.

Fig. C.2. NGC 1068. Left panel: HCO+(3−2) high resolution spectra along the major axis (from Imanishi et al. 2020 data). Right panel: dense disk cloud model high-resolution HCO+(3−2) spectra along the major axis.

Fig. C.3. NGC 1068. Dense disk cloud model high-resolution H13CN(3−2) spectra along the major axis.

The main difference between the HCO⁺(3−2) and HCN(3−2) spectra of the model are ≲1.5 times lower flux densities of the HCO⁺(3−2) emission and a broader and deeper absorption around 200 km s⁻¹. The model HCO⁺(3−2)/HCN(3−2) ratio is significantly smaller for the peaks at ∼80 km s⁻¹ than for the rest of the spectra.

Imanishi et al. (2020) presented high-resolution H¹³CN(3−2) and H¹³CO⁺(3−2) observations of NGC 1068. At their sensitivity limit (see Table 2), no emission of the two lines was detected within the inner 20 pc of NGC 1068. Only the H¹³CN(3−2) line was seen in absorption within the central resolution element (Fig. 10 of Imanishi et al. 2020). The depth of the absorption is about −20 K and the velocity width is about Δv ∼ 250 km s⁻¹. We included the H¹³CN(3−2) line emission into our model by assuming an isotope ratio of ¹²C/¹³C = 30 (see Fig. 2 of Halfen et al. 2017). As expected, the model H¹³CN(3−2) emission is much fainter than the HCN(3−2) emission. It is only seen in emission (∼6σ) in the resolution element at a distance of 1.4 pc northwest of the central black hole. The H¹³CN(3−2) line is seen in absorption in the central resolution element, as it is observed. The model absorption depth is narrower (∼130 km s⁻¹) and about two times deeper than the observed absorption depth.

The model H¹³CO⁺(3−2) spectra (Fig. C.4) do not show any emission, as it is observed. Contrary to the observations, the H¹³CO⁺(3−2) line is seen in absorption in the model.

Fig. C.4. NGC 1068. Dense disk cloud model high-resolution H13CO+(3−2) spectra along the major axis.

We conclude that the model only qualitatively resembles the observed HCN(3−2) and HCO⁺(3−2) spectra. Given the tiny amount of H¹³CN(3−2) emission and the absence of H¹³CO⁺(3−2) emission in the model spectra, the assumed isotope ratio ¹²C/¹³C = 30 is justified.

Appendix D: ALMA observations

During ALMA Cycle 2 (project code 2013.1.00014.S), NGC 1068 was observed in Band 7 centered at λ ∼ 866 μm. We used the longest baseline configurations available. The baseline range was 15 - 2270 m. Structures larger than 7.2″ are filtered out. Total time on source was 18.5 minutes. Source observations were interleaved with short scans of nearby bright calibrator sources. The phase calibrators was J0239+0416, located 4.3° from NGC 1068. The spectral windows were tuned to observe three continuum windows, 2 GHz total bandwidth each, and one tuned to observe HCN (J = 4 → 3; ν₀ = 354.50547 GHz) at barycentric redshift 1150 km s⁻¹ (Brinks et al. 1997). The channel-width of the HCN line measurements is ∼0.4 km s⁻¹, spanning about ±800 km s⁻¹ total bandwidth. Calibration and data reduction followed standard procedures for ALMA observations in CASA. NGC 1068 proved bright enough for self-calibration over scan intervals. Two rounds of phase-only self-calibration reduced sidelobe artifacts and improved the background rms. The photometric accuracy, ∼10 %, is limited by the photometric uncertainties of the flux calibrators (3C454.3, J0224+0659, J0238+166, and J2253+1608). The absolute astrometric accuracy is limited by the transfer of calibration solutions from the complex gain calibrators to NGC 1068. The expected accuracy for the observed phase calibrators is about 5% of the beam width (Asayama et al. 2017), corresponding to 7 mas (0.5 pc).

Images and data cubes were produced using multiscale CLEAN (Cornwell 2008) and Briggs weighting (Briggs 1995). Spectral line channels were combined to create cubes with 10 km s⁻¹ channel-width. The background rms noise levels are 0.078 mJy beam⁻¹ and 1.1 mJy beam⁻¹ for the continuum and HCN channel maps, respectively. The synthetic restoring beam is 0.15″×0.11″, PA 78°. The primary beams is about 20″. As emission is detected only in the inner 3″, and our focus is on the inner arcsecond, the data are not corrected for primary beam attenuation. The expected attenuation at 1.5″ from the pointing center (S1) is less than 1% at either band.

We discovered additional line detections during the Band 7 continuum processing. Both CS(7−6) and the two strongest lines of the CN(3−2) triplet appeared in bands tuned for continuum. These data were processed following the same procedure as used for the CO and HCN observations. The final data cubes were binned into 15 km s⁻¹ channels during imaging and CLEAN deconvolution.

Table D.1 summarizes the integrated emission line properties, which were derived from a Gaussian fit to the emission line profile. For the purposes of the summary, we fitted the two strongest lines of the CN(3−2) triplet.

All Tables

All Figures

	Fig. 1. Fraction of dissipated energy per cloud-cloud collision as a function of the distance to the central black hole. The dissipation energy fraction is shown for all clouds during the NGC 1068 model evolution presented in Fig. 3.
In the text

	Fig. 2. Evolution of the CND-like model. The times since the beginning of the simulations are indicated. Upper half: RA-Dec projection. Lower half: RA-LOS (line-of-sight) projection.
In the text

	Fig. 3. Evolution of the NGC 1068-like model. The times since the beginning of the simulations are indicated. Upper half: RA-Dec projection. Lower half: RA-LOS (line-of-sight) projection.
In the text

	Fig. 4. CS(2−1) emission distribution. Left panel: IRAM 30 m observations of the CND in the Galactic Center (Güsten et al., priv. comm.). Right panel: best-fit combined model (see Table 1). The levels are 13 to 260 K km s−1 in steps of 13 K km s−1. Both images have the same physical size.
In the text

	Fig. 5. HCN abundance xHCN = nHCN/nH as a function of density and temperature for different cosmic ray ionization fractions. Panel a: ζCR = 2 × 10−13 s−1, Panel b: 2 × 10−12 s−1, and panel c: 2 × 10−11 s−1.
In the text

	Fig. 6. Characteristics of individual model disk gas clouds across the CND at the time of interest. From upper left to lower right: cloud mass, size, density, H2 column density, gas and dust temperatures, optical depth, and brightness temperature.
In the text

	Fig. 7. CO(6−5) emission distribution. Upper panel: model; lower panel: APEX/CHAMP+ observations integrated between −150 and 150 km s−1 (Requena-Torres et al. 2012). The spatial resolution is 9.5″. The flux is given in units of K km s−1, in Tmb.
In the text

	Fig. 8. CO spectra. Upper panel: CO(6−5) model spectra of the southern and northern lobes. Lower panel: CO spectra toward the southern (left) and northern (right) lobe of the CND (from Requena-Torres et al. 2012). All spectra are convolved to the same angular resolution of 22.5″. The flux density is given in units of K, in Tmb.
In the text

Fig. 9. HCN(4−3) integrated intensity map. Upper panel: model. The contour levels are in steps of 53 K km s−1 from 27 to 875 K km s−1 The spatial resolution is 4.6″ × 3.0″. Lower panel: observations from Montero-Castaño et al. (2009). The rms is 5.9 K km s−1. The contour levels are in steps of 35 K km s−1 from 18 to 550 K km s−1, but for the highest contour level at 574 K km s−1. Sgr A* is marked with a star.

In the text

Fig. 10. CS(7−6) integrated intensity map. Upper panel: model. The contour levels are in steps of 22 K km s−1 from 11 to 275 K km s−1. The spatial resolution is 4.5″ × 3.1″. Lower panel: observations from Montero-Castaño et al. (2009). The rms is 2.5 K km s−1. The contour levels are in steps of 15 K km s−1 from 7 to 170 K km s−1, but for the highest contour level at 190 K km s−1. Sgr A* is marked with a star.

In the text

Fig. 11. HCN(3−2), HCN(4−3), HCO+(3−2), and HCO+(4−3) moment 0 maps. Panels a–d: model. Contour levels are from 100 to 500 K km s−1, spaced by 50 K km s−1. Panels e–h: observations (Mills et al. 2013). The spatial resolutions are: HCN(3−2): 23.6″, HCN(4−3): 17.7″, HCO+(3−2): 23.5″, HCO+(4−3): 17.6″. Contours are linearly spaced. (e) HCN(3−2), contours from 75 to 265 K km s−1, spaced by 27 K km s−1. (f) HCO+(3−2), contours from 58 to 140 K km s−1, spaced by 17 K km s−1. (g) HCN(4−3), contours from 60 to 262 K km s−1, spaced by 29 K km s−1. (h) HC0+(4−3), contours from 58 to 140 K km s−1, spaced by 17 K km s−1.

In the text

	Fig. 12. Model H13CN(3−2) (a), H13CN(4−3) (b), H13CO+(3−2) (c), and H13CO+(4−3) (d) moment 0 maps. The contour levels are those of Fig. 11 divided by a factor 4.
In the text

	Fig. 13. Characteristics of individual model intercloud gas across the CND at the time of interest. From upper left to lower right: cloud mass, size, density, H2 column density, gas and dust temperatures, optical depth, and brightness temperature.
In the text

	Fig. 14. Model mass accretion rate measured at R = 0.25 pc. Solid line: preferred model with ξ = 0.9, dashed line: preferred model without the massive infalling gas cloud, vertical line: time of interest.
In the text

	Fig. 15. NGC 1068 Q = 30 analytical model of a turbulent clumpy gas disk (Sect. 2.2). Characteristics of individual model disk gas clouds across the disk at the time of interest. From upper left to lower right: cloud mass, size, density, H2 column density, gas and dust temperatures, optical depth, and brightness temperature. IR pumping is included for the HCN and HCO+ lines.
In the text

Fig. 16. NGC 1068. Panel (a) observed CO(2−1) moment 0 map (I0 = 7923 K km s−1), (b) pv diagram (contours are Tb0 = 37.8 K) from García-Burillo et al. (2019), (c) dense disk cloud model CO(2−1) moment 0 maps, (d) pv diagrams along the major axis, (e) observed CO(3−2) moment 0 map (I0 = 7923 K km s−1), (f) pv diagram (Tb0 = 37.8 K) from García-Burillo et al. (2019), and (g) dense disk cloud model CO(3−2) moment 0 maps, (h) pv diagrams along the major axis. The model levels are the same as those of the observed maps.

In the text

	Fig. 17. NGC 1068. Panel (a) observed CO(6−5) moment 0 map (Gallimore et al. 2016) (I0 = 6657 K km s−1), (b) CO(6−5) pv diagram along the major axis (Tb0 = 38.0 K), (c) dense disk cloud model CO(6−5) moment 0 maps, and (d) pv diagrams along the major axis. The model levels are the same as those of the observed maps.
In the text

Fig. 18. NGC 1068. Left half: high resolution HCN(3−2) observations (Impellizzeri et al. 2019). Panel (a) moment 0 map (I0 = 17300 K km s−1), (b) pv diagram along the major axis, (c) pv diagram along the minor axis (Tb0 = 70 K), (d) moment 1 map (in km s−1), (e) continuum map (contour levels are (0.019, 0.037, 0.075, 0.15, 0.3, 0.6, 1, 2, 4, 8)×43 K). Right half: high resolution HCN(3−2) emission of the dense disk cloud model with IR pumping. Panel (f) moment 0 map, (g) pv diagram along the major axis, (h) pv diagram along the minor axis, (i) unsmoothed pv diagram along the major axis, (j) moment 1 map (in km s−1), (k) continuum map. The model levels are the same as those of the observed maps. The additional white contour is at 80 K.

In the text

Fig. 19. NGC 1068. Left half: high resolution HCO+(3−2) observations (Imanishi et al. 2020). Panel (a) moment 0 map (I0 = 18420 K km s−1), (b) moment 1 map (in km s−1), (c) pv diagram along the major axis, and (d) pv diagram along the minor axis (Tb0 = 75 K). Right half: dense disk cloud model. Panel (e) moment 0 map, (f) moment 1 map (in km s−1), (g) pv diagram along the major axis, and (h) pv diagram along the minor axis. The model levels are the same as those of the observed maps.

In the text

	Fig. 20. NGC 1068. Panel (a) observed HCO+(4−3) moment 0 map (I0 = 7217 K km s−1), (b) pv diagram (Tb0 = 47.5 K) from García-Burillo et al. (2019), (c) dense disk cloud model HCO+(4−3) moment 0 maps, and (d) pv diagrams along the major axis. The model levels are the same as those of the observed maps.
In the text

	Fig. 21. IR pumping. Left half: high resolution HCN(3−2) emission of the dense disk cloud model. Panel (a) moment 0 map, (b) moment 1 map, (c) pv diagram along the major axis, and (d) pv diagram along the minor axis. The model levels are the same as those of the observed maps (Fig. 18). Right half: model high-resolution HCN(3−2) spectra along the major axis.
In the text

	Fig. 22. NGC 1068 intercloud gas model. Characteristics of individual model intercloud gas across the disk at the time of interest. Q = 30, ζCR = 2 × 10−13 s−1, and full turbulent heating. Gas cloud characteristics. From upper left to lower right: cloud mass, size, density, H2 column density, gas and dust temperatures, optical depth, and brightness temperature.
In the text

	Fig. 23. NGC 1068 intercloud gas model. Left half: HCN(3−2); right half: HCO+(3−2). Panels (a) and (e) model moment 0 maps, (b) and (f) model moment 1 maps, (c) and (g) pv diagrams along the major axis, (d) and (h) pv diagrams along the minor axis. The model levels are the same as those of the observed maps (Figs. 18 and 19), except for the white contour at 80 K.
In the text

	Fig. 24. NGC 1068 intercloud gas model. Panels (a) and (b) HCO+(4−3) moment 0 maps and pv diagrams, (c) and (d) CO(2−1) moment 0 maps and pv diagrams, (e) and (f) CO(3−2) moment 0 maps and pv diagrams, (g) and (h) CO(6−5) moment 0 maps and pv diagrams. The model levels are the same as those of the observed maps (Figs. 16, 17, and 20).
In the text

	Fig. 25. Radial profiles of the molecular line emission from the dense clouds (thin lines) and the associated PDR regions (thick lines). The CR ionization rate is ζCR = 10−14 s−1.
In the text

	Fig. 26. Radial profiles of the molecular line emission from the intercloud gas (thin lines) and the associated PDR regions (thick lines). The CR ionization rate is ζCR = 10−14 s−1.
In the text

	Fig. 27. ALMA integrated spectra of the central source S1 of NGC 1068 (black line). Upper panel: HCN(4−3); middle panel: CN(3−2); lower panel: CS(7−6). The gray line corresponds to the second component of the CN(3−2) triplet. Colored lines: dashed: S1 spectrum from the dense disk cloud model; dotted: S1 spectrum from the diffuse gas model.
In the text

Fig. 28. NGC 1068 dense disk cloud model without absorption. Left panel: Q = 30 model high resolution HCN(3−2) emission of the dense disk cloud model. Right panel: same as left panel but with IR-pumping included in the model. Panels (a) and (e) moment 0 maps, (b) and (f) moment 1 maps, (c) and (g) pv diagrams along the major axis, (d) and (h) pv diagrams along the minor axis. The model levels are the same as those of the observed maps (Fig. 18). The additional white contours levels are 80 and 90 K.

In the text

Fig. A.1. Calibration of the analytical model (dashed lines) by the dynamical model (solid lines) at the time of interest. Panel (a) radial profile of the gas surface density of the CND-like model. Panel (b) radial profile of the gas velocity dispersion of the CND-like model. Panel (c) gas surface density cut along the line of sight of the NGC 1068-like model. Panel (d) radial profile of the gas velocity dispersion of the NGC 1068-like model.

In the text

	Fig. C.1. NGC 1068. Left panel: HCN(3−2) high resolution spectra along the major axis (from Impellizzeri et al. 2019 data). Right panel: high-resolution HCN(3−2) spectra along the major axis of the dense disk cloud model with IR pumping.
In the text

	Fig. C.2. NGC 1068. Left panel: HCO+(3−2) high resolution spectra along the major axis (from Imanishi et al. 2020 data). Right panel: dense disk cloud model high-resolution HCO+(3−2) spectra along the major axis.
In the text

	Fig. C.3. NGC 1068. Dense disk cloud model high-resolution H13CN(3−2) spectra along the major axis.
In the text

	Fig. C.4. NGC 1068. Dense disk cloud model high-resolution H13CO+(3−2) spectra along the major axis.
In the text