text:
# Toeplitz Determinants, Random Growth and Determinantal Processes

## 1. Introduction

Let $`\sigma `$ be a permutation from $`S_N`$. We say that $`\sigma (i_1),\mathrm{},\sigma (i_m)`$, $`i_1<\mathrm{}<i_m`$, is an increasing subsequence of $`\sigma `$ if $`\sigma (i_1)<\mathrm{}<\sigma (i_m)`$. The number $`m`$ is the length of the subsequence. The length of the longest increasing subsequence in $`\sigma `$ is denoted by $`\mathrm{}_N(\sigma )`$. If we pick $`\sigma `$ from $`S_N`$ uniformly at random $`\mathrm{}_N(\sigma )`$ becomes a random variable. Ulam’s problem, , is the study of the asymptotic properties as $`N\mathrm{}`$ of this random variable in particular its mean. It turns out that there is a surprisingly rich mathematical structure around this problem as we hope will be clear from the presentation below. It has been known for some time that $`𝔼[\mathrm{}_N]2\sqrt{N}`$ as $`N\mathrm{}`$, , . We refer to  for some background to the problem. A Poissonized version of the problem can be obtained by letting $`N`$ be an independent Poisson random variable with mean $`\alpha `$. This gives a random variable $`L(\alpha )`$ with distribution

$$[L(\alpha )n]=\underset{N=0}{\overset{\mathrm{}}{}}\frac{e^\alpha \alpha ^N}{N!}[\mathrm{}_Nn].$$
(1.1)
Since $`[\mathrm{}_Nn]`$ is a decreasing function of $`N`$, , asymptotics of the left hand side of (1.1) can be used to obtain asymptotics of $`[\mathrm{}_Nn]`$ (de-Poissonization).

The random variable $`L(\alpha )`$ can be realized geometrically using Hammersley’s picture, . Consider a Poisson process with intensity 1 in the square $`[0,\gamma ]^2`$, $`\gamma =\sqrt{\alpha }`$. An up/right path is a sequence of Poisson points $`(x_1,y_1),\mathrm{},(x_m,y_m)`$ in the square such that $`x_i<x_{i+1}`$ and $`y_i<y_{i+1}`$, $`i=1,\mathrm{},m1`$. The maximal number of points in an up/right path has the same distribution as $`L(\alpha )`$. A sequence of points realizing this maximum is called a maximal path. It is expected from heuristic arguments, see below, that the standard deviation of $`L(\alpha )`$ should be of order $`\gamma ^{1/3}=\alpha ^{1/6}`$. The proof that this is true, , and that we can also understand the law of the fluctuations is the main recent result that will be discussed below. Also, the deviations of a maximal path from the diagonal $`x=y`$ should be of order $`\gamma ^{2/3}`$. This last statement is proved in .

A generalization of the random variable $`L(\alpha )`$ can be defined in the following way. Let $`w(i,j)`$, $`(i,j)_+^2`$, be independent geometric random variables with parameter $`q`$. An up/right path $`\pi `$ from $`(1,1)`$ to $`(M,N)`$ is a sequence $`(1,1)=(i_1,j_1),(i_2,j_2),\mathrm{},(i_m,j_m)=(M,N)`$, $`m=M+N1`$, such that either $`i_{r+1}i_r=1`$ and $`j_{r+1}=j_r`$, or $`i_{r+1}=i_r`$ and $`j_{r+1}j_r=1`$. Set

$$G(M,N)=\underset{\pi }{\mathrm{max}}\underset{(i,j)\pi }{}w(i,j),$$
(1.2)
where the maximum is taken over all up/right paths $`\pi `$ from $`(1,1)`$ to $`(M,N)`$. Alternatively, we can define $`G(M,N)`$ recursively by

$$G(M,N)=\mathrm{max}(G(M1,N),G(M,N1))+w(M,N).$$
(1.3)
Some thought shows that if we let $`q=\alpha /N^2`$ then $`G(N,N)`$ converges in distribution to $`L(\alpha )`$ as $`N\mathrm{}`$, , so we can view $`G(N,N)`$ as a generalization of $`L(\alpha )`$. We can think of (1.2) as a directed last-passage site percolation problem. Since all paths $`\pi `$ have the same length, if $`w(i,j)`$ were a bounded random variable we could relate (1.2) to the corresponding first-passage site percolation problem, with a min instead of a max. in (1.2). The random variable $`G(M,N)`$ connects with many different problems, a corner growth model, zero-temperature directed polymers, totally asymmetric simple exclusion processes and domino tilings of the Aztec diamond, see ,  and references therein. It is also related to another growth model, the (discrete) polynuclear growth (PNG) model, ,  defined as follows. Let $`h(x,t)`$ denote the height above $`x`$ at time $`t`$. The growth model is defined by the recursion

$$h(x,t+1)=\mathrm{max}(h(x1,t),h(x,t),h(x+1,t))+a(x,t),$$
(1.4)
where $`a(x,t)`$, $`(x,t)\times `$, are independent random variables. If we assume that $`a(x,t)=0`$ whenever $`xt`$ is even, and that the distribution of $`a(x,t)`$ is geometric with parameter $`q`$, then setting $`w(i,j)=a(ij,i+j1)`$, we obtain $`G(i,j)=h(ij,i+j1)`$. The growth model (1.4) has some relation to the so called Kardar-Parisi-Zhang (KPZ) equation, , and is expected to fall within the so called KPZ-universality class. The exponents 1/3 and 2/3 discussed above are the conjectured exponents for $`1+1`$-dimensional growth models in this class.

## 2. Orthogonal polynomial ensembles

Consider a probability (density) on $`\mathrm{\Omega }^N`$, $`\mathrm{\Omega }=`$, $``$, $``$ or $`\{0,1,\mathrm{},M\}`$ of the form

$$u_N(x)=\frac{1}{Z_N}\mathrm{\Delta }_N(x)^2\underset{j=1}{\overset{M}{}}w(x_j),$$
(2.5)
where $`\mathrm{\Delta }_N(x)=_{1i<jN}(x_ix_j)`$ is the Vandermonde determinant, $`w(x)`$ is some non-negative weight function on $`\mathrm{\Omega }`$ and $`Z_N`$ is a normalization constant. We call such a probability an orthogonal polynomial ensemble. We can think of this as a finite point process on $`\mathrm{\Omega }`$. Let $`d\mu `$ be Lebesgue or counting measure on $`\mathrm{\Omega }`$ and let $`p_n(x)=\kappa _nx^n+\mathrm{}`$ be the normalized orthogonal polynomials with respect to the measure $`w(x)d\mu (x)`$ on $`\mathrm{\Omega }`$. The correlation functions $`\rho _{m,N}(x_1,\mathrm{},x_m)`$ of the point process are given by determinants, we have a so called determinantal point process, . In fact,

$`\rho _{m,N}(x_1,\mathrm{},x_m)`$ $`={\displaystyle \frac{N!}{(Nm)!}}{\displaystyle _{\mathrm{\Omega }^{Nm}}}u_N(x)𝑑\mu (x_{m+1})\mathrm{}𝑑\mu (x_N)`$ (2.6)
$`=det(K_N(x_i,x_j))_{1i,jm},`$
where the kernel $`K_N`$ is given by

$$K_N(x,y)=\frac{\kappa _{N1}}{\kappa _N}\frac{p_N(x)p_{N1}(y)p_{N1}(x)p_N(y)}{xy}(w(x)w(y))^{1/2}.$$
(2.7)
A computation shows that for bounded $`f:\mathrm{\Omega }`$,

$`𝔼[{\displaystyle \underset{j=1}{\overset{N}{}}}(1+f(x_j))]`$ $`={\displaystyle \underset{k=0}{\overset{N}{}}}{\displaystyle \frac{1}{k!}}{\displaystyle _{\mathrm{\Omega }^k}}{\displaystyle \underset{j=1}{\overset{k}{}}}f(x_j)det(K_N(x_i,x_j))_{1i,jk}d^k\mu (x)`$ (2.8)
$`=det(I+fK_N)_{L^2(\mathrm{\Omega },d\mu )},`$
where the last determinant is the Fredholm determinant of the integral operator on $`L^2(\mathrm{\Omega },d\mu )`$ with kernel $`f(x)K_N(x,y)`$. In particular, we can compute hole or gap probabilities, e.g. the probability of having no particle in an interval $`I\mathrm{\Omega }`$ by taking $`f=\chi _I`$, minus the characteristic function of the interval $`I`$. If $`x_{\text{max}}=\mathrm{max}x_j`$ denotes the position of the rightmost particle it follows that

$$[x_{\text{max}}a]=det(IK_N)_{L^2((a,\mathrm{}),d\mu )}.$$
(2.9)
As $`N\mathrm{}`$ we can obtain limiting determinantal processes on $``$ or $``$ with kernel $`K`$, i.e. the probability (density) of finding particles at $`x_1,\mathrm{},x_m`$ is given by $`det(K(x_i,x_m))_{1i,jm}`$. We will be interested in the limit process around the rightmost particle. This is typically given by the Airy kernel, $`x,y`$,

$$A(x,y)=\frac{\text{Ai }(x)\text{Ai }^{}(y)\text{Ai }^{}(x)\text{Ai }(y)}{xy}.$$
(2.10)
This limit, the Airy point process has a rightmost particle almost surely and its position has the distribution function

$$F_2(\xi )=det(IA)_{L^2(\xi ,\mathrm{})},$$
(2.11)
known as the Tracy-Widom distribution, . Hole probabilities as functions of the endpoints of the intervals satisfy systems of differential equations, , . For example, we have

$$F_2(\xi )=\mathrm{exp}(_\xi ^{\mathrm{}}(x\xi )u(x)^2𝑑x),$$
(2.12)
where $`u`$ solves the Painlevé-II equation $`u^{\prime \prime }=xu+2u^3`$ with boundary condition $`u(x)Ai(x)`$ as $`x\mathrm{}`$.

An example of a measure of the form (2.5) comes from the Gaussian Unitary Ensemble (GUE) of random matrices. The GUE is a Gaussian measure on the space $`_N^{N^2}`$ of all $`N\times N`$ Hermitian matrices. It is defined by $`d\mu _{\text{GUE},N}(M)=Z_N^1\mathrm{exp}(\text{tr }M^2)dM`$, where $`dM`$ is the Lebesgue measure on $`_N`$ and $`Z_N`$ is a normalization constant. The corresponding eigenvalue measure has the form (2.5) with $`w(x)=\mathrm{exp}(x^2)`$ and $`\mathrm{\Omega }=`$, . Hence the $`p_n`$:s are multiples of the ordinary Hermite polynomials. The largest eigenvalue $`x_{\mathrm{max}}`$ will lie around $`\sqrt{2N}`$. This is related to the fact that the largest zero of $`p_N`$ lies around $`\sqrt{2N}`$. The local asymptotics of $`p_N(x)\mathrm{exp}(x^2/2)`$ around this point, $`x=\sqrt{2N}+\xi /N^{1/6}\sqrt{2}`$, is given by the Airy function, $`\text{Ai }(\xi )`$. This asymptotics, some estimates, (2.9) and (2.11) give the following result,

$$_{\text{GUE},N}[\frac{\sqrt{2N}x_{\mathrm{max}}2N}{N^{1/3}}\xi ]F_2(\xi )$$
(2.13)
as $`N\mathrm{}`$.

## 3. Some theorems

The previous section may seem unrelated to the first but as the next theorems will show the problem of understanding the distribution of $`L(\alpha )`$ and $`G(M,N)`$ fits nicely into the machinery of sect. 2.

Theorem 3.1. . Take $`\mathrm{\Omega }=`$, $`MN`$ and $`w(x)=\left(\genfrac{}{}{0pt}{}{MN+x}{x}\right)`$ in (2.5). Then $`G(M,N)`$ is distributed exactly as $`x_{\mathrm{max}}`$.

The corresponding orthogonal polynomials are the Meixner polynomials, a classical family of discrete orthogonal polynomials, and we refer to the measure obtained as the Meixner ensemble. It is an example of a discrete orthogonal polynomial ensemble, . By computing the appropriate Airy asymptotics of the Meixner polynomials we can use (2.9) to prove the next theorem.

Theorem 3.2. . Let $`\gamma 1`$ be fixed and set $`\omega (\gamma ,q)=(1q)^1(1+\sqrt{q\gamma })^21`$ and $`\sigma (\gamma ,q)=(1q)^1(q/\gamma )^{1/6}(\sqrt{\gamma }+\sqrt{q})^{2/3}(1+\sqrt{q\gamma })^{2/3}`$. Then,

$$\underset{N\mathrm{}}{lim}[\frac{G([\gamma N],N)\omega (\gamma ,q)N}{\sigma (\gamma ,q)N^{1/3}}\xi ]=F_2(\xi ).$$
(3.14)
Thus $`G([\gamma N],N)`$ fluctuates like the largest eigenvalue of a GUE matrix.

As discussed above by setting $`q=\alpha /N^2`$ we can obtain $`L(\alpha )`$ as a limit of $`G(N,N)`$ as $`N\mathrm{}`$. By taking this limit in theorem 3.1, and using the fact that the measure has determinantal correlation functions, we see that $`L(\alpha )`$ behaves like the rightmost particle in a determinantal point process on $``$ given by the discrete Bessel kernel, , ,

$$B^\alpha (x,y)=\sqrt{\alpha }\frac{J_x(2\sqrt{\alpha })J_{y+1}(2\sqrt{\alpha })J_{x+1}(2\sqrt{\alpha })J_y(2\sqrt{\alpha })}{xy},$$
(3.15)
$`x,y`$. This gives

$$[L(\alpha )n]=det(IB^\alpha )_{\mathrm{}^2(\{n,n+1,\mathrm{}\})}.$$
(3.16)
Once we have this formula we see that all we need is the classical asymptotic formula $`\alpha ^{1/6}J_{2\sqrt{\alpha }+\xi \alpha ^{1/6}}(2\sqrt{\alpha }\text{Ai }(\xi )`$ as $`\alpha \mathrm{}`$ uniformly in compact intervals, and some estimates of the Bessel functions in order to get a limit theorem for $`L(\alpha )`$:

Theorem 3.3. . As $`\alpha \mathrm{}`$,

$$[\frac{L(\alpha )2\sqrt{\alpha }}{\alpha ^{1/6}}\xi ]F_2(\xi ).$$
(3.17)
Note the similarity with (2.13), just replace $`N`$ by $`\sqrt{\alpha }`$. This result was first proved in  by another method, see below. De-poissonizing we get a limit theorem for $`\mathrm{}_N(\pi )`$, see .

## 4. Rewriting Toeplitz determinants

The Toeplitz determinant of order $`n`$ with generating function $`fL^1(𝕋)`$ is defined by

$$D_n(f)=det(\widehat{f}_{ij})_{1i,jn},$$
(4.18)
where $`\widehat{f}_k=(2\pi )^1_\pi ^\pi f(e^{i\theta })e^{ik\theta }𝑑\theta `$ are the complex Fourier coefficients of $`f`$. Consider the generating function

$$f(z)=\underset{\mathrm{}=1}{\overset{M}{}}(1+\frac{a_{\mathrm{}}}{z})(1+b_{\mathrm{}}z),$$
(4.19)
where $`a_{\mathrm{}},b_{\mathrm{}}`$ are complex numbers. The elementary symmetric polynomial $`e_m(a)`$, $`a=(a_1,\mathrm{},a_M)`$ is defined by $`_{j=1}^M(1+a_jz)=_{|m|<\mathrm{}}e_m(a)z^m`$. A straightforward computation shows that when $`f`$ is given by (4.19) then

$$\widehat{f}_{ij}=\underset{m=0}{\overset{\mathrm{}}{}}e_{mj}(a)e_{mk}(b).$$

Insert this into the definition (4.18) and use the Heine identity,

$`{\displaystyle \frac{1}{n!}}{\displaystyle _{\mathrm{\Omega }^n}}`$ $`det(\varphi _i(x_j))_{1i,jn}det(\psi _i(x_j))_{1i,jn}d^n\mu (x)`$ (4.20)
$`=det({\displaystyle _\mathrm{\Omega }}\varphi _i(x)\psi _j(x)𝑑\mu (x))_{1i,jn},`$
to see that

$$D_n(f)=\underset{m_1>m_2>\mathrm{}>m_n0}{}det(e_{m_ij}(a))_{1i,jn}det(e_{m_ij}(b))_{1i,jn}.$$
(4.21)
Here we have removed the $`n!`$ by ordering the variables. These determinants are again symmetric polynomials, the so called Schur polynomials. Let $`\lambda =(\lambda _1,\lambda _2,\mathrm{})`$ be a partition and let $`\lambda ^{}=(\lambda _1^{},\lambda _2^{},\mathrm{})`$ be the conjugate partition, . Set $`m_i=\lambda _i^{}+ni`$, $`i=1,\mathrm{},n`$ and $`\lambda _i^{}=0`$ if $`i>n`$, so that $`\lambda ^{}`$ has at most $`n`$ parts, $`\mathrm{}(\lambda ^{})n`$, which means that $`\lambda _1n`$. Then the Schur polynomial $`s_\lambda (a)`$ is given by

$$s_\lambda (a)=det(e_{\lambda _i^{}i+j}(a))_{1i,jn}=det(e_{m_ij}(a))_{1i,jn},$$
(4.22)
the Jacobi-Trudi identity. Hence,

$$D_n(f)=\underset{\lambda ;\lambda _1n}{}s_\lambda (a)s_\lambda (b),$$
(4.23)
and we have derived Gessel’s formula, . If we let $`n\mathrm{}`$ in the right hand side we obtain $`_{i,j=1}^M(1a_ib_j)^1`$ by the Cauchy identity, . In the case when all $`a_i,b_j[0,1]`$, $`s_\lambda (a)s_\lambda (b)0`$, and we can think of

$$\underset{i,j=1}{\overset{M}{}}(1a_ib_j)s_\lambda (a)s_\lambda (b)$$
(4.24)
as a probability measure on all partitions $`\lambda `$ with at most $`n`$ parts, the Schur measure, . In this formula we can insert the combinatorial definition of the Schur polynomial, ,

$$s_\lambda (a)=\underset{T:\text{sh }(T)=\lambda }{}a_1^{m_1(T)}\mathrm{}a_M^{m_M(T)},$$
(4.25)
where the sum is over all semi-standard Young tableaux $`T`$, , with shape $`\lambda `$, and $`m_i(T)`$ is the number of $`i`$:s in $`T`$.

A connection with the random variables in section 1 is now provided by the Robinson-Schensted-Knuth (RSK) correspondence, . This correspondence maps an $`M\times M`$ integer matrix to a pair of semi-standard Young tableaux $`(T,S)`$ with entries from $`\{1,2,\mathrm{},M\}`$. If we let the random variables $`w(i,j)`$ be independent geometric with parameter $`a_ib_j`$ then the RSK-correspondence maps the measure we get on the integer matrix $`(w(i,j))_{1i,jM}`$ to the Schur measure (4.24). Also, the RSK-correspondence is such that $`G(M,M)=\lambda _1`$, the length of the first row. If we put $`a_j=0`$ for $`N<jM`$ and $`a_i=b_i=\sqrt{q}`$ for $`1iN`$ in the Schur measure and set $`x_j=\lambda _j+Nj`$, $`1jN`$, we obtain the result in Theorem 3.1.

In the limit $`M=N\mathrm{}`$, $`q=\alpha /N^2`$, in which case $`G(N,N)`$ converges to $`L(\alpha )`$, the Schur measure converges to the so called Plancherel measure on partitions, , . In the variables $`\lambda _ii`$ this measure is a determinantal point process on $``$ given by the kernel $`B^\alpha `$, (3.15). This result was obtained independently in , which also gives a descrption in terms of different coordinates. See also  for a more direct geometric relation between GUE and the Plancherel measure. In this limit the Toeplitz determinant formula (4.23) gives

$$[L(\alpha )n]=e^\alpha D_n(e^{2\sqrt{\alpha }\mathrm{cos}\theta }).$$
(4.26)
This variant of Gessel’s formula was the starting point for the original proof of Theorem 3.3 in . The right hand side of (4.26) can be expressed in terms of the leading coefficients of the orthogonal polynomials on $`𝕋`$ with respect to the weight $`\mathrm{exp}(2\sqrt{\alpha }\mathrm{cos}\theta )`$. These orthogonal polynomials in turn can be obtained as a solution to a matrix-valued Riemann-Hilbert problem (RHP), and the asymptotics of this RHP as $`\alpha \mathrm{}`$ can be analyzes using the powerful asymptotic techniques developed by Deift and Zhou, . This approach leads to the formula (2.12) for the limiting distribution.

Write $`f=\mathrm{exp}(g)`$ and insert the definition of the Fourier coefficients into the definition (4.18). By the Heine identity we obtain an integral formula for the Toeplitz determinant,

$`D_n(f)`$ $`={\displaystyle \frac{1}{(2\pi )^nn!}}{\displaystyle _{[\pi ,\pi ]^n}}{\displaystyle \underset{1\mu <\nu n}{}}|e^{i\theta _\mu }e^{i\theta _\nu }|^2{\displaystyle \underset{\mu =1}{\overset{n}{}}}e^{g(e^{i\theta _\mu })}d^n\theta `$ (4.27)
$`={\displaystyle _{U(n)}}e^{\text{tr }g(U)}𝑑U.`$
In the last integral $`dU`$ denotes normalized Haar measure on the unitary group $`U(n)`$ and the identity is the Weyl integration formula. The limit of (4.26) as $`\alpha \mathrm{}`$ is then a so called double scaling limit in a unitary matrix model, . The formula (4.26) can also be obtained by considering the integral over the unitary group, see .

Another way to obtain the Schur measure is via families of non-intersecting paths which result from a multi-layer PNG model, . The determinants in the measure then come from the Karlin-McGregor theorem or the Lindström-Gessel-Viennot method.

## 5. A curiosity

Non-intersecting paths can also be used to describe certain tilings, e.g. domino tilings and tilings of a hexagon by rhombi. By looking at intersections with appropriate lines one can obtain discrete orthogonal polynomial ensembles. In the case of tilings of a hexagon by rhombi, which correspond to boxed planar partitions, , the Hahn ensemble, i.e. (2.5) with $`\mathrm{\Omega }=\{0,\mathrm{},M\}`$ and a weight giving the Hahn polynomials, is obtained, . The computation leading to this result also gives a proof of the classical MacMahon formula, , for the number of boxed planar partitions in an $`abc`$ cube, i.e. the number of rhombus tilings of an $`abc`$-hexagon. In terms of Schur polynomials the result is

$$\underset{\mu ;\mathrm{}(\mu )c}{}s_\mu ^{}(1^a)s_\mu ^{}(1^b)=\underset{i=1}{\overset{a}{}}\underset{j=1}{\overset{b}{}}\underset{k=1}{\overset{c}{}}\frac{i+j+k1}{i+j+k2},$$
(5.28)
where the right hand side is MacMahon’s formula. (Here $`1^a`$ means $`(1,\mathrm{},1)`$ with $`a`$ components.) Comparing this formula with the formula (4.23) we find

$$(1)^{a+b}D_n(\underset{\mathrm{}=1}{\overset{a}{}}(1e^{i\theta })\underset{\mathrm{}=1}{\overset{b}{}}(1e^{i\theta }))=\underset{i=1}{\overset{a}{}}\underset{j=1}{\overset{b}{}}\underset{k=1}{\overset{c}{}}\frac{i+j+k1}{i+j+k2}.$$
(5.29)
It has been conjectured by Keating and Snaith, , that the following result should hold for the moments of Riemann’s $`\zeta `$-function on the critical line,

$$\underset{T\mathrm{}}{lim}\frac{1}{(\mathrm{log}T)^{k^2}}\frac{1}{T}_0^T|\zeta (1/2+it)|^{2k}𝑑t=f_{\text{CUE}}(k)a(k),$$
(5.30)
where $`a(k)`$ is a constant depending on the primes,

$$f_{\text{CUE}}(k)=\underset{n\mathrm{}}{lim}\frac{1}{n^{k^2}}_{U(n)}|Z(U,\theta )|^{2k}𝑑U=\underset{j=0}{\overset{k1}{}}\frac{j!}{(j+k)!}$$
(5.31)
and $`Z(U,\theta )=det(IUe^{i\theta })`$ is the characteristic polynomial of the unitary matrix $`U`$. If we take $`a=b=k`$ in (5.29) and use (4.27) we find

$$_{U(n)}|Z(U,\theta )|^{2k}𝑑U=\underset{i=1}{\overset{k}{}}\underset{j=1}{\overset{k}{}}\underset{\mathrm{}=1}{\overset{n}{}}\frac{i+j+\mathrm{}1}{i+j+\mathrm{}2}=\underset{j=0}{\overset{n1}{}}\frac{j!(j+2k)!}{(j+k)!^2}$$
(5.32)
as computed in  by different methods. Letting $`n\mathrm{}`$ we obtain the last expression in (5.31). Hence, we see that the formula (5.32) has a curious combinatorial interpretation via MacMahon’s formula.


format: text

text:
# Revisiting the location and environment of the central engine in NGC 1068

## 1 Introduction

Efforts in understanding the physics at work in active galactic nuclei (AGN) started in the early seventies, long after the first AGN or quasar had been discovered (Seyfert, 1943; Schmidt, 1959). Great progress has been made over the last decade, thanks to improved observing tools and modeling explorations. On the side of observations, one can acknowledge the impact of better spatial and spectral resolutions in the near-infrared (NIR), mid-infrared (MIR) and radio wavelengths, as well as the access to high energy photons (UV, X-rays). On the side of modeling, the availability of complex codes taking into account the large number of physical processes taking place in AGN, and their inter-relations, has permitted to bridge in a more realistic fashion observations and model predictions (see for a review Krolik, 1999).

It is now commonly accepted that the main source of energy in AGN is of gravitational origin. The various components building up an AGN and its surroundings have also been identified. However, their geometrical arrangement and their links remain to be understood in more detail, with the hope that this will bring also some clues about the formation process of such systems. The commonly accepted AGN model involves a central engine (black hole plus accretion disc, hereafter referred to as CE) surrounded by sets of dense matter clouds (the so-called broad-line and narrow-line regions, hereafter BLR and NLR) and by a dusty/molecular disc-like and thick structure (hereafter the torus) which funnels the emission of high energy photons and particles along privileged directions (the ionizing cone, the radio jet axis). This in turn gives rise to viewing-angle dependent effects which are often invoked to explain the various brands of AGN discovered so far. And last, the AGN, including its close environment, are embedded in the bulge stellar population of the AGN host galaxy.

The intrinsic physical scales of the different components mentioned above are quite different, with orders of magnitude such as: (a) CE (black hole plus accretion disc), $`10^{14}`$ cm, (b) BLR, $`10^{16}`$ cm, (c) NLR, $`10^{20}`$ cm, (d) torus, $`10^{18}`$ to $`10^{20}`$ cm, while the stellar system core hosting the AGN has a typical dimension of $`10^{20}`$ cm (but can be more extended). On the other hand, the radio jet and structures can show a huge extension, in some cases on kpc scales. The different AGN components and the components contributing along the line-of-sight close to the AGN (such as the stellar population) can be disentangled in several ways: through their relative flux contributions which vary largely among AGN, possibly reflecting the AGN evolutionary state, and through their different angular sizes, according to the intrinsic sizes mentioned above. However, even in the close-by AGN NGC 1068, which lies at a distance of 14.4 Mpc, a torus component with intrinsic size $`10^{19}`$ cm corresponds to an apparent angular size of 0.04″. So far, this is out of reach from an observational point of view, except in the radio range with interferometers, or using speckle techniques.

Therefore, the comparison between observations and model predictions, necessary to probe model parameters, is limited by the angular resolution of available data sets. Conversely, one can focus on physical signatures specific to each component, some of which are listed hereafter. High energy photons and particles arise essentially from the CE (physical processes related to matter accretion and emission at the surface of the accretion disc), while the radio maser emission is thought to be related to the inner walls of the dusty/molecular torus. The flux variations of the CE provide information about the size and structure of the accretion disc. The impact of flux variations on the ionization state of the surrounding BLR clouds makes the BLR line emission respond with a time-delay which reflects the BLR size. The molecular line emission in the NIR to millimeter range, as well as the MIR continuum emission come predominantly from the dusty/molecular component (torus or/and cold material in BLR/NLR clouds shielded from the intense UV radiation field). The contribution from the stellar population starts to fade off at wavelengths above 3.5$`\mu \mathrm{m}`$. In conclusion, a whole set of constraints can be used and this is why the multi-wavelength approach is a very powerful tool in AGN studies. Still, to assemble these pieces of information, we have to face the fact that they have been obtained at different angular resolutions and that they are most often lacking precise astrometric information. These are major sources of limitation which have indeed hampered so far a quantitative probe of model predictions.

The preferred target on which one could solve the issue is of course the close-by AGN, NGC 1068. Because it is a type 2 AGN, it contains all the components mentioned above; because of its proximity it offers two advantages, its brightness and one of the best intrinsic spatial resolution one can achieve at every wavelength. Thanks to these features, a very large number of observational studies have been performed on NGC 1068 since the early seventies (ADS returns 200 entries for refereed articles containing NGC 1068 in their title, from 1970 to 2002). It is therefore tempting to try assembling major pieces of information obtained so far on this object, with the goal of locating its CE and analyzing the AGN surroundings.

The first task is to register, on the basis of their absolute coordinates, high resolution maps available at different wavelengths and locate the CE with respect to other AGN components. This is performed in Sect. 2, where we present a new determination of the absolute positioning of the 2.2$`\mu \mathrm{m}`$ continuum emission peak using ISAAC/VLT data in the H2 line emission and PdB/IRAM interferometric measurements in the CO(2-1) line and where we compile and discuss the few data sets with astrometric positioning available so far. After having secured the best possible registration, we compare maps obtained at different wavelengths and on different scales, in order to extract quantitative information to be compared with model predictions. In Sect. 3, we discuss the spectral energy distribution (SED) of the unresolved AGN component (intrinsic size less than 8pc FWHM), while the origin of the MIR dust emission observed within the ionizing cone is investigated in Sect. 6. Finally, we discuss and model the excitation of the molecular material located in the plane perpendicular to the ionizing cone axis (the torus plane) in Sect. 5. Concluding remarks and prospective ideas are presented in Sect. 6.

## 2 Relative positioning of the emission peaks at different wavelengths, location of the CE

The precise location of the CE in the AGN of NGC 1068 is still subject to a hot debate. We review hereafter a number of published astrometric measurements at different wavelengths. Then, we present a new absolute positioning of the unresolved K-band continuum peak and we discuss the location of emission peaks at different wavelengths with respect to the CE.

### 2.1 The central engine

It is commonly admitted that radio observations can directly probe the CE. The 6 cm MERLIN map of the AGN shows a linear radio structure made of four compact components (Muxlow et al., 1996; Gallimore et al., 1996c, a). One of these components (source S1) is now identified with little doubt as the CE, because of its thermal inverted spectrum between 1.3 cm and 6 cm (Gallimore et al., 1996a). Moreover, a disc of H2O and OH maser emission, thought to be associated with the inner walls of a molecular torus around the CE, is observed with the VLBA at the same location (Gallimore et al., 1996b). The absolute astrometric precision of the MERLIN images is 20 mas. Muxlow et al. (1996) report the position of S1 at $`\alpha =02^\mathrm{h}42^\mathrm{m}40.7098^\mathrm{s},\delta =00^{}00\mathrm{}47.938\mathrm{}`$ (J2000). In the following, we consider that S1 is the CE.

### 2.2 Astrometric positionings

The objective way to register any map to the 6 cm MERLIN image is by performing absolute astrometry. Any other method is subjective in the sense that it will require some physical assumption.

Even though 0.1″ resolution HST UV-continuum and optical NLR images have been published before 1997 (Evans et al., 1991; Macchetto et al., 1994), they could not be precisely registered since the absolute positioning of the HST images was only possible with a precision of $`\pm `$0.5″. Based on a method developed by Lattanzi et al. (1997), Capetti et al. (1997) report for the first time on an absolute registration of HST images with an $`\pm `$0.08″ precision. They find that the optical (F547M filter) continuum peak (hereafter referred to as HST F547M continuum) is located at $`\alpha =02^\mathrm{h}42^\mathrm{m}40.711^\mathrm{s},\delta =00^{}00\mathrm{}47.81\mathrm{}`$ (J2000), that is 0.13″$`\pm `$0.1″ North and 0.02″$`\pm `$0.1″ East of S1. This corresponds to the peak of the NLR-cloud B. They also show that the UV (F218W) and F547M peaks are coincident within 25 mas.

On the photographic 103a-O plates they use to perform the absolute astrometry, Capetti et al. (1997) locate the visible photocenter at $`\alpha =02^\mathrm{h}42^\mathrm{m}40.727^\mathrm{s},\delta =00^{}00\mathrm{}47.69\mathrm{}`$ (J2000). The difference between the locations of the HST F547M continuum peak and the photographic plate visible photocenter is 0.27″. Such a large difference can arise from two effects. First the F547M filter cuts out the prominent \[OIII\] $`\lambda `$5007 emission lines, while the ground-based 103a-O photographic image does sense these lines, the photocenter of which may be slightly displaced from the continuum photocenter. However, such a displacement of the \[OIII\] peak with respect to the F547M peak only occurs at the scale of individual NLR clouds. Indeed Macchetto et al. (1994) show that the UV continuum emission peak in NLR-cloud B is shifted by 0.05″ towards the North with respect to the \[OIII\] peak. Actually, the main reason for the observed shift of 0.27″ is that the photographic image, which has a low spatial resolution (larger than 1″), mixes the different continuum and \[OIII\] emitting clouds. Indeed, one observes that the photographic plate photocenter is shifted towards the North-East with respect to the F547M peak. This is consistent with the overall shape of the NLR, which is extended towards the North-East. Considering the complexity of the NLR visible continuum and \[OIII\] emission, images at different spatial resolutions ought to produce emission peaks at different locations. This effect must be considered carefully when interpreting the relative astrometric measurements between the emission peaks in the optical and in the IR, since the bulk of the NIR emission is unresolved down to 0.03″ (Wittkowski et al., 1998), while the optical emission is made of HST-resolved clouds over a region of $``$0.5″ $`\times `$1″. Fig. 1 illustrates this effect by showing the 0.1″ HST \[OIII\] image taken from Macchetto et al. (1994), and the same image after convolution with 0.5″ and 1.0″ FWHM Gaussian profiles, simulating seeing effects. The shifts of the position of the \[OIII\] peak on the images respectively degraded by an 0.5″ and 1.0″ seeing effect are 0.35″ and 0.45″, with respect to the HST original image. This is of the same order of magnitude as the measured offset between the HST F547M peak and the 103a-O visible photocenter, sensing the continuum (and \[OIII\] lines for the 103a-O plates). As a consequence, the results of Braatz et al. (1993) and Marco et al. (1997) which position the N-band peak and K-band peak at 0.3″ South, 0.1″ West from respectively the optical peak (from an optical TV camera) and the I-band peak should be considered with caution: they both positioned the MIR or the NIR peaks with respect to the visible peak of a low resolution image. In order to gain in precision, we discuss in the following the positioning of the NIR peak with respect to the optical continuum peak of NLR-cloud B only.

### 2.3 Assumption-depending positionings

Positioning of the emission peaks at different wavelengths with respect to the CE have also been performed under various physical assumptions.

#### 2.3.1 Positioning of the K-band peak

The absolute astrometry of the K-band emission peak can be performed directly in two ways.

The first method is to compare the high resolution CO emission line map, for which absolute astrometry is available, and the H2 2.12$`\mu \mathrm{m}`$ emission line map, which appears to be very similar to the CO map and for which there is a relative position of the K-band peak. Tacconi et al. (1997) and Schinnerer et al. (2000) observed at $``$0.5″ resolution the CO(2-1) line emission in the millimeter wavelength range with the IRAM interferometer. They produced a CO map of the central 3″ region, finding that the CO emission is asymmetrically distributed, with a strong CO emission knot to the East of the CE. The CO emission peak corresponding to this knot is located at $`\alpha =2^\mathrm{h}42^\mathrm{m}40.767^\mathrm{s}`$ and $`\delta =0^{}00\mathrm{}47.95\mathrm{}`$ (J2000), with a precision of $`\pm `$0.05″. In parallel, Galliano & Alloin (2002) produced a map at similar resolution of the H2 2.12$`\mu \mathrm{m}`$ emission derived from VLT/ISAAC observations. The CO and H2 maps correlate very well (Fig. 4,e). Measurements on the data presented in Galliano & Alloin (2002) show that the eastern H2 knot is located 0.95″$`\pm 0.15\mathrm{}`$ East and 0.1″ $`\pm `$0.15″ North of the K-band continuum peak. The similarity of the two maps and the fact that, in the AGN context, CO and H2 emission are expected to be coincident (Maloney et al., 1996), allow us to assume that the CO and H2 peaks are at the same location within 0.1″. The absolute coordinates of the K-band peak can therefore be derived: $`\alpha =2^\mathrm{h}42^\mathrm{m}40.702^\mathrm{s}\pm 0.012^\mathrm{s}`$ and $`\delta =00^{}00\mathrm{}48.04\mathrm{}\pm 0.18\mathrm{}`$ (J2000). This position of the K-band emission peak is consistent with the position of S1, while it is only marginally consistent with the positions of the HST F547M continuum and the \[OIII\] peaks.

Another positioning of the K-band peak has been performed by Thompson et al. (2001). Using the similarity of the NLR features seen in the HST/WFPC2 Hα image and in the HST/NICMOS Pα image and using the fact that both Hα and Pα are Hydrogen recombination lines, and consequently should trace the same features, they have registered the position of the NIR continuum image with respect to the HST F547M continuum image. They find the location of the K-band peak to be consistent with both the center of the UV polarization pattern as determined by Kishimoto (1999, see Sect. 2.3), and the \[OIII\] peak (NLR-cloud B). This coincidence is claimed to be accurate to $`\pm 0.04\mathrm{}`$.

#### 2.3.2 Positioning of the center of the polarization patterns in the UV and in the IR

Capetti et al. (1995a, b) have presented HST imaging polarimetry of the NLR in the UV (F218W) and optical (F555W). The polarization vectors are found to be distributed in a circular pattern, as expected in the case of scattered radiation from a point-like source, assumed to be the CE. After a correction of their first result (Capetti et al., 1995a), they claimed to have determined the position of the center of polarization, at the location 0.45″ $`\pm `$0.015″ South, 0.06$`\pm `$0.015″ West of the optical (F555W) continuum peak. However, Kishimoto (1999) noticed that in the map of polarization position angles (PA) of Capetti et al. (1995a, b), clear deviations from the centrosymmetric pattern can be seen, although not discussed in that paper. Therefore, Kishimoto (1999) has extended this work, performing a very complete error budget, and has redetermined the location of the symmetry center of the polarization pattern in the UV within a convincing error circle. He finds a location for the CE which is significantly different from that of Capetti et al. (1995b). The revised location is at $`\alpha =\mathrm{\hspace{0.17em}02}^\mathrm{h}42^\mathrm{m}40.713^\mathrm{s}\pm 0.0006^\mathrm{s}`$ and $`\delta =00^{}00\mathrm{}47.93\mathrm{}_{0.11\mathrm{}}^{+0.14\mathrm{}}`$ (J2000).

Lumsden et al. (1999) have performed NIR and MIR imaging polarimetry at good spatial resolution ($``$0.5″ in the NIR). They detect NIR polarized emission in extended regions along PA$``$30° both towards the North and the South. In the optical the light from the southern part is absorbed, probably by matter in the plane of the host galaxy. They suspect that an additional mechanism, other than scattering, is contributing to the polarization at all wavebands close to the AGN. Therefore, they have masked the central 2″ region and derived the scattering center from the outer regions only. They find offsets between the polarization pattern center and the flux centroid in the three NIR bands (J,H,K) of respectively 0.11″ $`\pm `$0.18″ , 0.09″ $`\pm `$0.21″ and 0.34″ $`\pm `$0.39″ , offsets which are formally consistent with no offset. Unfortunately, they do not provide the PA of the offsets. If the PAs for the three offsets are different, then their measurements will actually be consistent with no offset, while on the contrary, if the three offsets were for example at PA$``$30°, then their result would be consistent with the UV polarization center found by Kishimoto (1999).

Recently, Simpson et al. (2002) have published NIR (1.9-2.1$`\mu \mathrm{m}`$) polarimetry from HST/NICMOS data set. They claim that the center of the NIR polarization pattern does correspond to the NIR peak within $`\pm 0.08`$″.

### 2.4 Summary and conclusion

We have compiled measured positions of the emission peaks at different wavelengths, of the centers of the polarization patterns in the UV and in the IR and of the positions of the radio sources in the region surrounding the AGN in NGC 1068. We find that the measurements converge to a unique picture, and that all the components are relatively positioned with a good precision. A sketch of the positions of the most important features is displayed in Fig. 2 while their absolute coordinates are listed in Table LABEL:table\_positions.

Thompson et al. (2001) found that the K-band peak is coincident with the position of the UV polarization pattern center and the peak of \[OIII\]$`\lambda `$5007 in the NLR-cloud B, at $`\pm 0.04\mathrm{}`$ precision. Simpson et al. (2002) found that, within $`\pm `$0.08″, the NIR (1.9-2.1$`\mu \mathrm{m}`$) polarization pattern is centered on the NIR peak. These three measurements converge to position the NIR peak at the position of S1 within $`\pm `$0.1″. We do not contemplate the offsets between visible peaks and N- and K-band peaks measured respectively by Braatz et al. (1993) and Marco et al. (1997), since these offsets are affected by the use of images at different spatial resolutions (Sect. 2.2). Regarding the UV polarization pattern center, we consider only the analysis of Kishimoto (1999), which is the most precise. There is no absolute positioning of the MIR maps available. A very reasonable assumption, supported by dusty torus models (Granato & Danese, 1994; Granato et al., 1997), is to consider that the MIR emission peak appears closer to the CE than the NIR emission peak. The reason for this is that the optical depth in the NIR is larger than in the MIR and the scattering efficiency is much higher in the NIR than in the MIR. Therefore, the NIR photons that we see might well be photons which were scattered by the dust along the axis of the torus, where the optical depth to the observer is smaller, and consequently, the center of light in the NIR may be displaced towards the scattering region, away from the CE. As the wavelength increases in the MIR, the optical depth to the CE decreases, as does the scattering efficiency, and the MIR light peak gets closer to the location of the CE.

Fig. 3 displays a superposition of the radio 6cm contours and the MIR 12.7$`\mu \mathrm{m}`$ contours over an image of the NLR in the \[OIII\]$`\lambda `$5007 line emission. The positions are set according to Sect. 2.4. The coincidence found between the NLR \[OIII\] emitting clouds and the MIR emission suggests that the same physical clouds are responsible for both emissions. With this registration, the radio source C is not coincident with the NLR-cloud B, as suggested in the preferred registration of Gallimore et al. (1996c). If one attempts to identify a particular NLR cloud responsible for the deflection of the radio jet, it would more likely be NLR-cloud C. The fact that NLR-cloud B does not affect the trajectory of the jet suggests that it lies on the foreground (or background) of the jet. The anti-correlation between the \[OIII\]+MIR images on one hand and the radio map on the other hand, is striking. It looks as if the \[OIII\] and MIR light were emitted principally in regions surrounding a cavity carved by the radio jet and localized at the edges of the ionizing cone. The matter inside the cone may have been swept away by the jet, and pushed to the actual location of NLR-cloud G, 1.7″ from S1 in the direction of the deflected radio jet. Gallimore et al. (1996c) originally proposed that the gaps in the NLR clouds may be caused by the jet, guessing an alignement between the MERLIN and the HST images, that turns out to be close to the actual alignement, as found in our analysis.

### 2.5 Map superpositions

Considering the absolute and relative positionings discussed above, Fig. 4 presents a series of superpositions of maps at different wavelengths and on different scale. The caption describes each map in detail.

Map (a) shows the excellent correlation between the large scale 6cm map and the 20$`\mu \mathrm{m}`$ image. At this resolution (0.6″), each 20$`\mu \mathrm{m}`$ peak has a counterpart on the 6cm map.

Map (b) reveals the good correlation existing between the low resolution (0.6″) 20$`\mu \mathrm{m}`$ image and the high resolution (0.1″) 12.3$`\mu \mathrm{m}`$ map. This demonstrates the very good consistency of the two datasets.

Map (c) shows the correlation between the \[OIII\] and the MIR emission and the anti-correlation between the \[OIII\] and the 6 cm radio emission on small scale. This is a strong evidence in favor of ionized clouds containing large amounts of cool dust.

Map (d) indicates that, as expected, the Br$`\gamma `$ emission is correlated with the \[OIII\] emission. It also shows that the Br$`\gamma `$ emission is strongly enhanced at the location of NLR-cloud E. Its particular excitation, the fact that it is not a strong MIR emitter (see Fig. 3) and its location with respect to the radio jet suggest that NLR-cloud E is undergoing a different excitation than the other NLR-clouds.

Map (e) highlights the very good correlation between the H2 2.12$`\mu \mathrm{m}`$ and the 12CO(2-1) maps. Still one notices the inverted brightness ratios in CO and H2 in the two western knots. On the contrary the eastern knot is very bright in both lines. This map also shows the spatial coincidence, at a 0.5″ scale, between the radio emission and the M-band emission. One even notices that the M-band emission bends together with the radio emission, at about 0.3″ North of the nucleus.

Map (f) summarizes the locations of the H2 2.12$`\mu \mathrm{m}`$ and Br$`\gamma `$ line emission, and the 12.3$`\mu \mathrm{m}`$ emission.

## 3 The unresolved torus emission

The observed IR SED of NGC 1068 and other AGN have been compared to predictions from radiative transfer models within dusty tori, in order to investigate the validity of the unified model and constrain possible geometries for the obscuring structure (e.g. Pier & Krolik, 1993; Granato & Danese, 1994; Efstathiou & Rowan-Robinson, 1995; Granato et al., 1997; Alonso-Herrero et al., 2001; Nenkova et al., 2002). A major problem is that it is still quite difficult to assess the truly nuclear SED to be compared with torus models, even for the best studied AGN, NGC 1068. In the NIR bands, even small aperture photometry is likely to include a non-negligible contribution from starlight. For this object, the best option so far is to combine the nuclear fluxes by Rouan et al. (1998) in the J-,H-,K-bands (corrected for star contamination) and by Marco & Alloin (2000) in the L-,M-bands with the nuclear MIR fluxes provided in Tomono et al. (2001) since they are all derived from high spatial resolution images. The resulting nuclear SED is shown on Fig. 5 and the flux values reported in Table. LABEL:n1068\_fluxes.

In the absence of precise physical ideas concerning the structure of the obscuring torus, several geometries – and associated free-parameters – are plausible and indeed have been investigated in the papers quoted above (flared discs, tapered discs, cylinders, with and without substantial clumping etc..).

An usually overlooked problem in this kind of study, is that dust optical properties around AGN are likely to be to some extent peculiar. As a matter of fact, Maiolino et al. (2001a, b) presented evidence for ’anomalous’ properties of dust in AGN. Here the term ’anomalous’ is relative to the standard properties of dust grains thought to be responsible for the average Milky Way extinction law and cirrus emission. As pointed out by Maiolino et al. (2001a, b), it is not surprising that the properties of dust grains might be very different in the dense and extreme environment of an AGN. In particular, they invoke a dust distribution biased in favor of large grains, ranging up to 1-10 $`\mu `$m, much larger than the usual cut at $`0.3\mu `$m of Mathis et al. (1977) type models (e.g. Draine & Lee, 1984; Silva et al., 1998). This solution has been introduced as the best, among the many explored, to explain the well established decrease of $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ and $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ ratios in the line-of-sight of AGN. It is worth noticing that the tendency for larger $`R_VA_V/(A_BA_V)`$ found in dense galactic regions are commonly explained in a similar way (Draine, 2001; Maiolino & Natta, 2002, and references therein).

Thus, we have performed SED fitting, adopting several different geometries, but also investigated models in which the size distribution of grains extends to radii, $`a_{max}`$, larger than the standard value. Of course, the uncertainty on optical properties of dust further limits the already moderate success of this approach in constraining the geometry (see below).

Due to the relatively long computing time required to produce a model, the SED fitting is performed through comparison with libraries of models. Each library consists of several hundreds of models belonging to a given ”geometry class”, in which typically 4-5 parameters are varied assigning to them 3-4 different values over a quite wide range. In particular we have compared the NGC 1068 nuclear SED with both anisotropic flared discs and tapered discs (for definitions, see Efstathiou & Rowan-Robinson, 1995). The former geometry consists of a structure whose height above the equatorial plane $`h`$ increases linearly with the distance in the equatorial plane $`r`$ (Fig. 1a in Efstathiou & Rowan-Robinson, 1995). Therefore the dust-free region (ignoring here dust possibly surviving in the NLR clouds) is exactly conical, with half-opening angle $`\mathrm{\Theta }_h`$. We also introduce (at variance with respect to Efstathiou & Rowan-Robinson, 1995) a dependence of the dust density $`\rho `$, and therefore of the radial optical depth, on the polar angle $`\mathrm{\Theta }`$. The general form used for this dependence and that on the spherical coordinate $`r`$ is:

$$\rho (r,\mathrm{\Theta })r^\beta \mathrm{exp}(\alpha \mathrm{cos}^2\mathrm{\Theta })$$

where $`\alpha `$ and $`\beta `$ are model parameters.

In tapered discs, $`h`$ increases with $`r`$ only up to a max value $`h_{\mathrm{max}}`$. In this case the dust density within the disc has been always kept independent of $`\mathrm{\Theta }`$. The two additional parameters used to fully characterize the models in both cases are the ratio between the outer and inner (i.e. sublimation) radii $`r_{max}/r_{min}`$, and the equatorial optical depth to the nucleus at $`9.7\mu `$m, $`\tau _{e,9.7}`$

The parameter values used in the libraries are:

* for both flared and tapered discs:

$$a_{max}=0.3,0.5,1.0\mu \text{m}$$

$$\beta =1,0,0.5,1,2$$

$$\tau _{e,9.7}=0.1,0.3,1,3,10$$

$$r_{max}/r_{min}=30,100,300$$
* for flared discs only:

$$\mathrm{\Theta }_h=20,30,40\text{°}$$

$$\alpha =0,1,3,6$$
* for tapered discs only:

$$\mathrm{\Theta }_h=30\text{°}$$

$$h_{max}/r_{max}=0.2,0.3,0.4$$

The impacts of changing various parameters on the predicted SED have been explored in Granato & Danese (1994) and Granato et al. (1997). Here we recall only that, roughly speaking, $`r_{max}/r_{min}`$ is related to the width of the IR bump, whilst $`\tau _{e,9.7}`$ controls mainly the NIR slope of the SEDs, as observed from obscured directions, as well as its anisotropy.

All possible combinations of parameter values are considered in the libraries. For each corresponding model we compute the SED for viewing angles from 0° to 90° in steps of 10 degrees. Then we search through these libraries which models provide a reasonable fit to the observed SED. The employed merit function is $`\chi ^2`$. We have assigned to each data point an error of 30%. This could well be an underestimate of the true uncertainty of the nuclear flux, especially in NIR bands, due to the difficulty of disentangling the starlight contribution. For instance, the nuclear K-band flux (corrected for stellar light) by Alonso-Herrero et al. (2001) in 3″ aperture is greater than our value by as much as a factor 3.

Indeed, the main result of our analysis is that several rather different combinations of geometry and parameters of the torus are perfectly compatible with the observed IR SED of NGC 1068 (examples of library fits are shown in Fig. 5 for the tapered models). In other words, as expected, the nuclear SED alone does not allow to put strong constraints on geometrical parameters. Several previous analysis have been too optimistic on this aspect, due essentially to an insufficient exploration of the parameter space. There is some tendency though, for the models with $`a_{max}`$ larger than the standard value to allow more combinations of the other parameters and still get an acceptable fit. We also confirm that some constraint can be put on the torus extension and on the optical thickness, and we find that good fits are obtained both with tapered (fits shown in Fig. 5) and flared disc models, provided the tori are moderately thick, in the sense discussed by Granato et al. (1997). The development of interferometric observations should, in the future, provide more observational constraints.

## 4 The extended MIR emission associated with the ionizing cone.

MIR high resolution images of NGC 1068 show elongation perpendicularly to the putative dusty torus-like structure (Alloin et al., 2000; Bock et al., 2000; Tomono et al., 2001). This is interesting in two ways: (i) the extension of the emission, of the order of 70 pc, exceeds the estimated size of the torus, which is not expected to be larger than 20 pc (Granato & Danese, 1994; Efstathiou & Rowan-Robinson, 1995; Granato et al., 1997), and (ii) it traces approximately the two ionization cones where dust grains are less likely to survive. Also, this emission is more prominent on the North-East side of the CE (which is believed to trace the cone closer to our line of sight) than on the opposite South-West side. The decline of the surface brightness with distance from the CE is relatively shallow and the observed ”local” SEDs of the emission (see Fig. 6) show interesting trends: in the North-East side, the SEDs along the axis have a steep decline towards shorter wavelengths, while in the South-West side, the SEDs appear to be flatter.

### 4.1 Diffuse dust in the cone?

Is it possible to explain the observations with diffuse dust distributed in the cone? To investigate this hypothesis, we have run several sets of models with diffuse dust in the cone extending up to $`1`$″, with different density and grain size distributions, and could not find a good match with the observations. Such models are affected by either, but usually both, of the following problems: (i) the emitted specific intensity decreases more steeply than the observed MIR maps show, and/or (ii) the local SEDs in the cones are too flat. The dust volume emissivity can be enhanced either by increasing the dust density or its temperature. The first possibility works only up to a certain point above which the dust in the cone becomes optically thick to optical-UV photons and, as a consequence, is less heated by the CE. We have tested with our model that this “turnover” is reached before the observed emission can be produced (see also the analytical estimate below). Moreover, an excess density of diffuse dust in the cone would obscure the CE and the BLR from any direction, and would also require a different mechanism than direct photo-ionization from the CE to produce the conical NLR.

Let us perform as well an analytical calculation in order to evaluate when this “turnover” is reached. Let us assume a distribution of homogeneous diffuse dust in the cone, optically thin in IR (not necessarily in optical-UV). The observed specific intensity $`I_{\nu ,18}`$ at 18 $`\mu `$m is:

$$I_{\nu ,18}=l\alpha _{\nu ,18}B_{\nu ,18}(T)=\tau _{\nu ,18}(l)B_{\nu ,18}(T)$$

$`\tau _{\nu ,18}(l)`$ is the optical depth along the line of sight crossing the cone (l), $`\alpha _{\nu ,18}`$ the extinction coefficient and $`B_\nu `$ is the Planck function at the frequency $`\nu `$. Specific quantities are all evaluated at 18$`\mu \mathrm{m}`$. Given the geometry, we can evaluate the dust temperature $`T`$ assuming that the heating is only from the CE and that the distance from the central source, $`r`$, is almost constant along the line-of-sight. Also, $`rl`$. Thus we can identify the optical depth $`\tau _{\nu ,18}(l)`$ with the optical depth to the central source $`\tau _{\nu ,18}(r)`$. For graphite grains of typical size, the equilibrium dust temperature $`T`$ can be estimated as (Barvainis, 1987):

$$T=1650\left(\frac{L_{UV,46}}{r_{pc}^2}\right)^{1/5.6}\mathrm{exp}\frac{\tau _{UV}(r)}{5.6}$$

Setting $`r=70`$ pc (i.e. 1″, where the observed $`I_{\nu ,18}`$ is $`10`$ Jy arcsec-2) and $`L_{UV}=0.02\times 10^{46}`$ $`\mathrm{erg}\mathrm{s}^1`$ we have

$$T=180\mathrm{exp}\frac{\tau _{UV}(r)}{5.6}$$

On the other hand, for standard Galactic dust mixture, $`\tau _{UV}80\times \tau _{\nu ,18}`$. This leads to:

$$I_{\nu ,18}=\frac{\tau _{UV}(r)}{80}B_{\nu ,18}\left(180\mathrm{exp}\frac{\tau _{UV}(r)}{5.6}\right)$$

This function of $`\tau _{UV}`$ has a maximum of $`10`$ Jy arcsec-2 (i.e. more or less the observed level) at $`\tau _{UV}1`$, but as long as the dust is optically thin in the cones at optical-UV wavelengths – the condition for photoionization – the predicted 18$`\mu \mathrm{m}`$ brightness is significantly lower than observed (for instance by a factor 4 for $`\tau _{UV}0.1`$, or by a factor 100 for $`\tau _{UV}0.01`$).

To increase the dust temperature, it is necessary to invoke some additional heating source for the dust in the cone, such as anisotropic emission from the CE, young stars, or even shock heating (see also discussion in Villar-Martín et al., 2001, and Sect. 4.3). While all these are rather likely to occur, they are not satisfactory explanations in this respect: indeed the dust temperature in the cone, under the assumption of optically thin dust heated only by an isotropic central source, is already too high ($``$ 200-300 K at $`0.5`$″, depending on the dimension and composition of grains), to produce a SED declining as steeply as observed towards short wavelengths. A higher T, as could result from an additional heating source, would worsen the problem, independently of the nature of this source.

Therefore, we are led to conclude that a distribution of diffuse dust within the ionizing cone cannot account for the MIR observed SEDs and flux levels.

### 4.2 Optically thick dust clouds in the cone

The ”local” SEDs of the emission displayed in Fig. 6 provide another key-element for understanding the MIR emission: the SEDs in the North-East side have a steep decline towards shorter wavelengths and show silicate 9.7 $`\mu `$m feature in absorption. Both the SED steepness and the depth of the silicate feature increase with distance from the CE. On the contrary, in the South-West side, where the SEDs appear to be flatter, the 9.7$`\mu \mathrm{m}`$ feature is seen in emission. In other words, observations suggest MIR optically thick dust emission in the North-East cone and optically thin emission in the South-West cone. The most natural way to explain this behavior is through the presence of dust within the conical regions, organized in individual optically thick (in the MIR) clouds, and having a small covering factor. Indeed, in such a configuration each cloud would see directly the CE, and would have a much higher dust temperature on its side facing the nucleus than on its opposite side. As a consequence, the overall emission from each individual cloud would be highly anisotropic. In the North-East cone (approaching us) we observe the clouds mostly from their colder surface, while the contrary holds true in the South-West cone (receding from us). Therefore, we expect to see exactly the trend depicted in Fig. 7: the North-East side shows SEDs similar to that labeled $`\alpha =0\mathrm{°}`$, while the South-West SEDs are similar to those of the kind labeled $`\alpha =180\mathrm{°}`$.

A picture where the dust is organized in MIR optically thick clouds gives a natural explanation to both the observed emission levels, as well as to the SED shapes, including their slope and silicate feature. The surface brightness of the extended MIR emission can be explained by this model, without invoking additional heating sources. Yet, the observations do not discard either the presence of other sources of heating, which we discuss in the coming section.

### 4.3 Possible additional heating mechanisms

In the former section, we discarded a picture where diffuse dust inside the ionizing cone, heated by the CE, was responsible for the observed MIR emission. On the contrary, the local SEDs between 7 and 20$`\mu \mathrm{m}`$ strongly favor a model in which optically thick dusty clouds, with small covering factor, produce the MIR emission. In this context, the radiation from the CE is the major source of heating for these clouds, but the possibility and evidence of additional heating mechanisms ought to be discussed here.

First, the radiation from the CE might be beamed. This would add a continuum source in the direction of the beam enhancing the dust emission along this region.

A second possibility would be additional heating of the dust by a starburst. But in the case of NGC 1068, the observed starburst ring is too far away (10″-20″) from the MIR structures discussed here. An alternative possibility could be the presence of local star formation along the jet (jet induced star formation). Our superposition of the small scale MIR and radio maps (Fig 4) shows however that the radio emission is much more spatially focused than the MIR emission, and consequently, the hypothesis of jet-induced star formation is not well supported by observations.

Finally, interactions between radio and optical structures generate shocks that can heat the gas to very high temperatures ($`10^6`$ K). This gas then becomes an efficient source of UV and X-ray photons able to ionize the surrounding material. Among others, an expected effect is extra heating of the dust. There are indeed some arguments in favor of a contribution of heating by shocks.

– the subarcsec NIR continuum emission (Rouan et al., 1998; Marco & Alloin, 2000) is elongated along an axis reminiscent of the radio emission axis, and particularly, the 3.5$`\mu \mathrm{m}`$ emission shows a change of direction which follows that of the radio jet (see Fig. 4.e).

– the subarcsec resolution Chandra (0.25-7.50 keV) image published by Young et al. (2001) reveals a bright nucleus and bright extended emission up to 5″ to the North-East, probably coincident with the extended North-East radio lobe observed at 6cm. Moreover, this image shows that the X-ray emission immediately to the South of the nucleus has a more or less North-South extension, just like the radio emission. Could this X-ray emission be contributed by the shocked gas?

Axon et al. (1998) find signs of interaction between the radio and optical structures. For example, they find that emission lines are split into two velocity components separated by $``$1500 $`\mathrm{km}\mathrm{s}^1`$ at the position of the jet axis. Outside this region, the lines are quiescent. Supporting the shock model is the finding by these same authors of an excess of local optical continuum associated with the radio jet. The existence of an excess of UV continuum is supported by the increase in excitation along the jet axis (outside the jet, the gas excitation is much lower). This continuum excess should naturally contribute to the dust heating, as shown by Villar-Martín et al. (2001). Clues that shocks are present also come from the H2 2.12$`\mu \mathrm{m}`$ and Br$`\gamma `$ high resolution spectroscopy of Galliano & Alloin (2002). At the position of the North-Br$`\gamma `$-H2 knot, corresponding to the position of NLR-cloud E, the H2 emission is double peaked ($``$ 400$`\mathrm{km}\mathrm{s}^1`$ velocity difference), showing a kinematical perturbation at this location, possibly induced by shocks. Br$`\gamma `$ is also very broad at this location, displaying a FWHM of $``$1000 $`\mathrm{km}\mathrm{s}^1`$, against $``$200 $`\mathrm{km}\mathrm{s}^1`$ for H2. Can we explain this difference with shocks? In the absence of entrainment processes, the cores of the clouds impacted by the jet will be accelerated to a velocity $`V`$ given by:

$$VV_s\sqrt{\frac{n_h}{n_w}}$$

where $`V_s`$ is the velocity propagation of the shock in the hot phase (that can be identified to the jet advance speed) , and $`n_h`$ and $`n_w`$ are the densities of the clouds (ionized or molecular) and of the hot phase ($`10^7`$ K). This equation comes from applying the balance between the pressure in the shocked inter-cloud medium ($`P_hn_h\times V_s^2`$), and the pressure behind the shock in the cloud ($`P_wn_w\times V^2`$, Klein et al., 1994).

Molecular clouds will be accelerated to lower velocities than ionized clouds, due to the higher density (the shock finds higher resistance). Comparing the two velocities:

$$\frac{V_{B_\gamma }}{V_{H_2}}=\sqrt{\frac{n_{B_\gamma }}{n_{H_2}}}$$

Spectroscopy of the ionized regions (i.e. Br$`\gamma `$ emitting gas, but not shocked) not interacting with the radio jet implies densities of a few 100 cm-3 (Axon et al., 1998). With the molecular gas being $``$ 100 times denser, then $`V_{Br_\gamma }/V_{H_2}\sqrt{100}=10`$. Therefore, the ionized clouds can easily be accelerated to speeds $``$10 times higher than the molecular clouds. This would explain the observed difference in the spectral profiles of H2 and Br$`\gamma `$ in the North-Br$`\gamma `$-H2 knot. An alternative possibility is that only molecular clouds impacted by slow shocks can survive, while those enduring high velocity shocks are destroyed. In this case, the H2 emission would come from more quiescent clouds.

Therefore, it is likely that shocks, witnessed by the kinematics of the molecular material, play a role in the heating of the dust, although this remains to be quantified.

### 4.4 Conclusion

The contribution of shocks to the heating of dust and molecular material looks plausible and promising to explain detailed kinematical features. However, we find that the main source of extended MIR emission to the North-East and South-West are MIR optically thick dust clouds distributed within the ionizing cone, with small covering factor.

## 5 The molecular emission perpendicular to the cone axis

The material in the plane perpendicular to the ionization cone axis is visible essentially through molecular line emission of the warm/cold material. No conspicuous NIR or MIR continuum emission is detected along the East-West direction. In particular, the strong CO(2-1) and H2 2.12$`\mu \mathrm{m}`$ emission lines witness the presence of molecular material up to 100 pc along the East-West direction (Schinnerer et al., 2000; Galliano & Alloin, 2002). The structures revealed through these lines correlate remarkably well.

The 0.5″ resolution H2 2.12$`\mu \mathrm{m}`$ emission line map of Galliano & Alloin (2002) shows that the H2 emission is distributed in two unequally bright knots, located at 1″ (70 pc) from the CE, the brightest to the East and the weakest to the West. Each of the two knots displays a North-South extension which even extends further like an arc or spiral. The East-H2 knot is 3 times more intense than the West-H2 one. The CO emission line map has a similar aspect but presents a smaller East-knot to West-knot intensity ratio (1.5 for CO against 3 for H2). The fact that the CO/H2 ratios are different in the East knot and the West knot suggests that extinction plays a role.

There is another reason to believe that differential extinction is likely to occur in NGC 1068, from the observational facts reported on hereafter. Several studies (Schinnerer et al., 2000; Baker, 2000) argued, on the basis of the CO kinematics, about the presence of a warp in the molecular disc surrounding the AGN of NGC 1068. Moreover, the maps at different scales (see Fig. 4) show that the PA of the main axis of the observed North-South structures change with scale, which supports the idea of a warp. X-ray observations show that the material towards the central engine is Compton thick, which implies a column density $`10^{24}\mathrm{cm}^2`$. This Compton thick material (which we will call the X-ray absorber hereafter) is likely to be located at the very central parts of the AGN, and might be e.g. a corona of hot gas such as the atmosphere of the accretion disk (Collin-Souffrin et al., 1996). Since the excitation of molecular emission is dominated by the X-rays from the CE (Maloney, 1997), a tilt between the orientations of the X-ray absorber and the molecular disk, which would naturally occur in the presence of a warp, would break the axisymmetry and hence could significantly affect the observed molecular line intensity distribution.

### 5.1 Differential extinction between the Eastern and Western knots

If one assumes that both the East and West knots have intrinsically the same intensity, and that the western emission is dimmed by dust extinction only, a column density towards the West-H2 knot of the order of $`N_{West}=\mathrm{0.5\hspace{0.17em}10}^{22}\mathrm{cm}^2`$ (for $`A_K/A_V=0.11`$ and $`A_V/N_H=\mathrm{0.19\hspace{0.17em}10}^{22}\mathrm{mag}(\mathrm{nuclei}\mathrm{cm}^2)^1`$) is required. This is of the order of the column density through typical molecular clouds in our Galaxy. Dust extinction affects the H2 line photons, but do not affect the CO line photons at 1.3 mm. Still, both CO and H2 photons can be absorbed by molecules in intervening clouds. Large and dense molecular clouds can in fact screen completely the H2 and CO photons. The estimation of the optical depth through an $`N_H=10^{22}\text{cm}^2`$ cloud with $`n_{CO}/n_{\text{H}\text{2}}=10^4`$ results in $`\tau >10`$ for photons with $`\nu (2\sigma _{turb})\nu \nu (2\sigma _{turb})`$, where $`\nu `$ is the frequency of the crossing photon, $`\sigma _{turb}`$ is the velocity dispersion of molecules in the cloud, and $`\nu (2\sigma _{turb})`$ is the frequency of the molecular transition for a molecule at radial velocity of $`2\sigma _{turb}`$.

In the following, we show quantitatively that a small tilt between the molecular disc and the X-ray absorber can lead to a notable intensity difference in the H2 and even, but to a much lesser extent, in the CO emission, between the East and West knots.

### 5.2 The model of X-ray excitation

The excitation of the molecular emission is computed using the results of the model of X-ray irradiated molecular gas developed by Maloney et al. (1996). This model computes the physical and chemical state of dense neutral gas exposed to an intense X-ray flux. We use the numerical predictions of Maloney et al. (1996) for a power-law X-ray source of index $`\alpha =0.7`$ ($`F_\nu \nu ^\alpha `$) and consider molecular clouds of density $`n=10^5\mathrm{cm}^3`$ through which the column density is $`N=10^{22}\mathrm{cm}^2`$. The most important parameter that controls the physical conditions of the molecular material is the ratio of the attenuated X-ray flux to density. This parameter is expressed in terms of an effective ionization parameter:

$$\xi _{eff}100\frac{L_{44}}{n_5r_{pc}^2N_{22}^{0.9}}$$

where the 1-100 keV X-ray luminosity of the central source is $`L_X=10^{44}L_{44}`$$`\mathrm{erg}\mathrm{s}^1`$, the distance to the X-ray source is $`r_{pc}\mathrm{pc}`$, the gas density is $`n=n_5\mathrm{\hspace{0.17em}10}^5\mathrm{cm}^3`$ and the attenuating column density to the X-ray source is $`N_{att}=10^{22}N_{22}\mathrm{cm}^2`$. The predictions of Maloney et al. (1996) for the emergent surface brightness of a cloud of density $`10^5\mathrm{cm}^3`$ in the 2.12$`\mu \mathrm{m}`$ line and CO rotational lines (all lines together) are reproduced on Fig. 8.

In the case of NGC 1068, $`L_{44}=1`$ (Maloney, 1997; Matt et al., 1997). Since the maxima of H2 emission occur at 1″(70 pc) from the CE, we know that the value $`\xi _{eff}=0.026`$ (corresponding to the value of $`\xi _{eff}`$ at the maximum of $`I/N_{22}`$ for H2 in Fig. 8) corresponds to $`r_{pc}`$=70. This gives us access to the column density between the CE and the H2 knots, using the expression of $`\xi _{eff}`$:

$$N_{22}=\left(\frac{100L_{44}}{\xi _{eff}n_5d_{pc}^2}\right)^{\frac{1}{0.9}}=0.75$$

Therefore, this value is much smaller than the column density of $`10^{24}\mathrm{cm}^2`$ to the CE, deduced from X-ray observations. This supports the hypothesis that the central X-ray absorber and the emitting molecular material are not coplanar. Hence, we have developed a new model along this line of arguments, which is presented in the next section.

### 5.3 Molecular disc model

The configuration sketched in Fig. 9 consists of a molecular disc and of a central axisymmetrical X-ray absorber, tilted with respect to the molecular disc. In Fig. 9, the grey rectangle represents a cut through the edge-on view of the molecular disc. The central X-ray absorber axis is tilted with respect to the axis of the disc. The region filled with the “hexagon” patterns are the regions where the X-rays are deeply attenuated, and consequently from which no strong emission can emerge. The circle centered on the CE traces the radius at which, in the absence of a central X-ray absorber, the value of $`\xi _{eff}`$ would correspond to the maximum of emission. The regions from which the emission will actually arise are shown with the hatched circles. The observer sees directly the “non obscured emitting region”, while it sees the “obscured emitting region” through the molecular disc. This sketch provides a qualitative explanation to the fact that the molecular emission is offset from the CE, and how it can naturally appear brighter on one side.

We have built a simple numerical model to simulate such a configuration. It consists of the following two components:

(1) A disc distribution of molecular clouds around the central X-ray source. All clouds are identical to the densest cloud studied by Maloney et al. (1996), i.e. their density is $`10^5\mathrm{cm}^3`$ and their column density is $`N=10^{22}\mathrm{cm}^2`$. Complex cloud number density distributions describing flared or warped discs were computed as well, but a simple distribution (with less free parameters) is sufficient for getting interesting results. In an (x,y,z) reference frame, where z is the axis of the disc, we have chosen the following distribution:

$$d_c(z)=D_ce^{(|z|/Z)^\lambda }$$

where $`d_c`$ is the cloud number density, $`D_c`$ is the density at z=0, and $`(Z,\lambda `$) controls the dependence of the density on the height above the equator of the disc. A flat rotation curve was assigned to the disc (plateau at 130 $`\mathrm{km}\mathrm{s}^1`$), in order to estimate the opacity of the clouds in the lines.

(2) A central axisymmetrical X-ray absorber, for which the equatorial attenuation is very high ($`10^{24}\mathrm{cm}^2`$). We have set the azimuthal dependence of its attenuating column density to the following form:

$$N_{att}(\mathrm{\Theta })=N_{att}(0)\mathrm{\hspace{0.17em}10}^{sin^\beta (\mathrm{\Theta })},$$

where $`\mathrm{\Theta }`$ is the azimuthal angle and $`N_{att}(0)`$ the equatorial thickness of the X-ray absorber, that is $`N_{att}\mathrm{\hspace{0.17em}10}^{24}\mathrm{cm}^2`$.

These two components can be freely and independently rotated in space.

For each cell of a 3D grid we have computed the value of the attenuation towards the X-ray source, as the sum of the attenuations due to both components. The value of $`\xi _{eff}`$ in each cell is then computed.

The attenuation of the emission line photons on their way to the observer is computed taking into account two effects. (1) the dust extinction for H2 only, (2) the high opacity of some intervening molecular clouds in both H2 and CO lines (see Appendix A).

### 5.4 Results

Fig. 10 and Fig. 11 display the results of an illustrating model for a seeing of 0.4″.

In Fig. 10, the top two figures (a) and (b) present an edge-on view of the disc. For a plane perpendicular to the disc, passing through the center, image (a) shows the cloud density distribution, and image (b) the values of $`\xi _{eff}`$ for an X-ray absorber tilted with respect to the molecular disc by 15° in the direction of the x axis. The subsequent simulated images (c,d,f) are not seen edge-on anymore. Images (c) and (d) respectively show the actual orientations of the chosen model for the molecular disc and the X-ray absorber.

Fig. 11 compares model predictions with observations: (a) displays the H2 emission map and (b) the model prediction for the H2 emission. (c) displays the CO(2-1) emission map and (d) the model prediction for the CO (all lines) emission.

The parameter values corresponding to this model are: $`D_c=100\mathrm{pc}^3`$, $`Z=20\mathrm{pc}`$, $`\lambda =1.15`$, $`N_{att}(0)=10^{26}\mathrm{cm}^2`$, $`\beta =1.5`$, $`\sigma _{cin}=100`$$`\mathrm{km}\mathrm{s}^1`$ and $`\sigma _{turb}=20`$$`\mathrm{km}\mathrm{s}^1`$

The orientation of the two components is illustrated by the ellipses of Fig. 10 (c) and (d).

This two component model reproduces quite well the most characteristic features of the observations of the following H2 2.12$`\mu \mathrm{m}`$ and CO emission: The distance of the observed emission knots to the CE, the extended arc or spiral structure, the intensity ratios between the East and West knots in H2 and CO (3.2 for H2 and 1.25 for CO) and is consistent with the Compton thickness of the absorber along the line-of-sight. Yet, the surface brightnesses of the simulated emission knots are too large compared to the observations. Choosing a smaller cloud number density lowers the intensity, but lowers as well the intensity difference between the East and West knots. This suggests that there is probably a difference in the amount of material between the East and the West knots. If we consider such a difference, we can then fit more easily the surface brightnesses and the intensity (all lines together) ratios at the same time. Still, in that case a configuration with selective extinction is required.

Another point is that, for the chosen model, the predicted CO emission knots are more extended than observed. This can be resolved by adjusting the dependance of the cloud density distribution on the distance to the central engine. The fact that we predicted the total CO rotational emission intensity explains part of this discrepancy. We also point out that, because CO is still emitted at lower $`\xi _{eff}`$ values than H2, in the configuration we present, the H2 and the CO peaks do not fall exactly at the same position, the CO emission being slightly rotated counter-clockwise with respect to the H2 emission.

An interesting observational feature, which is not reproduced by this model, is the increase of CO emission towards the South in the western knot, a feature which is not observed for the H2 emission. In the kinematical configuration proposed by Galliano & Alloin (2002) to fit the H2 line profiles, both a rotating and an outflowing kinematical components had to be introduced. One notices that it is precisely in the South-West and North-East regions that the projected velocities of both components have the greatest difference (where the profiles are the most asymmetric). In the western region, the outflowing component would be more absorbed by the rotating disc to the North (the projected velocities tend to be more equal), than to the South. Hence, this would qualitatively explain a rise of CO emission to the South of the western knot.

## 6 Summary and conclusion

We have discussed in this paper several features of the central 100 pc region around the CE in the AGN of NGC 1068.

The second section was dedicated to revisiting the relative positioning of maps at different wavelengths. We provide a new estimate of the absolute positioning of the unresolved peak in the K-band. After compiling the most important observations published, and after elucidating the apparent disagreement between different results, we could register with a high precision the radio, NIR, optical and UV continuum maps, as well as the submillimiter, NIR, optical and UV emission line maps. We produced a series of map overlays, which constitute an essential tool for giving a general view of the complexity of the central region of NGC 1068. These superpositions show that the MIR is very well correlated with the \[OIII\] emitting NLR clouds, revealing the survival of dust grains in the UV-irradiated ionization cone. The final registration of the small scale high resolution radio maps with respect to the HST optical maps show a surprising anti-correlation: no NLR cloud is found along the radio jet trajectory. This suggests a picture where the radio jet sweeps away material within and along the axis of the ionizing cone.

The third section deals with the modeling of the NIR-MIR unresolved peak emission in NGC 1068. First, we could show that the SED fitting approach is somewhat limited since several acceptable combinations of torus geometries and dust parameters are compatible with the observed IR SED of NGC 1068. We conclude that no strong constraint on the geometry of the torus can be set using this approach, while its extension and thickness can to some extent be estimated: the torus in NGC 1068 is smaller than 10-20 pc and moderately thick in the MIR.

Then, a discussion is provided about the origin of the observed extended MIR emission, mainly along the axis of the ionizing cone. We find that the most natural way to match the observations is to consider individual optically thick clouds, distributed in the NLR and having a small covering factor. An alternative picture, consisting in diffuse dust distributed within the cone, is confidently discarded. Possible additional sources of heating are discussed, and we show in particular that the likely presence of shocks can explain some of the observed kinematical features such as the difference between the FWHM of the molecular emission lines and of the ionized gas emission lines, in the North-Br$`\gamma `$-H2 knot in particular.

The last section presents a model of the excitation and distribution of the H2 and CO line emission along the East-West direction (perpendicular to the ionizing cone axis). We argue that the observed intensity differences between the eastern and the western regions, with respect to the CE, can be partly due to selective extinction produced by a slight tilt between the orientations of the Compton thick X-ray absorber in the very central region of the AGN and the molecular material distribution at larger scales.

###### Acknowledgements.

We would like to thank Daigo Tomono, Linda Tacconi and Eva Schinnerer, who kindly provided their data in the MIR and in the CO emission line respectively. We also aknowledge interesting discussions with Laura Silva, about the extended MIR emission in NGC 1068. Finally, we are gratefully indebted to the referee J. Gallimore for valuable comments. E.G. aknowledges the support of an ESO studentship and of the french ministery of Foreign Affairs, as a cooperant.

## Appendix A Attenuation of H2 and CO photons by intervening molecular clouds

We computed that a molecular cloud of column density $`N=10^{22}\text{cm}^2`$ is completely opaque ($`\tau 10`$) to CO and H2 photons, if the radial velocity difference between the emitting molecule and the molecular cloud is $`2\sigma _{turb}`$, where $`\sigma _{turb}`$ is the velocity dispersion inside the cloud (of the order of 10-20$`\mathrm{km}\mathrm{s}^1`$). Let us estimate the corresponding attenuation effect.

We consider two cells of the grid aligned along the line of sight, $`cell_a`$ and $`cell_b`$. $`Cell_a`$ contains one emitting cloud called $`cloud_a`$ and $`cell_b`$ contains one non-emitting cloud called $`cloud_b`$. $`Cell_b`$ is located between $`cell_a`$ and the observer. We can estimate a cross section $`S_{eff}`$ of radius $`r_{eff}`$ for $`cloud_a`$, so that we can consider that, when the center of $`cloud_b`$ lies inside $`S_{eff}`$, $`cloud_a`$ is completely covered, and when the center of $`cloud_b`$ lies outside, then $`cloud_a`$ is completely uncovered. The quantity $`r_{eff}`$ is such that the overestimation in covering when the center of $`cloud_b`$ does not coincide with the center of $`cloud_a`$ (but still lies within $`S_{eff}`$) is balanced by the underestimation in covering when the center of $`cloud_b`$ is outside $`S_{eff}`$ (but still covers part of $`cloud_a`$). We find $`r_{eff}`$=0.89$`r_{cloud}`$.

Then we can simply compute the absorption in terms of probability for the center of $`cloud_b`$ to lie within the cross section of $`cloud_a`$ in the area of a cell $`S_{cell}`$. This prbability is $`P_1=S_{eff}/S_{cell}`$. In P1$`\times `$100% of the cases, $`cloud_b`$ can be considered to be hiding completely $`cloud_a`$, and in (1-P1)$`\times `$100% of the cases, $`cloud_b`$ does not hide $`cloud_a`$ at all.

We can define as well $`V_{eff}`$ (the analogous to $`r_{eff}`$ in the velocity space) such that, if $`V_{eff}V_a`$-$`V_bV_{eff}`$, we can consider that the photon, in average will be absorbed, where $`V_a`$ and $`V_b`$ are the respective velocities of $`cloud_a`$ and $`cloud_b`$ along the line of sight. We find $`V_{eff}=1.2\sigma _{turb}`$.

Let $`\sigma _{cin}`$ be the velocity dispersion between the cloud velocities. If $`\mathrm{\Delta }V=abs(V_aV_b)`$ is the difference between the radial velocities of $`cloud_a`$ and $`cloud_b`$, then the probability for $`cloud_b`$ to have the same velocity as $`cloud_a`$ is: $`P_2=V_{eff}(\sigma _{cin}\mathrm{\Delta }V)/\sigma _{cin}^2`$.

Finally, in $`(P1\times P2)`$% of the cases, the photon emitted by $`cloud_a`$ will be absorbed by $`cloud_b`$. Thus, we approximate the attenuation by molecules in intervening clouds, by a factor $`(1P_1\times P_2)`$ per cloud in the cells located along the line of sight to the observer.