Modeling the Vertical Structure of Nuclear Starburst Discs: A Possible Source of AGN Obscuration at $z\sim 1$

TL;DR

This study models the vertical structure of nuclear starburst discs in active galaxies at intermediate redshifts, revealing their potential to obscure AGNs and produce various observed types depending on viewing angles.

Contribution

It introduces a 2D modeling approach for NSDs, improving upon previous 1D models by providing detailed vertical solutions and statistical predictions of AGN obscuration.

Findings

2D models predict less atmospheric expansion than 1D models.

Approximately 56% of models suggest Type 1 AGN appearance.

Distribution of column densities can vary widely with viewing angle.

Abstract

Nuclear starburst discs (NSDs) are star-forming discs that may be residing in the nuclear regions of active galaxies at intermediate redshifts. One dimensional (1D) analytical models developed by Thompson et al. (2005) show that these discs can possess an inflationary atmosphere when dust is sublimated on parsec scales. This make NSDs a viable source for AGN obscuration. We model the two dimensional (2D) structure of NSDs using an iterative method in order to compute the explicit vertical solutions for a given annulus. These solutions satisfy energy and hydrostatic balance, as well as the radiative transfer equation. In comparison to the 1D model, the 2D calculation predicts a less extensive expansion of the atmosphere by orders of magnitude at the parsec/sub-parsec scale, but the new scale-height $h$ may still exceed the radial distance $R$ for various physical conditions. A total of 192 NSD models are computed across the input parameter space in order to predict distributions of a line of sight column density $N_H$. Assuming a random distribution of input parameters, the statistics yield 56% of Type 1, 23% of Compton-thin Type 2s (CN), and 21% of Compton-thick (CK) AGNs. Depending on a viewing angle ($\theta$) of a particular NSD (fixed physical conditions), any central AGN can appear to be Type 1, CN, or CK which is consistent with the basic unification theory of AGNs. Our results show that $\log[N_H(\text{cm}^{-2})]\in$ [23,25.5] can be oriented at any $\theta$ from 0$^\circ$ to $\approx$80$^\circ$ due to the degeneracy in the input parameters.

Figures

Tables

Table 1. Table 1: Description of symbols which are used frequently in the paper.

Symbol	Description
$c$	Speed of light
$c_s$	Total speed of sound
$f\text{<sub><span>CK</span></sub>}$	Fraction of Compton-thick AGNs
$f\text{<sub><span>CN</span></sub>}$	Fraction of Compton-thin AGNs
$f\text{<sub><span>g,out</span></sub>}$	Gas fraction at the outer radius
$G$	Gravitational constant
$h$	Surface scale-height
$\kappa$	Opacity
$k_B$	Boltzmann constant
$m$	Radial Mach number
$M\text{<sub><span>bh</span></sub>}$	A black hole mass
$m_p$	Mass of proton
$N_H$	Column density along a line of sight
$\Omega$	Keplerian angular frequency of a disc
$P\text{<sub><span>gas</span></sub>}$	Gas pressure
$P\text{<sub><span>rad</span></sub>}$	Radiation pressure
$P\text{<sub><span>turb</span></sub>}$	Turbulence pressure
$P$	Total pressure
$Q$	Toomre parameter
$\rho$	Density
$R$	Radial distance at $z=0$
$R\text{<sub><span>cric</span></sub>}$	Distance where the largest inflation occurs in a disc
$R\text{<sub><span>out</span></sub>}$	Size of a disc
${\sigma }_{\text{sb}}$	Stefan-Boltzmann constant
$\Sigma$	Vertically mass integrated column density
$\Sigma \text{<sub><span>mp</span></sub>}$	Vertically mass integrated total column density to mid-plane
$\sigma$	Velocity dispersion
$\dot{\Sigma }\text{<sub>*</sub>}$	Star-formation rate density
$T$	Temperature
$\tau$	Optical depth
$T\text{<sub><span>eff</span></sub>}$	Effective temperature
$\theta$	Viewing angle
$v\text{<sub><span>turb</span></sub>}$	Turbulence speed
$z$	Vertical height

Table 2. Table 2: We list all the bins in the domain of each input parameter ( M ​ bh 𝑀 bh M\textsubscript{\text{bh}} , R ​ out 𝑅 out R\textsubscript{\text{out}} , f ​ g,out 𝑓 g,out f\textsubscript{\text{g,out}} , m 𝑚 m ). After taking all the possible combination of these bins, a total of 192 models are computed. The Mach number m 𝑚 m is chosen such that gas can accrete to the parsec scale and the other parameters represent a possible range of physical conditions for NSDs.

log($M\text{<sub><span>bh</span></sub>}/M_{☉}$)	$R\text{<sub><span>out</span></sub>}$(pc)	$f\text{<sub><span>g,out</span></sub>}$	$m$
6.5	240	0.2	0.1
7.0	180	0.4	0.3
7.5	120	0.6	0.5
8.0	60	0.8	–

Table 3. Table 3: A total 99 models out of 192 show the starburst phenomenon causing the large atmosphere expansion at the parsec scale. Here, we show a distribution of these models across the input parameter space. This shows that a smaller disc size, a larger gas fraction and Mach number are in favor of a burst of star-formation.

	log($M\text{<sub><span>bh</span></sub>}/M_{☉}$)	6.5	7.0	7.5	8.0	Total (99)
$R\text{<sub><span>out</span></sub>}$(pc)	60	7	6	9	10	32
	120	5	6	6	8	25
	180	5	4	6	6	21
	240	5	5	5	6	21
$f\text{<sub><span>g,out</span></sub>}$	0.2	0	0	1	3	4
	0.4	5	6	7	8	26
	0.6	8	8	9	9	34
	0.8	9	7	9	10	35
$m$	0.1	1	1	2	3	7
	0.3	9	10	11	13	43
	0.5	12	10	13	14	49

Table 4. Table 4: 2D interpolation of n H subscript 𝑛 𝐻 n_{H} and T 𝑇 T . The resolution among temperature points is taken to be 0.01 except for ranges of [3.13,3.26] and [3.243,3.2431] which has resolutions of 0.005 and 0.0001, respectively.

Region	I	A	B	C	D
$\log \left(n_H\right)$	Independ.	[3.0,6.0]	[6.0,13.0]	[3.0,13.0]	[-7.0,3.0]
$\Delta$ $\log \left(n_H\right)$	N/A	0.25	0.5	0.5	0.25
$\log \left(T\right)$	[1.0,3.0]	[3.0,4.1]	[3.0,4.1]	[4.1,6.1]	[3.0,6.1]

Equations

$$\displaystyle\Omega(R)=\Big{(}\frac{GM_{\text{bh}}}{R^{3}}+\frac{2\sigma^{2}}{R^{2}}\Big{)}^{1/2},$$

$$\displaystyle M_{\text{bh}}=2\times 10^{8}\Bigg{(}\frac{\sigma}{200\text{ km s}^{-1}}\Bigg{)}^{4}M_{\astrosun}.$$

$$\dot{M}(R)={\dot{M}}_{\text{out}}-{\int }_R^{R_{\text{out}}}2\pi R^{\prime }{\dot{\Sigma }}_{\ast }(R^{\prime })d{R}^{\prime }$$

$$\dot{M_{\ast }}(R)=\pi R^2{\dot{\Sigma }}_{\ast }(R),$$

$$\displaystyle h_{\text{ph}}\approx h_{\text{mid}}\Big{(}\frac{\kappa_{\text{ph}}}{\kappa_{\text{mid}}}\Big{)},$$

$$\mu \frac{\delta I_{\nu }(z,\mu )}{\delta z}={\eta }_{\nu }(z,\mu )-{\chi }_{\nu }(z,\mu )I_{\nu }(z,\mu ).$$

$$\displaystyle T^{4}(\tau)=\frac{3}{4}T^{4}_{\text{eff}}\Big{[}\tau+q(\tau)\Big{]}$$

$$\tau (\Sigma )$$

$$z(\Sigma )={\int }_{\Sigma \text{textsubscript}\text{mp}}^{\Sigma }-\frac{1}{\rho ({\Sigma }^{\prime })}d{\Sigma }^{\prime },$$

$$T=(T_{\ast }^4+T_{\text{vis}}^4{)}^{1/4}.$$

$$\kappa =\kappa (\rho ,T)$$

$$\frac{d^{2}P(\Sigma )}{d{\Sigma }^2}=-\frac{c_s^2(\Sigma ){\Omega }^2}{P(\Sigma )}-4\pi G,$$

$$P\text{textsubscript}\text{mp}(\Sigma \text{textsubscript}\text{mp})=P\text{textsubscript}\text{mp}\text{and}{P}_1({\Sigma }_1)=P_1.$$

$$c_s^2(\Sigma )=\frac{P(\Sigma )}{\rho (\Sigma )}.$$

$$P(\Sigma )=P\text{textsubscript}\text{gas}(\Sigma )+P\text{textsubscript}\text{rad}(\Sigma )+P\text{textsubscript}\text{turb}(\Sigma ).$$

$$P\text{textsubscript}\text{gas}(\Sigma )=\frac{\rho (\Sigma )k_{B}T(\Sigma )}{m_p}$$

$$P\text{textsubscript}\text{rad}(\Sigma )=\frac{4{\sigma }_{\text{sb}}}{3c}{T}^4(\Sigma ),$$

$$P\text{textsubscript}\text{turb}(\Sigma )$$

$$=\frac{1}{2}\rho (\Sigma )v\text{textsubscript}{\text{turb}}^2,$$

$$Q(\Sigma \text{textsubscript}\text{mp})=\frac{c_s(\Sigma \text{textsubscript}\text{mp})\Omega }{2\pi G\Sigma \text{textsubscript}\text{mp}}=\mathrm{1.0.}$$

$$c_s\to c_s+\delta c_s\text{where}\delta c_s=\frac{1}{2}\frac{\delta P\text{textsubscript}\text{E}}{c_s\rho }.$$

$$N_H\approx {\int }_{R_{\text{in}}}^{R_{\text{out}}}\rho (\theta ,R){\cos }^{-1}(\theta )dR/m_p.$$

$$\displaystyle r_{s}=1.3\Bigg{(}\frac{L_{\text{UV}}}{10^{46}\text{erg s\textsuperscript{-1}}}\Bigg{)}^{1/2}\Bigg{(}\frac{T_{s}}{1500\text{K}}\Bigg{)}^{-2.8}\text{pc}.$$

$$10\text{K}\le T\le 10^{6.1}\text{K}$$

$$10^{-7}\text{cmtextsuperscript-3}\le n_H\le 10^{13}\text{cmtextsuperscript-3}$$

$$\text{\%Error}\equiv \frac{100\times (\text{Predicted value-Actual value)}}{\text{Actual value}}.$$

Full text

X-rays:diffuse background

††pagerange: Modeling the Vertical Structure of Nuclear Starburst Discs: A Possible Source of AGN Obscuration at $z\sim 1$–A††pubyear: 2017††pagerange: Modeling the Vertical Structure of Nuclear Starburst Discs: A Possible Source of AGN Obscuration at $z\sim 1$–A

1 Introduction

Active galactic nuclei (AGNs) are extremely luminous regions, more than a billion times the solar luminosity, in galaxies. They emit energy across a large portion of the electromagnetic spectrum from radio waves to X-rays through accretion onto a supermassive black hole (Balbus, 2003). The majority of AGNs discovered in the thousands by Chandra and XMM-Newton show the presence of obscuration by dust and gas (Brandt & Hasinger, 2005; Gilli et al., 2007; Alexander & Hickox, 2012; LaMassa et al., 2016). Infrared (Thatte et al., 1997; Burtscher et al., 2013) and X-ray based observations (Elvis et al., 2004; Risaliti et al., 2005) point the location of this obscuration at the parsec scale. The simple unification scheme of AGNs suggests that all the variance in the obscuration are due to the difference in viewing angle ($\theta$) of AGNs with respect to an observer, but they are intrinsically the same objects (Antonucci, 1993; Netzer, 2015). The mechanisms and physical conditions responsible behind this obscuration is still not well understood, but they seem to be dependent on AGN luminosity (Lawrence & Elvis, 2010) and, perhaps, redshift (Ballantyne et al., 2006; Hasinger, 2008).

To explain the obscuration at the parsec scale, many possible mechanisms have been proposed. The detection of polarized broad lines in Seyfert 2 galaxies (e.g., Antonucci & Miller, 1985) led to the modeling of a simple uniform toroidal absorber (Krolik & Begelman, 1986, 1988; Pier & Krolik, 1992). In addition, many groups have proposed a geometrically thick torus supported by infrared (IR) radiation pressure (Krolik, 2007; Dorodnitsyn et al., 2011; Dorodnitsyn & Kallman, 2012; Chan & Krolik, 2016; Dorodnitsyn et al., 2016). Another possible model is a torus with turbulence pressure due to supernovae or stellar winds (Wada & Tomisaka, 2005; Watabe & Umemura, 2005). Moreover, many authors have suggested a warped/tilted-disc to explain a range of obscuration (Nayakshin, 2005; Caproni et al., 2006; Lawrence & Elvis, 2010). In recent years, there has been some work done in modeling a clumpy torus (Hönig & Beckert, 2007; Nenkova et al., 2008; Hönig & Kishimoto, 2010) partly due to the observed column density variability seen in X-rays (Risaliti, 2002). Another possible source for the observed obscuration is nuclear starburst regions (Fabian et al., 1998; Wada & Norman, 2002; Thompson et al., 2005; Ballantyne, 2008; Hopkins et al., 2016). In this paper, we explore the effect of nuclear starburst discs (NSDs) on the AGN obscuration.

More than 50 years has passed since the discovery of the cosmic X-ray background (CXB) (Giacconi et al., 1962) and it is still not entirely resolved into discrete sources. The CXB is the spectral energy distribution (SED) which is characterized by a power law with the photon index of $\Gamma=1.4$ in the 2-10 keV band and a peak in the $10-30$ keV band (Gruber et al., 1999). The SED up to a few keV energies can be explained by integrating point-like sources (Worsley et al., 2005) and many of these sources are observed to be AGNs (e.g., Mushotzky et al., 2000; Bauer et al., 2004). Understanding the distribution of column density along a line of sight $N_{H}$ is an important aspect in modeling the CXB spectrum (Gilli et al., 2007) since the relative fraction of unobscured (Type 1, $N_{H}$ $<$ 1022 cm-2), Compton-thin Type 2s (CN, 1022 cm-2 $\leq$ $N_{H}$ $<$ 1024 cm-2), and Compton-thick (CK, $N_{H}\geq$ 1024 cm-2) AGNs determine the shape of the spectrum. In particular, the observed peak of the CXB spectrum requires a significant number of CK AGNs at moderate redshifts (Comastri et al., 1995; Gilli et al., 2001; Churazov et al., 2007; Moretti et al., 2009). Many studies have shown that 10% to 25% of CK AGNs are required in order to produce the observed peak of the CXB near 30 keV (Draper & Ballantyne, 2009; Akylas et al., 2012; Ueda et al., 2014). However, Akylas et al. (2012) shows that the fraction of CK AGNs ($f\textsubscript{\text{CK}}$) can range from 5% to 50% due to the degeneracy of the input parameters in the modeling of the CXB spectrum. If NSDs are the potential source for producing the CXB peak, we can place a constraint on $f\textsubscript{\text{CK}}$ by modeling the 2D structure of NSDs.

Many reasons point toward NSDs as a potential source of obscuration. Young stellar populations are found in the inner regions of nearby AGNs (e.g., González Delgado et al., 2001; Gu et al., 2001; Cid Fernandes et al., 2004; Storchi-Bergmann et al., 2005; Ruschel-Dutra et al., 2017) including Seyfert galaxies (Davies et al., 2007). Davies et al. (2007) highlights a possible causality between star formation (SF) and AGN activity from studies of nine nearby Seyfert galaxies: the AGN activity occurs later in time separated by 50-200 Myrs from the peak of the SF rate. Such a causality suggests a possible fueling of gas toward a central black hole by stellar winds and supernovae (e.g., Vollmer et al., 2008; Hopkins, 2012), perhaps indicating a strong coupling of SF, stellar winds/supernovae, and the AGN activity. With a simple 1D model of NSDs, Thompson et al. (2005) shows a region near dust sublimation can inflate to $h\sim R$ which can obscure the incoming AGN irradiation. Such a mechanism can build a bridge among star-forming regions, AGN activity, and the obscuration. These reports provide further motivation to study NSDs in detail. Multi-dimensional modeling of NSD is important in order to resolve the dependency of $N_{H}$ on the orientation angle $\theta$ of AGNs with respect to an observer.

In order to place a proper constraint on the column density $N_{H}$ and the fractions of CN Type 2s ($f\textsubscript{\text{CN}}$) and CK AGNs, the vertical structure of NSDs is required. In the past, the hydrostatic structure of a disc has been computed for accretion discs using an iterative method (Hubeny, 1990). Using similar procedures as described by Hubeny (1990) and Hubeny & Hubeny (1998), we compute 2D structure of NSDs under various physical conditions. Later, we study the dependency of $N_{H}$ on $\theta$ and also place a constraint on $f\textsubscript{\text{CN}}$ and $f\textsubscript{\text{CK}}$ using these models. These are outcomes of the NSD theory which can be used to test whether these discs can plausibly obscure and fuel Seyfert galaxies at $z\sim 1$. (Starburst discs may be a more prominent source for AGN obscuration when a large gas fraction is available in galaxies (Ballantyne, 2008)-conditions that are more favorable at intermediate redshifts.)

Sect. 2 reviews the one-dimensional theory of NSDs. Sect. 3 describes the methodology of vertical structure. In Sect. 4, we provide and discuss the results including caveats and future work. Then, the paper is concluded in Sect. 5. Appendix A provides details on the opacity calculation used in the modeling.

2 A Brief Overview of Radial Structure

The one dimensional (radial structure) model of NSDs in the radiation dominated regime has been computed by Thompson et al. (2005) for ultra luminous infrared galaxies (ULIRGs). These models require four input parameters: size of the disc $R\textsubscript{\text{out}}$, the black hole mass $M\textsubscript{\text{bh}}$, the gas fraction at the outer radius $f\textsubscript{\text{g,out}}$, and the Mach number $m\equiv v_{r}/c_{s}$, where $v_{r}$ is the radial velocity and $c_{s}$ is the total speed of sound. (Variables which are frequently used in the paper are listed in Table 1.) The 1D model has a NSD rotating with the Keplerian frequency $\Omega$ around a central black hole with the mass of $M\textsubscript{\text{bh}}$. The $\Omega$ includes the potentials of the black hole and the bulge,

[TABLE]

where $G$ is the gravitational constant and $R$ is the radial distance. The dispersion $\sigma$ is computed using the $M_{\text{bh}}-\sigma$ relationship (Ferrarese & Merritt, 2000; Gebhardt et al., 2000; Tremaine et al., 2002),

[TABLE]

The stability of a disc against gravitational collapse is always governed by Toomre’s parameter $Q$. (The parameter $Q$ is maintained to be 1 because there is an observational evidence in infrared ultraluminous galaxies (Downes & Solomon, 1998) and local spiral galaxies (Martin & Kennicutt, 2001; Leroy et al., 2008; Westfall et al., 2014) including the Milky Way (Binney & Tremaine, 1987; Rafikov, 2001) for $Q\simeq$1.) A gas fraction is injected at the outer radius of the disc which is driven towards the black hole at the constant Mach number through a presumed global torque (which can be due to a spiral instability or a bar; Goodman (2003)). The amount of available gas decreases with the radius as gas depletes into the formation of stars in order to maintain the Toomre stability criteria $Q=1$. Then, the accretion rate ($\dot{M}$) and the SF rate ($\dot{M_{*}}$) are given by

[TABLE]

and

[TABLE]

where $\dot{\Sigma}_{*}$ is the SF rate density. The radiation pressure from this SF provides a vertical support for the dusty atmosphere against gravity. At some inner radius $R_{\text{in}}$, the energy flux production from viscosity takes over the SF energy flux and the $R_{\text{in}}$ is considered as a transition radius between a NSD and an accretion disc (AD).

In this 1D model, the vertical structure is not computed except the scale-height of a photosphere $h_{\text{ph}}$ is approximated by scaling the mid-plane scale-height $h_{\text{mid}}$ with the ratio of $\kappa_{\text{ph}}$ and $\kappa_{\text{mid}}$:

[TABLE]

where $\kappa_{\text{ph}}$ and $\kappa_{\text{mid}}$ are the opacities at the photosphere and mid-plane, respectively. $\kappa_{\text{ph}}$ is computed by assuming that the temperature at the surface is the same as the effective temperature. With this theory, Thompson et al. (2005) showed that a large expansion of an atmosphere, $h_{\text{ph}}\sim R$, is possible on parsec scales for certain conditions of ULIRGs. Using the same theory, Ballantyne (2008) reaches similar conclusions for “Seyfert-like" conditions.

Here, we relax these assumptions and compute the exact vertical structure of NSDs by solving the coupled equations of hydrostatic balance, radiative transfer, and energy balance.

3 2D Structure

For time-independent conditions, radial and vertical structures are not possible to solve simultaneously due to a lack of prior knowledge on the distribution of any physical quantities. Therefore, the distributions of total energy flux and mass-integrated column density at the mid-plane ($\Sigma\textsubscript{\text{mp}}$) are computed using the Thompson et al. (2005) model. These distributions are used to compute the vertical structure at every annulus of radial distance. (This means the solutions of each annulus do not communicate with each other.) We assume discs are symmetric around the $z$-axis and $z=0$ plane. NSDs have two energy sources: viscosity and SF. The energy transport is solely by the radiation in the vertical direction. In other words, the conduction and convection heat transports are excluded in modeling. The total radiation flux is parameterized by an effective temperature ($T\textsubscript{\text{eff}}$) since the energy is transported losslessly to the surface and then radiated away.

3.1 Vertical Structure of an Annulus

One can ask what is the vertical structure for a given column of gas with total energy flux at a radial distance $R$ away from a central black hole of mass $M\textsubscript{\text{bh}}$. In order to obtain the hydrostatic solution for such a column of gas, the coupled one dimensional time independent radiative transfer equation (RTE), hydrostatic equation (HSE), and energy balance equation (EBE) must be solved.

The RTE governs the distribution of the radiation field within a material which is illustrated in Eq. 6 for a 1D slab:

[TABLE]

Here, $\nu$ is the frequency of light, $\mu$ is the cosine angle, $\eta_{\nu}$ is the emissivity, and $\chi_{\nu}$ is the extinction coefficient. To solve the RTE, we make the following approximations:

a homogeneous plane-parallel geometry (azimuthal symmetry) where properties of a material only depends on the $z$-direction, i.e. $\chi_{\nu}$ is isotropic, 2. 2.

isotropic absorption and emission where a source function $S_{\nu}$ is independent of $\mu$, $S_{\nu}(\mu,z)=S_{\nu}(z)$, 3. 3.

local thermodynamical equilibrium (LTE), 4. 4.

grey approximation (frequency integrated radiation field and opacity), and 5. 5.

the Eddington approximations: $f_{k}(\tau)\equiv\frac{K(\tau)}{J(\tau)}$ and  $f_{H}(\tau=0)\equiv\frac{H(\tau=0)}{J(\tau=0)}$ where $J$, $H$, and $K$ are the first, second, and third angular moments of a radiation field.

Then, the temperature profile computed from the RTE for an energy source is given by

[TABLE]

where $q(\tau)$ is the Hopf function111The form of $q(\tau)$ can be obtained from Eq. 3.8 in the paper by Hubeny (1990). and $T\textsubscript{\text{eff}}$ is the total effective temperature. The optical depth $\tau$ is defined as

[TABLE]

where $\kappa$ is the opacity and $\Sigma_{1}$ is the mass-integrated column density of the surface layer. For numerical convenience, $\Sigma$ (which increases monotonically from the surface layer of a slab to its mid-plane) is chosen as the independent variable for all the equations. Then, the vertical distance of a layer from the mid-plane, scale-height $z$, is related to $\Sigma$ with the following equation:

[TABLE]

where $\rho$ is the density. Since a NSD has two energy sources, an overall temperature can be computed by adding the fluxes of all the sources, $F_{\text{tot}}=F_{*}+F_{\text{vis}}$, which results in

[TABLE]

Here, $T_{*}$ and $T_{\text{vis}}$ are temperatures corresponding to the flux produced by viscosity and SF, respectively. Under LTE, $\kappa$ depends on $\rho$, $T$,

[TABLE]

and the composition of the slab. The calculation of $\kappa$ is described in Appendix A.

For the condition of hydrostatic balance, the second order differential form of the HSE,

[TABLE]

is used with the following boundary conditions222Subscript ‘1’ and ‘mp’ represent quantities at the surface and mid-plane, respectively.:

[TABLE]

Here, $P$ is the total pressure and $\Omega$ is the Keplerian frequency. $P_{1}$ is computed assuming a constant temperature for $\Sigma<\Sigma_{1}$ (Hubeny, 1990, section IV(b)). $c_{s}$ is the total speed of sound defined as:

[TABLE]

The HSE also includes an approximated self-gravity term (Paczynski, 1978) and it is solved in the logarithmic space in order to resolve the structure accurately near the surface. The solution of this boundary value problem can be achieved using the finite-difference method.

Finally, the EBE, Eq. 14, balances the total pressure with the gas pressure $P\textsubscript{\text{gas}}$, radiation pressure $P\textsubscript{\text{rad}}$, and the turbulence pressure $P\textsubscript{\text{turb}}$:

[TABLE]

The gas and radiation pressures are computed as

[TABLE]

and

[TABLE]

while $P\textsubscript{\text{turb}}$ is approximated (Thompson et al., 2005) as

[TABLE]

where $E_{0}=$ 1051 ergs is the energy released from a supernova. The turbulence speed $v\textsubscript{\text{turb}}$ is approximated to be a constant in the vertical direction using the density computed from the radial structure if $\dot{\Sigma}\textsubscript{*}$ is only allowed to occur at the mid-plane. $k_{B}$, $m_{p}$, $\sigma_{\text{sb}}$, and $c$ are the Boltzmann constant, mass of the proton, Stefan-Boltzmann constant, and the speed of light, respectively. We also require a column of gas to satisfy the Toomre stability criteria at the mid-plane which is parameterized by $Q$:

[TABLE]

Since the only known information about these coupled equations are the boundary conditions, an iterative method is used to solve them.

3.2 Iterative Method

To solve the coupled equations described in the previous section, an iterative algorithm is developed which is non-trivial due to the difficulties in achieving convergence. However, the major iterative steps for constructing the 2D structure of a NSD are as follow:

Using the Thompson et al. (2005) model, compute the radial distribution of mass-integrated column density at the mid-plane $\Sigma_{\text{mp}}(R)$ and the effective temperature $T\textsubscript{\text{eff}}(R)$. 2. 2.

Divide the radial distance (from $R\textsubscript{\text{in}}$ to $R\textsubscript{\text{out}}$) into 50 annuli which are spaced logarithmically. 3. 3.

For a given annulus, we set 300 vertical logarithmic grid points in terms of $\Sigma$ from the surface layer $\Sigma_{1}=$ 10-3 g cm*-2* to the mid-plane layer $\Sigma\textsubscript{\text{mp}}$. 4. 4.

Compute the initial solution (isothermal profiles) for the energy source embedded only at the mid-plane using a similar procedure to the one described in section IV(a) of Hubeny (1990). 5. 4.1.

Assume $\kappa(\Sigma_{1})=1.0$  cm2 g*-1* and $\rho_{\text{mp}}$ from the radial density profile. Using the analytical isothermal solution, compute $z_{1}$, $\tau_{1}$, and $T_{1}$ in subsequent order. Compute $\kappa_{1}$ for a given $\rho_{1}$ and $T_{1}$ using the procedure described in Appen. A; thereafter, compute the remaining physical quantities (i.e., pressures and speed of sound). 6. 4.2.

Repeat step 4.1 for the rest of the layers to complete the initial vertical profiles. If there is any numerical issue in computing the initial solution, it can be fixed by changing $\rho_{\text{mp}}$ or $\Sigma_{1}$. For example, the $\Sigma_{1}$ may be adjusted for a low density environment. 7. 5.

Update the pressure by solving Eq. 12, where $P\textsubscript{\text{mp}}$, $P_{1}$, and $c_{s}$ are used from a previous iteration. 8. 6.

Using the updated pressure (or the speed of sound), re-compute $\rho(\Sigma)$, $\kappa(\Sigma)$, $\tau(\Sigma)$, and $T(\Sigma)$ in subsequent order. 9. 7.

Since, the RTE and HSE are not solved simultaneously, the pressure computed from the HSE does not necessarily match with the EBE which induces a discrepancy in pressure, $\delta P\textsubscript{\text{E}}$. If Eq. 13 is not consistent with the HSE, it causes another discrepancy in pressure, $\delta P\textsubscript{\text{c}}$. To reduce $\delta P\textsubscript{\text{E}}$ and $\delta P\textsubscript{\text{c}}$, iterate the previous solution by updating $c_{s}$ using the induced $\delta P\textsubscript{\text{E}}$ as follows,

[TABLE]

Once $\delta P\textsubscript{\text{E}}$ and $\delta P\textsubscript{\text{c}}$ reach their minimum values, move to the next step. 10. 8.

If $\delta P\textsubscript{\text{E}}$ and $\delta P\textsubscript{\text{c}}$ are within the set discrepancy of 0.05, skip this step. Otherwise, go back to the step 5. 11. 9.

So far the achieved solution does not necessary satisfy the Toomre criteria ($Q=1$). To maintain this condition at the mid-plane, repeat steps 4-8, but change $P\textsubscript{\text{mp}}$ while solving the HSE. $P\textsubscript{\text{mp}}$ corresponding to $Q=$1 can be found using the Newton-Raphson method333$P\textsubscript{\text{mp}}=0.9P\textsubscript{\text{rad}}$ from the initial solution is a good initial guess for the Newton-Raphson method where $P\textsubscript{\text{rad}}$ is the radiation pressure at the mid-plane.. Finally, the self-consistent profiles are computed for a given annulus. 12. 10.

Repeat steps 1 to 9 for all the radii of a disc to construct the 2D profiles where $\theta$ is the viewing angle measured from the $z=0$ plane.

Using this developed numerical scheme, we explore NSDs under various conditions (input parameter space) that are shown in Table 2. The Mach number $m$ is chosen on the order of a tenth because the gas is not accreted to a parsec scale for a very low $m$ since all the gas is converted into stars in the outer region of the disc and the high Mach number (close to unity or higher) requires an inclusion of additional physics (i.e., shock waves). The size of the disc, the gas fraction at the outer radius, and the mass of the black hole are chosen such that they cover a reasonable physical range of NSDs. Taking all the possible combination of these four input parameters, we compute in total 192 models.

4 Results & Discussion

Out of 192 models, 99 discs show a large expansion of an atmosphere (comparable to $R$ or even higher) on parsec/sub-parsec scale. Here, we present the explicit profiles for one of these discs (hereafter Model A) as an example and its vertical profiles at the critical radius $R\textsubscript{\text{cric}}$ (where the largest expansion occurs). Later, we compute a distribution of $N_{H}$ based on a random selection from sets of $\theta$ and 192 models.

4.1 Model A

Model A is the nuclear starburst disc whose $f\textsubscript{\text{g,out}}$, $R\textsubscript{\text{out}}$, $m$, and $M\textsubscript{\text{bh}}$ are set to 80%, 180 pc, 0.5, and $10^{7.5}M_{\astrosun}$, respectively. As gas accretes toward the black hole, the gas is depleted into star-formation in order to maintain the stability criteria (Fig. 1 (left)). In addition, each column of gas becomes denser as gas gets closer to the black hole due to an increase in the vertical gravitational component. This increase in $\Sigma$ enhances the optical depth. As the gas becomes more optically thick, its temperature rises since $T$ is directly related to $\tau$ (which is the consequence of the RTE solution). When the gas crosses the dust sublimation temperature which are 1400K and 1750K for silicate and graphite grains, respectively (Kishimoto et al., 2013), starbursts occur near $\sim$4.0 pc and $\sim$0.884 pc as shown in Fig. 1 (left). The loss of the dust opacity induces a large opacity gradient. This opacity gradient causes a large expansion of the atmosphere which is illustrated in the right panel of Fig. 1. This atmosphere is supported by radiation pressure from starbursts (Fig. 1 (left)). The right panel of Fig. 1 compares the approximated surface scale-height (Eq. 5) from the Thompson et al. (2005) model to the explicit calculation. The dashed-red curve assumes that SF occurs at the mid-plane, while the dashed-dotted-green curve represents the model with a vertical distribution of SF (which is discussed further in detail in Sect. 4.5). The exact calculation predicts a lower scale-height ($\approx$ 4 orders of magnitude lower at $R\textsubscript{\text{cric}}=0.884$ pc) than the simple approximation except at the outer part of the disc. This discrepancy is due to the inaccuracy in assumptions made under Eq. 5: (1) the surface temperature is the same as $T\textsubscript{\text{eff}}$, (2) a constant $\kappa_{\text{mid}}$ from $z=0$ to $h_{\text{mid}}$, and (3) $\kappa$ increases linearly with $z$. The expansion of the atmosphere at $R\textsubscript{\text{cric}}$ is larger than the radial distance by a factor of 21.5 and covers 97% of the sky (observed from the central black hole).

Fig. 2 shows several important vertical profiles of a slab at $R\textsubscript{\text{cric}}=0.884$ pc for model A. As $\Sigma$ increases from the surface layer ($\Sigma=$10-2 g cm*-2* ) to the mid-plane ($\Sigma\textsubscript{\text{mp}}$), density and temperature also increase to achieve hydrostatic balance except when dust sublimates near $\Sigma=1.0$ g cm*-2* which is shown in Fig. 2 (a). The sublimation of grains results in the immediate drop in the opacity (solid-blue line in Fig. 2 (b)). This opacity gradient increases the total speed of sound (Fig. 2 (c)) which, in turn, expands the dusty atmosphere (dashed-green line in Fig. 2 (b)). Panel (d) shows that radiation pressure (dashed-red curve) dominates throughout the slab over the gas pressure (solid-black curve) and the turbulence pressure (dashed-dotted-blue curve). In such slabs, the density inversion is expected locally since grains are sublimated; however, its effect on the overall scale-height $h$ is insignificant since only a small fraction ($<1\%$) of $\Sigma\textsubscript{\text{mp}}$ is inverted. In addition, the gas pressure curve in panel (d) shows a slight inversion near the mid-plane. This is due to another density inversion caused by the boundary condition of our systems (outburst of energy flux at the mid-plane).

Once the 2D profiles of a disc is calculated successfully, a subsequent task is to compute the column density along a line of sight $N_{H}$ which is given by

[TABLE]

For model A, the distribution of $N_{H}$ as a function of $\theta$ is shown in Fig. 3 and shows that $N_{H}$ decreases with $\theta$. The peaks in $N_{H}$ are due to the density inversion phenomenon which occurs in slabs where dust is sublimated (Fig. 2a). The relative size of these peaks trace the offset of the large opacity gradient among these slabs. The figure illustrates that $N_{H}$ appears to be CK for $\theta\la$30°and also near $\theta$ where the peaks occur ($\approx$72°and 87°). For any other $\theta$ between $\approx$ 30°and $\theta_{\text{max}}=87.3$, a central AGN is obscured by the CN medium. For $\theta$>87.3°, the central AGN appears to be Type 1. A NSD therefore is consistent with the simple AGN unification model where the line-of-sight into the nucleus determines the observed obscuration properties. Although Model A would imply a Type 2/Type 1 ratio of $\sim$7:1, far higher than the observed ratio $\sim$4:1 (Gilli et al., 2007)), depending on the physical conditions $\theta_{\mathrm{max}}$ can be as low as 40°, giving a Type 2 to Type 1 ratio of 0.7. The distribution of the obscured to unobscured sky covering factors of all 99 NSD models with an inflationary atmosphere at parsec scale is shown in Fig. 4. The plot shows that, while large ratios are more common in our parameter space, Type 2/Type 1 ratios in line with the observations are found in 31% of the models. These are the models whose $h/R$ at the critical radius is less than 2 and their physical conditions are scattered across the input parameter space. Therefore, the fraction of sky obscured by the NSD can vary significantly from galaxy to galaxy.

4.2 All Models: scale-height $h$

A total of 99 out of 192 models (52%) show the inflation of an atmosphere at a critical radius which is illustrated in Fig. 5. It shows the dependency of $h/R$ on the black hole mass and the critical radius. The black, red, green and white data points at the parsec scale show that the critical radius gets closer to a black hole as $M\textsubscript{\text{bh}}$ decreases. This is expected since the gas has to accrete more inward in order to reach to the dust sublimation temperature for the lower $M\textsubscript{\text{bh}}$. Depending on the physical conditions, $h/R$ at the sub-parsec scale can range on the order of few tenths to tens. These large expanded vertical structure can obscure the incoming AGN light. Especially, since these atmospheres are dusty, the UV/optical photons will be reprocessed into IR photons. For a given black hole mass, the $R\textsubscript{\text{cric}}$ increases as a function of $f\textsubscript{\text{g,out}}$ since the optical depth reaches the sublimation temperature quicker due to the availability of gas. The amount of expansion at these critical radii are mainly governed by the magnitude of starbursts occurring at the mid-plane. Table  3 shows the distribution of 99 models across the input parameter space. In order to acquire an expanded atmosphere at the parsec scale, a smaller size disc, larger gas fraction, and higher Mach number are favorable conditions which is in the agreement with the work of Ballantyne (2008).

4.3 $N_{H}$ Distribution

To estimate the distribution of the column density observed along a line of sight ($N_{H}$) from a large sample of NSDs, the domain of the viewing angle $\theta\in$ [0°,90°] is divided into 30 bins. Afterward, the distribution of $N_{H}$ is computed by randomly selecting $N_{H}$ 10,000 times from the sets of 192 models and 30 bins of $\theta$. The top panel of Fig. 6 shows the histogram of the $N_{H}$ fraction in a given bin. Based on the random selection, the fraction of Type 1, CN Type 2s, and CK are 0.56, 0.23, and 0.21, respectively. For obscured AGNs, the fraction of AGNs peaks near $\log[N_{H}(\text{cm}\textsuperscript{-2})]\approx 23.5$. The bottom panel of Fig. 6 represents the scattering of $\theta$ for obscured AGNs when they are observed randomly. The solid data points correspond to CN AGNs and the starred ones are for CK AGNs. The data points with $\theta\la 20$°  near $N_{H}$=1022.5 cm-2  are gas dominated NSDs where starbursts did not occur. All the models appeared to be CK when they are viewed edge-on. This is expected since dense gas occupies the mid-plane regions of the discs. The figure illustrates that a given $N_{H}$ of CK or a heavily obscured CN AGN spans almost the entire domain of $\theta$ (from 0° to $\approx$80°) dependent on physical conditions. That means that, under certain conditions, the AGNs can be obscured by near face-on CK dusty gas.

The 2D NSD theory predicts $f\textsubscript{\text{CN}}$ and $f\textsubscript{\text{CK}}$ to be 23% and 21%, respectively, which is shown in the top panel of Fig. 6. This is in reasonable agreement with the values mentioned in the literature. For instance, the modeling of CXB suggests that $f\textsubscript{\text{CK}}$ can range from 5% to 50% (Akylas et al., 2012). Brightman & Ueda (2012) found $f\textsubscript{\text{CK}}$ $\approx 0.25$ around redshift 1 by studying the sample of Chandra Deep Field South survey. In the local universe, $f\textsubscript{\text{CK}}$ is $\approx 20\%$ based on the hard X-ray (Burlon et al., 2011), the optical (Akylas & Georgantopoulos, 2009), and IR (Brightman & Nandra, 2011a, b) samples. Ueda et al. (2014) constraints the ratio of CK to CN AGNs to be $\approx 0.5-1.6$ in order to produce the 20-50 keV part of the CXB spectrum which is in the agreement with the prediction of the NSD theory, $0.9$. The outcomes of our work is in good agreement with observations and suggests that NSDs could be potential sources of AGN obscuration at intermediate redshifts.

4.4 Inversion Phenomena

Many of our hydrostatic models of NSDs show a density inversion phenomenon. For example, the phenomenon is observed in model A whose vertical profiles at $R\textsubscript{\text{cric}}$ are shown in Fig. 2. There are two causes for the density inversion: (1) sublimation of grains and (2) demand of efficient energy flux transport. The latter case has been widely studied by theorists (Chitre & Shaviv, 1967; Joss et al., 1973; Schwarzschild, 1975; Tuchman et al., 1978; Harpaz, 1984). This inversion phenomenon can lead to a possible Rayleigh-Taylor instability (RTI). The development of RTI in radiation pressure dominated dusty environment has been also observed in hydrodynamical numerical simulations (e.g., Krumholz & Thompson, 2012; Davis et al., 2014). We expect the density inversion due to the sublimation of grains to have an insignificant effect on $N_{H}$ distributions (thus also $f\textsubscript{\text{CN}}$ and $f\textsubscript{\text{CK}}$) since the inverted column density of gas is a small fraction (< 1%) of the total column density.

The density inversion in gas pressure and turbulence pressure dominated layers causes a local pressure inversion. The development of a pressure inversion is expected in thermally unstable regions when temperature is monotonically increasing with $\tau$ (e.g., Rozanska & Czerny, 1996). Such a development of inverted pressure has been already reported in the modeling of astrophysical atmospheres in the past (Achmad et al., 1997; Asplund, 1998; Helling & Winters, 2001; Różańska et al., 2002). In order to analyze the effect of the inverted gas on a NSD, a time-dependent calculation becomes necessary which is beyond the scope of this paper. However, our comprehensive conclusion that NSDs are plausible sources for AGN obscuration is still absolute since the loss of dust opacity can cause an inflation to a parsec/sub-parsec atmosphere by orders of magnitude which is shown in the panel (b) of Fig. 2.

4.5 Vertical Distribution of $\dot{\Sigma}\textsubscript{*}$

Star-formation (SF) does not necessarily have to occur only at the mid-plane. Therefore, it would be interesting to study how the vertical distribution of SF affects the scale-height $h$ and $N_{H}$ distribution of NSDs. The Kennicutt-Schmidt law (which is also called the star-formation law) relates star-formation rate density $\dot{\Sigma}\textsubscript{*}$ to the surface gas density $\Sigma_{\text{gas}}$ by $\Sigma_{\text{gas}}\propto\dot{\Sigma}^{\alpha}_{*}$ (Schmidt, 1959; Kennicutt, 1998; Kennicutt & Evans, 2012), where $\alpha$ is expected to be $\approx 1$ for radiation pressure dominated regions (Ballantyne et al., 2013). Hence, we model the linear vertical distribution of $\dot{\Sigma}\textsubscript{*}$ on $\Sigma$ using the prescription presented by Hubeny & Hubeny (1998). Implementation of the vertical distribution assumes: (1) a surface scale-height at $\tau=1$ and (2) $\tau_{\text{mp}}$>1.0 (For optically thin column of gas, SF is only allowed to occur at the mid-plane.). With this description, we present the result for model A with the vertical distribution of $\dot{\Sigma}\textsubscript{*}$ (hereafter, model $\bar{\text{A}}$).

Fig. 1 (dashed-dotted-green line) shows the radial distribution of the surface scale-height $\bar{h}$ (for model $\bar{\text{A}}$) and Fig. 7 plots the ratio of $\bar{h}$ to $h$ (for model A) as a function of $R$. For $R\leq 5$ pc, the atmospheres are entirely dominated by radiation pressure. On sub-parsec scale ($<$1.5 pc), grains are sublimated and the opacity gradient induced due to the dust sublimation mainly governs the surface scale-height. Therefore, $h$ is approximately $\bar{h}$ in this region. For $1.5<R<5$ pc, the region is still radiation dominated, but the maximum temperature is below the sublimation temperature of dust. Here, model $\bar{\text{A}}$ possesses a lower density region than model A which results in a slightly higher $\bar{h}$ than $h$. For $R>5$ pc, the region is mostly turbulence pressure dominated for model A and gas pressure dominated for model $\bar{\text{A}}$. In these regions, the surface scale-height is primarily set by the pressure inversion phenomena. Model $\bar{\text{A}}$ introduces higher pressure inversion than model A. Increasing pressure inversion causes increase in density (Helling & Winters, 2001) which results in decreasing $\bar{h}$ in comparison to $h$. At the outer radii, (last three slabs), a column of gas is optically thin; therefore, $\bar{h}/h=1.0$. Another source of a difference in $\bar{h}$ and $h$ is a difference in treatment of turbulence pressure in both models. For model A, $P\textsubscript{\text{turb}}$ was approximated by assuming a constant $v\textsubscript{\text{turb}}$, while, for model $\bar{\text{A}}$, $P\textsubscript{\text{turb}}$ was computed exactly using Eq. 17. In summary, the vertical distribution of SF has a small effect on the surface scale-height. We also checked the $N_{H}$ distributions from models A and $\bar{A}$ and their difference is negligible. In general, one expect to have a minimum effect of $\dot{\Sigma}\textsubscript{*}(\Sigma)\propto\Sigma$ on the vertical structure. That is because a distribution of energy source only affects the Hopf function $q(\tau)$ which is introduced in the temperature solution of RTE and $T(\tau)$ has a weak dependance on the linear distribution of $\dot{\Sigma}\textsubscript{*}$.

4.6 Caveats & Directions for Future Work

One of the important tasks in modeling a physical system is to be aware of possible caveats which can be introduced due to assumptions, approximations, or an exclusion of relevant mechanisms. In the presented framework of a NSD, due to the assumptions of diffusion approximation and grey problem, results from modeling may be inaccurate for an optically thin limit ($\tau<<1$). However, this does not affect our conclusion that NSDs posses an inflationary atmosphere on parsec/sub-parsec scale which can potentially obscure the AGN irradiation since $h\sim R$ is achieved before $\tau$ reaches the optically thin limit. Another implicit assumption due to statistics based on a random selection is that there is not any internal/natural biased in the sample of input parameters. However, this assumption may be imprecise; for instance, the black hole mass distribution function is not evenly weighted for a given redshift. Another caveat in the modeling is that it misses outflows. In star-forming regions, outflows are possible to occur due to stellar feedback (Leitherer et al., 1992), supernovae feedback (Chevalier & Clegg, 1985), and radiation pressure on dusty atmosphere (Murray et al., 2005). Moreover, outflows are observed in many star-forming galaxies with the speed of $>1000$km s-1 at $z\sim 0.6$ (Tremonti et al., 2007; Diamond-Stanic et al., 2012). In the sample of local starbursts, outflows with the speed of 1500 km s-1 are observed by Heckman et al. (2011). In the future, we may study our sampled models using Bayesian statistics in order to remove an intrinsic biased and compare its conclusion with statistics based on a random drawings.

Another missing ingredient in the modeling is the effect of AGN (if they are present in galaxies) on the NSD structure. To study its importance on the presented modeling, we compute the optical depth along line of sights ($\tau_{\text{los}}$) for $\theta$ above the dust sublimation layer at the critical radius and check if $\tau_{\text{los}}$ is greater than the unity. Only 16 out of 99 models are optically thin and being transparent to the ultraviolet (UV) radiation. Depending on the physical conditions of the NSD, the UV energy can be deposited near the surface or in the deeper region along a line of sight since $\tau_{\text{los}}$ can be as large as 104. Moreover, we can compute the dust sublimation radius ($r_{\text{s}}$) with the following equation (Barvainis, 1987):

[TABLE]

Here, $L_{\text{UV}}$ is the UV luminosity of the AGN and $T_{s}$ is the sublimation temperature. Using the sublimation temperature of graphite grains (1750 K) and a typical bolometric luminosity (which gives the upper limit on $r_{s}$) range of Seyfert galaxies (1042 to 1044 ergs s*-1*), the approximated $r_{s}$ ranges from 0.008 to 0.08 pc, whereas Fig. 5 shows that the critical radius (where the largest expansion occurs in a disc) ranges from 0.1 to 2.3 pc. This implies that the AGN heating does not completely destroy the grains which are residing at the surface of the NSD at parsec scale. Based on this rough estimate, the AGN heating may not have a significant effect on the presented hydrostatic structure. However, a proper treatment of the AGN irradiation is required to be conclusive which will be done in detail in future work. Finally, NSDs may be present in galaxies with no central AGNs; therefore, this model may have wider applicability beyond AGNs.

Based on Hubble Space Telescope and ground based telescopes, nuclear star clusters (NSCs) are observed to be a common phenomena occurring in many types of galaxies: 70% in spheroidal galaxies (Côté et al., 2006) and 75% in late-type (Sc-Sd) disc galaxies (Böker et al., 2002; Walcher et al., 2005). Their sizes are at the parsec scale, 2-5 pc (Geha et al., 2002; Böker et al., 2004; Côté et al., 2006). Two possible scenarios are proposed for the formation of NSCs: (1) merger of dense clusters via dynamical friction (Andersen et al., 2008; Capuzzo-Dolcetta & Miocchi, 2008) and (2) gas accretion at the parsec-scale (Mastropietro et al., 2005; Seth et al., 2006). NSDs fall under the second scenario which can be a potential progenitor of NSCs. The NSD theory predicts a compact star-forming region with intense SF rate ($>3M_{\astrosun}$year-1) on parsec scale. Moreover, the 2D calculation of NSDs show these regions possess an environment whose scale-height is approximately larger by an order of magnitude than the radial distance. These characteristics of NSDs pose an interesting question: can nuclear star clusters be remnants of nuclear starburst discs? We will be addressing this question in the future investigation.

5 Conclusion

To summarize, we have successfully developed an iterative algorithm to compute the 2D hydrostatic structure of a nuclear starburst disc. These discs can be a potential source in obscuring the AGN radiation when an abundant amount of gas is available in galaxies at intermediate redshift fueling the central black hole. This modeling allows us to confirm that NSDs can possess the radiation pressure supported dusty structure on sub-parsec/parsec scale covering a major part of the sky (observed from a central back hole). Our results show that the NSD with fixed physical conditions can appear to be obscured by CN or CK gas depending on a viewing angle which supports the basic unification theory of AGNs. Below we summarize our main findings:

• •

A starburst phenomenon is more likely to occur in a disc with a smaller size, larger gas fraction, and higher Mach number.

• •

52% (99/192) of models show the large expansion of an atmosphere on parsec/sub-parsec scales. This indicates that a large expansion at the parsec/sub-parsec scale is a common phenomenon in dusty star-forming regions.

• •

Within the input parameter space, the atmosphere can be expanded from 0.2 to $\sim$ 30 $h/R$ at the critical radius ($R\textsubscript{\text{cric}}$) ranging from $\sim$0.2 to $\sim$2 pc. For a lower black hole mass, the inflation occurs closer to the black hole as indicated in Fig. 5.

• •

Based on the random sample of input parameters, the NSD theory predicts 56% of Type 1, 23% of CN Type 2s, and 21% of CK AGNs. Within the sample of obscured AGNs, the distribution of CN Type 2s peaks near 1023.5 cm-2. These predictions are consistent with observational evidence (e.g., Akylas & Georgantopoulos, 2009; Burlon et al., 2011; Brightman & Nandra, 2011a; Ueda et al., 2014).

• •

Based on the 2D NSD theory, $N_{H}$ along a line of sight varies from 1022.5 cm-2  to 1029.5 cm-2. From an edge-on view, all NSD models appear to be CK which is expected since dense gas resides at the mid-plane.

• •

Our results show a heavily obscured CN and CK AGNs span almost the entire range of $\theta$ from 0°  to $\approx$80°. The given AGN-NSD system (fixed input parameters) appears to be CN to CK for $\theta$ less than $\theta_{\text{max}}$; otherwise, the central AGN would appear to be Type 1. This supports the basic unification theory of AGNs. The restrictions on $N_{H,\text{max}}$ and $\theta_{\text{max}}$ of this system is determined by its physical conditions (input parameters).

• •

A linear vertical distribution of $\dot{\Sigma}\textsubscript{*}$ has an insignificant effect on $N_{H}$ distributions since the distribution of SF only affects the temperature solution of RTE (Eq. 7) and its effect is negligible.

Acknowledgement

This work is supported by NSF award AST 1333360. The authors are thankful to the referee for useful comments.

Appendix A Rosseland Mean Opacity

Modeling of radiation pressure dominated regions, such as NSDs, requires one to take into account the interaction between the radiation field and matter. This interaction can be quantified with a macroscopic quantity, opacity. A computation of opacity is not a trivial task since each chemical element interacts differently with a different frequency of light. Moreover, the inclusion of dust in the astrophysical composition increases the complexities due to the variability (e.g., size, shape, and distribution of grains) of dust opacity. However, for an optically thick medium (optical depth $\tau>>1$), the opacity of a medium can be well treated as the Rosseland mean opacity ($\kappa_{R}$) (Hubeny & Mihalas, 2014) where the radiation pressure can be well described by the diffusion approximation (Hubeny & Mihalas, 2014). Under the condition of LTE, $\kappa_{R}$ depends on the composition of matter, its temperature $T$, and the number density $n_{H}$. For NSDs, we choose Orion abundances to compute opacity under LTE using CLOUDY (Ferland et al., 2013). CLOUDY is a one dimensional photo-ionization code which includes grain physics (Weingartner & Draine, 2001; van Hoof et al., 2004; Weingartner et al., 2006) with a realistic size distribution of various grain species (e.g. Mathis et al., 1977); for instance, the size and shape based absorption and scattering processes of grains. For NSDs, we explore a wide range of temperature $T$ and density $n_{H}$:

[TABLE]

The curves in Fig. 8 show $\kappa_{R}$ as a function of $T$ for $n_{H}$ of 103 cm-3 and 106 cm-3. In low temperature regime, the $\kappa_{R}$ is increasing with the temperature for a given $n_{H}$. This is expected since dust is the dominant

source of opacity in this region and the dust cross-sectional area is proportional to the wavelength $\lambda$ of incoming electromagnetic wave by $\lambda^{-\alpha}$, where $\alpha$ approaches to 2 in the Rayleigh limit (Pollack et al., 1985). (As temperature increases, the peak of radiation power shifts toward a lower $\lambda$). The $\kappa_{R}$ has a small drop near 1400K (Kishimoto et al., 2013) due to the sublimation of silicate grains and the second deep drop occurs near 1750K (Kishimoto et al., 2013) due to the sublimation of graphite grains. The increase in $\kappa_{R}$ after the deep drop is due to the H- scattering opacity then subsequently bound-free and free-free opacity. Hydrogen starts ionizing near 5000K, and further ionized by $\sim 80\%$ and $\sim 100\%$ near 104K and 1.5$\times$104K, respectively. Using Saha’s equation, one can show an ionization fraction decreases (higher bound-free opacity) with density. This is responsible for an overall decreased in $\kappa_{R}$ with $n_{H}$ in the intermediate temperature regime since hydrogen is the abundant element in the composition. After hydrogen is fully ionized, the opacity is fairly independent of density since Thompson scattering becomes dominant. As temperature increases, gas absorbs more energy from incoming radiation which increases overall opacity of the medium.

Based on the complexities and behaviors of the Rosseland mean opacity as shown in Fig. 8, the two dimensional space of density and temperature is divided into five regions to have an optimal algorithm for interpolation. The logarithmic grids and its resolution for each region are described in Table 4 and were chosen such that it optimizes the computational time and accuracy of the interpolation algorithm. We choose a linear interpolation among logarithmic grid points. For region I, $\kappa_{R}$ is an independent of density as seen in Fig. 8. Therefore, the $\kappa_{R}$ of $\log[n_{H}(\text{cm}\textsuperscript{-3})]$ = 3.0 is taken for whole density range in interpolation function and only temperature grid points are interpolated.

To check the accuracy of the interpolation algorithm, opacities from the algorithm were predicted for various temperature and density (off the grids) points and the discrepancy against the real data from CLOUDY was computed as

[TABLE]

Examples of the error analysis is illustrated in Fig. 9. The percentage error in regions I and C are less than 4%. For regions A and B, the spacing of grid points is not linear because of the deep drop in $\kappa_{R}$ due to the sublimation of graphite grains. To minimize the percentage error in region B, various resolutions are chosen for subspaces. In region A, the interpolation function gives the percentage error that is mostly less than 10% except near $\log[T(\text{K})]\sim$ 3.3 and 3.6, where the error reaches 50% and 30%, respectively. Region B shows less than 4% error except it can reach up to 25% near 103.8K. In summary, every vertical slab may have few layers with an inaccurate opacity (yet less than 50% uncertainty). However, this will not have a significant effect on the overall structure for the purpose of this project since the number of these layers ($\approx$ 2-5) is a small fraction of total grid points (300+).