Extracting the Quantum Hall Conductance from a Single Bulk Wavefunction

Ruihua Fan
Rahul Sahay
Ashvin Vishwanath
Department of Physics, Harvard University, Cambridge, MA 02138, USA

We propose a new formula that extracts the quantum Hall conductance from a single (2+1)D gapped wavefunction. The formula applies to general many-body systems that conserve particle number, and is based on the concept of modular flow: i.e., unitary dynamics generated from the entanglement structure of the wavefunction.
The formula is shown to satisfy all formal properties of the Hall conductance: it is odd under time-reversal and reflection, even under charge conjugation, universal and topologically rigid in the thermodynamic limit.
Further evidence for relating the formula to the Hall conductance is obtained from conformal field theory arguments.
Finally, we numerically check the formula by applying it to a non-interacting Chern band where excellent agreement is obtained.

Since the discovery of the quantum Hall effect, the investigation of gapped topological phases has been at the forefront of quantum many-body physics <cit.>.
Such phases cannot be characterized by any local order parameters and all of their unusual properties emerge from the entanglement of the ground state wavefunction <cit.>.
Quantifying this connection has been a persistent and interesting challenge for the past few decades.

For example, entanglement entropy in a ground state wavefunction reveals the total quantum dimension of anyons  <cit.> while access to multiple ground states on a torus allows for a determination of modular matrices characterizing the anyons <cit.>. Further information, such as the entanglement spectrum, reveals the protected edge states
Yet, how to extract the quantum Hall conductance, the first discovered topological quantity, from the entanglement pattern of a gapped wavefunction alone remains an open question.
Here, by gapped wavefunction, we mean the ground state wavefunction of a gapped, local Hamiltonian. Besides its foundational interest, such a formula could help unravel the entanglement structure of local quantum systems and thus characterize limits to the efficient representations of quantum states by tensor networks or stabilizers <cit.>.

There are several approaches to extracting Hall conductance that do not invoke entanglement.
For free fermion systems, there are the TKNN formula <cit.> and Fredholm index formula <cit.>.
While the latter idea has been generalized to the interacting case, its formalism uses quasi-adiabatic evolution and is more convenient to work with if the actual Hamiltonian, rather than a single wavefunction, is known <cit.>.
We also note earlier works on extracting the many-body Chern number from a single wavefunction based on connections to topological quantum field theory and surgery <cit.>. This quantity coincides with the  Hall conductance in free-fermion systems but is a distinct quantity in generic interacting systems (see <cit.> for details).

In this letter, we conjecture a new formula, with analytical and numerical support, which uses entanglement to compute the Hall conductance from a single wavefunction directly. The formula is proposed for generic interacting systems without assuming translational invariance [The close relation between entanglement spectrum (ES) and protected edge states suggests a possible strategy based on the ES. However, typically, translation symmetry needs to be invoked to fully extract information present in the ES.].

The Hall conductance relies on charge conservation, which the entanglement spectrum of a state is blind to [One can in principle use the $U(1)$ resolved ES and study the spectral flow, which, however, requires a continuous family of states and is not convenient to work with.].

Therefore, more data is necessary, motivating our use of the modular Hamiltonian and modular flow in the formula presented here.

In a recent advance, the modular Hamiltonian has appeared in a proposed formula for the chiral central charge <cit.>, which can be viewed as entanglement transport under the modular flow <cit.>.
(See also <cit.> for related works.)

Here we are interested in the Hall conductance $\sigma_{xy}$, a distinct topological invariant, which requires an additional conserved $U(1)$ charge $Q$.

Proposed Formula for Hall Conductance—
Consider a lattice system on the plane with a $U(1)$ symmetry that is generated by the total charge operator $Q = \sum_i n_i$, where $i$ labels the lattice sites.
Let $\ket{\psi}$ be a $U(1)$ symmetric gapped wavefunction, $\rho_D = \Tr_{\bar{D}} \ket{\psi} \bra{\psi}$ the reduced density matrix of any given region $D$.
The modular Hamiltonian is defined as $K_D := -\ln \rho_D$, which is a Hermitian operator with a lower-bounded spectrum and conserves the total number of charges $Q_D$ in that region, i.e., $[Q_D,K_D]=0$.
The modular Hamiltonian generates a unitary evolution $U_D(s) := e^{-\i sK_D}$, called modular flow, where the modular time $s$ is dimensionless <cit.>.
The Hall conductance is encoded in a charge response under modular flow at the linear order in $s$.

Specifically, divide the plane into four regions $A,B,C$ and $D$, whose linear sizes are larger than the correlation length, and every three regions meet at a point once and only once <cit.>.
Let $K_{\bullet}$ and $Q_{\bullet}$ denote the modular Hamiltonian and charge of the corresponding region.
The response function we focus on is:

\begin{equation}
\label{eq:sigmaxy formula 1}
\begin{gathered}
\includegraphics[width=2cm]{ABCdemo.pdf}
\end{gathered}\quad
\begin{aligned}
\Sigma(\psi;\, & A,B,C) \\
&= \frac{\i}{2} \braket{\psi|[K_{AB},Q_{BC}^2]|\psi}\,,
\end{aligned}
\end{equation}

i.e., the charge response in $BC$ under the modular flow on $AB$.
The modular Hamiltonian itself is non-universal and ultraviolet divergent, but the expectation value of the commutator will be argued to be topological and insensitive to details of the wavefunction.
Our proposal for Hall conductance $\sigma_{xy}$ is:

\begin{equation}
\label{eq:conjecture}
\Sigma(\psi;A,B,C) = \sigma_{xy}\,.
\end{equation}

Note that the left-hand side containing $Q_{BC}^2$ is consistent with the fact that the unit of Hall conductance is $e^2/\hbar$.

The rest of the letter presents evidence for this formula from various angles.
We first show that the formula shares the same general properties as the Hall conductance.
We then establish their quantitative connection via the bulk-edge correspondence and conformal field theory (CFT) arguments.
Our discussion, which uses the defect operator, provides a unifying viewpoint on the two different formulas for the chiral central charge and Hall conductance.
Finally, we provide numerical evidence using both lattice and continuum free-fermion models.

General properties—
Any formula for the Hall conductance needs to meet the following minimal requirements.
First of all, $\Sigma(\psi;A,B,C)$ should be real and additive under stacking, i.e., on combining independent systems:

\begin{equation}
\begin{aligned}
\Sigma(\psi_1\otimes\psi_2;&A,B,C) = \\
&\Sigma(\psi_1;A,B,C) + \Sigma(\psi_2;A,B,C)\,,
\end{aligned}
\end{equation}

which are clear from its definition.
Although additivity may seem trivial, other topological quantities e.g., the many-body Chern number do not possess this property.
We now show that $\Sigma(\psi;A,B,C)$ also satisfies other requirements that are not manifest in the construction.
It has the same $CRT$ transformation rules as the Hall conductance, namely, it (i) is even under charge conjugation:

\begin{equation}
\label{eq:charge conjugation}
\calC: \quad \Sigma(\psi;A,B,C) = \Sigma(\calC\psi;A,B,C)\,,
\end{equation}

and (ii) is odd under the reflection

\begin{equation}
\label{eq:reflection}
\calR: \quad \Sigma(\psi;A,B,C) = -\Sigma(\psi;B,A,C) \,,
\end{equation}

and (iii) is also odd under the time-reversal

\begin{equation}
\label{eq:time reversal}
\calT: \quad \Sigma(\psi;A,B,C) = - \Sigma(\calT\psi; A,B,C)   \,.
\end{equation}

Furthermore, (iv) it is topological, in the sense that:

\begin{equation}
\label{eq:topological}
\Sigma(\psi;A,B,C) = \Sigma(\psi;A',B',C')\,,
\end{equation}

where $A',\,B',\,C'$ are smoothly deformed subregions with the same topology as $A, \, B,\,C$.
Finally, (v) it is universal

\begin{equation}
\label{eq:deform Hamiltonian}
\Sigma(\psi;A,B,C) = \Sigma(\psi';A,B,C)\,,
\end{equation}

where $\psi'$ is equal to $\psi$ deformed by a $U(1)$ symmetric local operator.

The justification of these properties uses three ingredients.
The first one is a conversion formula that follows the Schmidt decomposition

\begin{equation}
\label{eq:modular duality}
K_D \ket{\psi} = K_{\bar D} \ket{\psi}\,,
\end{equation}

where $K_D$ and $K_{\bar D}$ are the modular Hamiltonians of the state $\ket{\psi}$ for an arbitrary region $D$ and the complement $\bar{D}$, respectively.
The second is the clustering property, the correlation function of two operators $O_x$ and $O_y$ factorizes when their distance $|x-y|$ is larger than the correlation length $\xi$

\begin{equation*}
\braket{O_x O_y}_\psi = \braket{O_x}_\psi \braket{O_y}_\psi + \calO(e^{-|x-y|/\xi})\,,
\end{equation*}

where $\braket{\cdot}_\psi$ is the expectation value in the state $\ket{\psi}$.
The third one is an assumption on the decomposition of modular Hamiltonians when they act on the state:

\begin{equation*}
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\draw (0,0) -- (2,0) -- (2,0.8) -- (0,0.8) -- (0,0);
\draw (0.65,0) -- (0.65,0.8);
\draw (2-0.65,0) -- (2-0.65,0.8);
\node at (0.325,0.4) {\scriptsize $A$};
\node at (0.975,0.4) {\scriptsize $B$};
\node at (2-0.325,0.4) {\scriptsize $C$};
\end{tikzpicture}\quad
\begin{aligned}
(K_{AB} & + K_{BC}) \ket{\psi}  \\
&\approx (K_{ABC} + K_B)\ket{\psi}\,,
\end{aligned}
\end{equation*}

where $A$ and $C$ do not meet directly <cit.>.
We expect the decomposition to hold when the linear sizes of the subregions are larger than the correlation length.
The proof of the $CRT$ transformation rules only uses eq:modular duality.
The argument of the topological rigidity and universal nature uses the other two ingredients as well. The third assumption was also used in <cit.> to argue the rigidity of the modular commutator formula for the chiral central charge as well.

Charge conjugation
The charge conjugation sends the charge operator from $Q_D$ to $N_D - Q_D$, where $N_D$ is the maximal number charge of the region $D$. Here, we have assumed a finite Hilbert space per site.
Thus, the right-hand side of eq:charge conjugation is essentially

\begin{equation*}
\begin{aligned}
\langle[K_{AB}, & (N_{BC} - Q_{BC})^2] \rangle_\psi \\
\braket{[K_{AB}, -2 N_{BC}Q_{BC} + Q_{BC}^2]}_\psi\,,
\end{aligned}
\end{equation*}

which can be shown to equal to the left-hand side by using $\braket{[K_{AB}, Q_{BC}]}_\psi = 0$ (see app:first law).

For a symmetric choice of $A,B$ and $C$, reflection effectively interchanges $A$ and $B$.
In fact, eq:reflection holds more generally even when $A,B,C$ have unequal sizes and different shapes.
Nevertheless, we still call eq:reflection the reflection transformation. Proving it is equivalent to showing the vanishing of the following quantity

\begin{aligned}
\langle[K_{AB},& Q_{BC}^2 + Q_{AC}^2] \rangle_\psi \\
& =\langle[K_{AB}, Q_{BC}^2 + Q_{AC}^2 - Q_{ABC}^2]\rangle_\psi\,,
\end{aligned}

where we have used $\braket{[K_{AB}, Q_{ABC}^2]}_\psi=0$ in the second step, i.e., modular flow in $AB$ does not change the charge in $ABC$.
Expanding the charge operators in terms of $Q_A,Q_B$ and $Q_C$ gives

\langle[K_{AB}, Q_{C}^2 - 2 Q_{A} Q_{B}]\rangle_\psi \stackrel{\text{Eq.}\eqref{eq:modular duality}}{=} 0,

which finishes the proof.

Time reversal
The time reversal keeps the charge operator invariant. eq:time reversal follows from $\braket{\calT\psi| \calT \phi} = \braket{\phi | \psi}$ and the hermiticity of $K_AB$ and $Q_BC$.

\begin{figure}
\begin{tikzpicture}[scale=0.9,baseline={(current bounding box.center)}]
\draw[dashed, gray] (-1.5,-1.5) -- (-1.5,1.5) -- (1.5,1.5) -- (1.5,-1.5) -- (-1.5,-1.5);
\filldraw[fill=red!10] (-0.7,0.65) circle (0.3);
\node at (-0.8,0.75) {\scriptsize $b$};
\filldraw[fill=white] (0,0) circle (1);
\draw (0,0) -- (0,1);
\draw (0,0) -- (0.85,-0.5);
\draw (0,0) -- (-0.85,-0.5);
\node at (1.1,-1.2) {\scriptsize $D$};
\node at (0.4,0.4) {\scriptsize $A$};
\node at (-0.4,0.4) {\scriptsize $B$};
\node at (0,-0.55) {\scriptsize $C$};
\node at (-1.2,1.2) {\scriptsize (a)};
\end{tikzpicture}\hspace{15pt}
\begin{tikzpicture}[scale=0.9,baseline={(current bounding box.center)}]
\draw[dashed, gray] (-1.5,-1.5) -- (-1.5,1.5) -- (1.5,1.5) -- (1.5,-1.5) -- (-1.5,-1.5);
\filldraw[fill=red!10] (-0.8,-0.6) circle (0.3);
\node at (-0.9,-0.7) {\scriptsize $b$};
\filldraw[fill=white] (0,0) circle (1);
\draw (0,0) -- (0,1);
\draw (0,0) -- (0.85,-0.5);
\draw (0,0) -- (-0.85,-0.5);
\node at (1.1,-1.2) {\scriptsize $D$};
\node at (0.4,0.4) {\scriptsize $A$};
\node at (-0.4,0.4) {\scriptsize $B$};
\node at (0,-0.55) {\scriptsize $C$};
\node at (-1.2,1.2) {\scriptsize (b)};
\end{tikzpicture}
\caption{Geometric deformation. Undeformed regions A,B,C with: (a) Deformation "b" is away from the triple contact points. (b) Deformation is close to one triple contact point.}
\label{fig:geometric deformation}
\end{figure}

\paragraph{Topological}
We first deform the boundary between regions far from any triple contact point.
One typical example is \figref{fig:geometric deformation}~(a), where a small blob $b$ is removed from $D$ and attached to $B$. We write the change of $Σ$ as the sum of three terms
\begin{equation*}
%\label{eq:geometric deformation}
\begin{aligned}
2\i\Delta\Sigma=& \braket{[K_{ABb}, Q_b^2]}_\psi \\
& + \braket{[K_{ABb}, 2Q_b Q_{BC}]}_\psi \\
& + \braket{[K_{ABb} - K_{AB}, Q_{BC}^2]}_\psi\,,
\end{aligned}
\end{equation*}
and show that each vanishes individually.
We use \eqnref{eq:modular duality} to replace $K_ABb$ with $K_CD∖b$, the modular Hamiltonian on the complement region.
The vanishing of the first term becomes clear.
The vanishing of the second term requires the clustering property, i.e.,
\begin{aligned}
\braket{[K_{CD\setminus b},Q_b Q_{BC}]}_\psi
\approx  \braket{Q_b}_\psi \braket{[K_{CD\setminus b}, Q_{BC}]}_\psi = 0
\end{aligned}
where we used that the support of $[K_CD∖b, Q_BC]$ is far from the blob $b$.
The vanishing of the third term follows from the decomposition $K_CD∖b →K_CD + K_D∖b - K_D$, where we have to assume that the blob $b$ is larger than the correlation length.

Next, consider deformations that are close to one of the triple contact points. One typical example is shown in \figref{fig:geometric deformation}~(b).
We write the change of $Σ$ as
\begin{equation*}
\begin{aligned}
2\i\Delta\Sigma = & \langle [K_{ABb}, Q_{BbC}^2] \rangle - \langle [K_{AB}, Q_{BC}^2] \rangle \\
=& \langle [K_{ABb}, Q_{BbC}^2 - Q_{ABbC}^2] \rangle - \langle [K_{AB}, Q_{BC}^2] \rangle
\end{aligned}
\end{equation*}
where we imitate the proof of reflection transformation to add $Q_ABbC^2$ in the second step. The vanishing of the above expression then follows the decomposition $K_ABb →K_AB+K_Bb - K_B$.

\paragraph{Universal} We deform the state $|ψ⟩$ to $|ψ_x'⟩ ∝|ψ⟩+ O_x |ψ⟩$ by a $U(1)$ symmetric local operator at position $x$, and show
\eqnref{eq:deform Hamiltonian}.
If $O_x$ is inside a region and far from any boundary, the clustering property implies that the reduced density matrices and thus modular Hamiltonians of other regions do not change.
One can then prove \eqnref{eq:deform Hamiltonian} by using \eqnref{eq:modular duality} and the clustering property repetitively (see \appref{app:rigidity}).
The case where $O_x$ is close to a boundary can be justified by exploiting the topological rigidity of $Σ(ψ;A,B,C)$, similar to the argument in \cite{Kitaev:2005dm}.

Although the results thus far are compatible with identifying $Σ∝σ_xy$ in \eqnref{eq:conjecture}, the constant of proportionality must be determined. In particular we need to eliminate the trivial possibility that \eqref{eq:conjecture} always vanishes. We turn to this task next.

\emph{Connection to Hall conductance}---
The universal property \eqnref{eq:deform Hamiltonian} allows us to deform the state to a nicer form.
Therefore, without loss of generality, we can consider gapped states whose edge can be described by a (1+1)D conformal field theory (CFT) with the central charge $c$ and $c̅$.
The modular Hamiltonian of a simply connected region $D$ with a smooth boundary can be approximated by the same CFT Hamiltonian supported on the boundary $∂D$~\cite{Kitaev:2005dm,HaldaneLi}
\begin{equation}
\label{eq:KD of 2d bulk}
\rho_D = \frac{1}{Z_D} e^{-K_D}\,,\quad
K_D = \beta H_{CFT}\,,
\end{equation}
where $β$ is a non-universal constant. $H_CFT$ is defined by the Virasoro generators $L_0$ and $L̃_0$
\begin{equation*}
H_{CFT} = \frac{2\pi}{\ell_D} (L_0 + \tilde{L}_0 - \frac{c+\bar{c}}{24})\,,
\end{equation*}
where $ℓ_D$ is the circumference of the region $D$.
Thus, $K_D$ only generates non-trivial modular flow along $∂D$, where the excitation moves with the speed $v=β$. In the case with $U(1)$ symmetry, the CFT is augmented by a holomorphic and anti-holomorphic current $J,J̃$ at level $k_L, k_R$. The Hall conductance is $σ_xy = (k_L-k_R)/2π$.

It is more instructive to consider the \emph{$U(1)$ defect operator} $e^ıμQ_D$ and obtain $Q_D^2$ via Taylor expansion in $μ$.
This operator creates a line defect along the boundary of $D$, and thus its expectation value $Z_D(μ) := ⟨ψ|e^ıμQ_D|ψ|$⟩ satisfies an area-law

\begin{equation}
\label{eq:area law}
\ln Z_{D}(\mu) = -\alpha \ell_{D} + \gamma + \cdots
\end{equation}

where $\alpha$ is the line tension, $\gamma$ is universal, and the ellipsis represents terms that vanish in the limit $\ell_D \rightarrow \infty$ <cit.>. Our focus is the area-law coefficient, which is dictated by the charged Cardy formula <cit.>

\begin{equation}
\label{eq:area law coefficient}
\alpha = \frac{(k_L+k_R)\mu^2}{4\pi \beta}\,.
\end{equation}

This result is related to the spectral flow of the CFT ground state energy under a $U(1)$ twist <cit.>.
Its linearity in $k_L+k_R$ is because the levels physically quantify the number of charged modes. The $\mu^2$ dependence will play an important role in our following derivation.

We reinterpret the result with a quasi-particle picture. There are chiral and anti-chiral modes sitting near the boundary $\partial D$, the numbers of which are proportional to $k_L$ and $k_R$, respectively. Their correlation across the boundary contributes to the line tension. We can then separate $\alpha$ into two terms that account for the contribution from the two types of modes

\begin{equation}
\label{eq:line tension}
\alpha_{\text{chiral}} = \frac{k_L \mu^2}{4\pi \beta}\,,\quad
\alpha_{\text{anti-chiral}} = \frac{k_R \mu^2}{4\pi \beta}\,.
\end{equation}

The idea is sketched in fig:ABC partition (a), where the chiral and anti-chiral modes are the red and blue dots respectively. Under the modular flow generated by $K_D$, the chiral and anti-chiral modes inside the region $D$ move in opposite directions, as shown by the arrows.

(a) Red and blue dots are the chiral and anti-chiral modes, the bonds represents the correlation. The arrows designate their motion under the modular flow $U_D(s)$. (b) The anti-chiral degrees of freedom are suppressed for the clarity of the figure. The arrows designate the motion under the modular flow $U_{AB}(s)$. The yellow shaded disks emphasis modes that make contribution to the change of $Z_{BC}(\mu;s)$ at the linear order in $s$.

Corresponding to eq:sigmaxy formula 1, we apply this picture to compute the response of $e^{\i\mu Q_{BC}}$ under the modular flow generated by $K_{AB}$

\begin{equation*}
Z_{BC}(\mu;s):=\braket{\psi|e^{\i sK_{AB}} e^{\i\mu Q_{BC}} e^{-\i sK_{AB}}|\psi}\,.
\end{equation*}

It suffices to focus on the chiral modes, the contribution from the anti-chiral ones is the opposite.
See fig:ABC partition (b). Only the motion of chiral modes that are near the $ABC$ and $ABD$ triple-contact points can affect $Z_{BC}(\mu;s)$ at the linear order in $s$. Near the $ABC$ triple-contact point, chiral modes in $B$ that are correlated with $C$ will enter $A$ and increase the line tension. The $ABD$ triple-contact point can be analyzed similarly. After summing up the two contributions, we have

\begin{equation*}
\frac{d}{ds} \ln Z_{BC}(\mu;s) \big|_{s=0} = -2(\alpha_{\text{chiral}} - \alpha_{\text{anti-chiral}}) v
\end{equation*}

where the velocity is $v=\beta$ as introduced before.
By plugging in eq:line tension, we have

\begin{equation}
\label{eq:2+1d result}
\frac{d}{ds} \ln Z_{BC}(\mu;s) \big|_{s=0} = -\frac{k_L-k_R}{2\pi} \mu^2 = -\sigma_{xy} \mu^2\,,
\end{equation}

Here the non-universal factor $\beta$ is completely canceled, which yields a universal answer. This is the main result of this work.
Expanding $e^{\i\mu Q_{BC}}$ on the left-hand side above to quadratic order in $\mu$ gives the right-hand side of eq:sigmaxy formula 1 which verifies our identification in eq:conjecture.
In app:anomaly, we apply the same idea to purely 1+1D CFTs and show that it detects the perturbative chiral anomaly, which is consistent with the spirit of the bulk-edge correspondence.

(a) $\Sigma(\psi,A,B,C)$ across the phase diagram. Blue, red and green curves correspond to the disk geometry at the corresponding radii $r$. The dashed line is the value of the quantized Hall conductance. (b) Finite size scaling. The blue dots are the data points and the red line is the fitting result. We choose $L_x=L_y=40$. For the annulus geometry, the inner and outer radius is $7$ and $14$. For (b), the inverse correlation length is $\xi^{-1} = 0.398$.

We provide numerical support to our conjecture by simulating free fermion systems. In this section we report on a lattice model and focuses on $\Sigma(\psi;A,B,C)$ in different phases, its finite size scaling and topological rigidity.
app:details of numerics contains other aspects as well as results on a continuum model.
In particular, we numerically check the remarkable $\mu^2$ dependence in the area-law coefficient eq:area law coefficient in lattice systems and find excellent agreement.

We simulate the $\pi$-flux model on a square lattice. The Hamiltonian, $H = \sum_{ij}t_{ij} c_i^\dag c_j$, consists of nearest-neighbor hopping terms that implement $\pi$-flux per plaquette and next-nearest-neighbor hoppings that open up a gap and give rise to various phases.
We choose a 1-dimensional path inside the entire phase diagram that is parameterized by $\eta\in\bbR$.
The system is time-reversal symmetric for $\eta\leq 0$ and breaks the symmetry for $\eta>0$.
There is a topological phase at $0<\eta<2$ with a Chern number $C=1$ and two trivial phases elsewhere. In the numerics, we also add weak on-site disorders to make the situation more general.

The result of $\Sigma(\psi;A,B,C)$ across the phase diagram is shown in fig:phase transition (a).
When $A,B,C$ form a complete disk, as required by the definition of the formula, we obtain the correct Hall conductance of each phase.
Our formula indeed vanishes identically in the time-reversal symmetric region $\eta<0$.
Near the two phase boundaries, the curve becomes steeper for larger subsystem sizes.
In the topological phase, we perform a more detailed finite size scaling. Let $\Delta := 1- 2\pi \Sigma(\psi;A,B,C)$ be the difference between the theoretical and numerical value, and the result is shown in fig:phase transition (b). The deviation decays exponentially with the linear size of the subsystem and the decaying exponent is twice of the correlation length.
As a sanity check, if $A,B,C$ form an annulus (replace a disk neighborhood of the $ABC$ triple-contact point by region $D$) we find a null result. This is consistent with our quasi-particle picture.

We also provide numerical evidence on the topological rigidity as independent support from our analytical argument. For example, we deform the boundaries near the triple-contact points, depicted in fig:topological rigidity (a1)-(a3), and compute $\Sigma(\psi)$ for the deformed geometry.
As shown in fig:topological rigidity (b), its value remain the same across the entire phase diagram even though the added blobs are not much larger than the correlation length.
The results are the same for other deformations.

Topological rigidity against deformation near the triple-contact points. (a1)-(a3) Attach a blob to region $A,B$ and $C$ near the triple-contact point, respectively. (b) Evaluate $\Sigma(\psi,A,B,C)$ for the deformed geometry across the phase diagram. We choose $L_x=L_y=30$ and the radius of the disk $ABC$ to be $r=10$. The blob has a linear size of $4$.

Discussion—The present results can be understood as linear response extensions that go beyond `static' entanglement features by utilizing the dynamics generated by the modular (or entanglement) Hamiltonian. Here, the charged Cardy formula evaluates the `static' features of the $U(1)$ twist operator $e^{i\mu Q}$ <cit.>, and its linear response under the modular flow is shown to give the Hall conductance.
This perspective also suggests new potential connections between conventional responses and quantum entanglement. For instance, in free fermion systems with single-charge fermion, we may anticipate an entanglement version of the Wiedemann–Franz (WF) law, which indeed can be obtained by combining the result in <cit.> and our formula:

\begin{equation}
\braket{[K_{AB}, K_{BC} - \frac{\pi^2}{3} Q_{BC}^2]}_\psi = 0\,.
\end{equation}

Understanding the violation of this entanglement WF law and finding other such relations can shed new light on the study of entanglement of many-body systems.

Entanglement entropy suffers from UV divergence in generic Lorentz invariant quantum field theories (QFTs). Distilling universal information from entanglement requires one to carefully define  entanglement in QFTs and further distill the universal part of it.
Our formula, shown to be universal, is an interesting quantity to compute especially in topological quantum field theories.
Formulating the calculation in a meaningful way requires a better understanding of multipartite entanglement in QFTs, which is an interesting and nontrivial question in its own right [We thank Daniel Jafferis and Xi Yin for discussions on this point.].

Acknowledgements—We thank John Cardy, Yiming Chen, Meng Cheng, Joel Moore, Huajia Wang, Xi Yin, Carolyn Zhang, and especially Yingfei Gu, Daniel Jafferis, Ari Turner and Michael Zaletel for many helpful discussions and comments.

RF and AV are supported by a Simons Investigator award (AV) and by the Simons Collaboration on Ultra-Quantum Matter, which is a grant from the Simons Foundation (651440, AV).

RF is supported by the DARPA DRINQS program (award D18AC00033).

This work is funded in part by a QuantEmX grant from ICAM and the Gordon and Betty Moore Foundation through Grant GBMF9616 to Ruihua Fan.

RS acknowledges support by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Department of
Energy Computational Science Graduate Fellowship under Award Number DESC0022158.

§ FURTHER DETAILS IN THE PROOF OF GENERAL PROPERTIES

§.§ Expectation value of the commutator

In the main text, we have repeatedly used the following equality in the proof of general properties of the formula $\Sigma(\psi; A,B,C)$

\begin{equation}
\label{eq:commutator}
\braket{\psi|[K_{AB}, Q_{BC}]|\psi} = 0
\end{equation}

where $A,B,C$ are three disjoint spatial subregions, $K_{AB}$ and $Q_{BC}$ are the modular Hamiltonian and total charge operator for the associated region, respectively.
In this section, we explain eq:commutator by doing direct manipulations.

Note that $Q_{BC} = Q_B + Q_C$, and we can break the left-hand side of eq:commutator into two pieces

\braket{\psi|[K_{AB}, Q_{BC}]|\psi} = \braket{\psi|[K_{AB}, Q_{B}]|\psi} + \braket{\psi|[K_{AB}, Q_{C}]|\psi}\,.

The second piece vanishes because $K_{AB}$ and $Q_C$ are supported on different subregions.
The vanishing of the first term follows from the conversion formula eq:modular duality.
Let $D$ denote the complement of the subregions $A,B$ and $C$, then we have $K_{AB}\ket{\psi} = K_{CD} \ket{\psi}$ and

\braket{\psi|[K_{AB}, Q_{B}]|\psi} = \braket{\psi|[K_{CD}, Q_{B}]|\psi}\,,

which obviously vanishes.
More generally, for any operator $O$, we can have

\begin{equation}
\braket{\psi|[K_{AB},O]|\psi} = 0\,,
\end{equation}

whenever $[K_{CD},O]=0$ even though $[K_{AB},O]$ does not have to be zero.

§.§ Rigidity

We fill in details on the rigidity of $\Sigma(\psi;A,B,C)$ when the Hamiltonian, or equivalently the state, is locally deformed. We change the state $\psi$ by a local operator $O_x$ supported near $x$, which preserves the $U(1)$ symmetry $[O_x,Q]=0$. For simplicity, we assume that the operator has a zero expectation value and write the deformed state as

\begin{equation}
\ket\psi \mapsto \ket{\psi_x'} = \frac{1}{\sqrt{Z}}\left(\ket{\psi} + O_x \ket{\psi} \right)
\end{equation}

where $Z = 1 + \braket{O_x^\dag O_x}_\psi$ is the normalization factor of the wavefunction.
The equality we want to show is

\begin{equation}
\label{eq:deform Hamiltonian2}
\Sigma(\psi;A,B,C) = \Sigma(\psi_x';A,B,C)
\end{equation}

By the clustering property, the correlation functions of local opeartors that are far separated from $x$ are not changed, so are the reduced density matrices and modular Hamiltonians.

Suppose the operator $O_x$ is inside the bulk of region $C$ and is far from any boundary. The right-hand side is

\begin{equation*}
\begin{aligned}
\Sigma(\psi_x';A,B,C) =& \frac{-i}{2} \braket{\psi'|[K_{AB}', Q_{BC}^2]|\psi'} \\
\approx & \frac{-i}{2Z} \braket{\psi|(1+O_x^\dag)[K_{AB}, Q_{BC}^2](1+O_x)|\psi}
\end{aligned}
\end{equation*}

Note that the commutator $[K_{AB}, Q_{BC}^2]$ is supported in region $B$ and is far separated from $x$ so that we can use the clustering property, it follows from which that eq:deform Hamiltonian2 is true.

When the operator $O_x$ is inside the bulk of region $B$, we can apply the conversion formula for the modular Hamiltonian to first convert $K_{AB}'$ to $K_{CD}$ and repeat the above argument.

§ RELATION TO THE U(1) CHIRAL ANOMALY

A nonzero Hall conductance in the (2+1)D bulk corresponds to the perturbative $U(1)$ chiral anomaly on the (1+1)D edge.
It is natural to expect eq:sigmaxy formula 1 to be able to detect the anomaly in a pure (1+1)D CFT. We confirm this intuition by imitating the analysis of the defect operator in the main text.

A generic (1+1)D CFT with a global $U(1)$ symmetry is characterized by a holomorphic and anti-holomorphic current $J, \tilde{J}$. The $U(1)$ charge density $\rho$ is given by $2\pi\rho(x) = J(x) + \tilde J(\bar{x})$.
The two currents are spin-1 primary fields with the following operator-product expansion on the complex $z$-plane <cit.>

\begin{equation}
\begin{aligned}
J(z) J(0) \sim & \frac{k_L}{z^2}\,,\quad \tilde J(\bar{z}) \tilde J(0) \sim  \frac{k_R}{\bar{z}^2} \,.
\end{aligned}
\end{equation}

where $k_{L}$ and $k_R$ are called the levels. When $k_L \neq k_R$, the system has a perturbative chiral anomaly and cannot be defined on manifolds with a boundary.

Consider a system on a circle of length $L=2\pi$. We assume that the conformal vacuum $\ket{0}$ can be realized as the ground state of a critical lattice system and has a Schmidt decomposition.
Let $A,B,C$ be three adjacent disjoint intervals shown below.
By conformal invariance, the modular Hamiltonian $K_{AB}$ is a local integral of the stress-energy tensor

\begin{equation}
\label{eq:KD CFT}
\begin{tikzpicture}[scale=0.8,baseline={(current bounding box.center)}]
\draw[gray!50,thick] (0,0) circle (1);
\draw[fill=black] (1,0) circle (0.05) node[right] {\scriptsize$x'=0$ \text{or} $L$};
%\node[below left] at (1,0) {$x'=L$};
\draw[fill=black] (-0.174,0.984) circle (0.05) node[above] {\scriptsize$x'=y$};
\draw[fill=black] (-0.985,-0.174) circle (0.05) node[left] {\scriptsize$x'=\ell$};
\draw[fill=black] (0.1736,-0.9848) circle (0.05) node[below] {\scriptsize$x'=x$};
\draw[<->,>=stealth,line width=1.5,Acolor!90!black] (1,0) arc (0:98:1 and 1) node at (0.9,0.9) {$A$};
\draw[<->,>=stealth,line width=1.5,Bcolor!90!black] (-0.19,0.984) arc (100:188:1 and 1) node at (-0.95,0.8) {$B$};
\draw[<->,>=stealth,line width=1.5,Ccolor!90!black] (-0.985,-0.185) arc (190:278:1 and 1) node at (-0.85,-0.95) {$C$};
\end{tikzpicture}\qquad
K_{AB} = 4\pi \int_{AB} \frac{\sin \frac{\ell-x}{2} \sin \frac{x}{\pi}}{\sin \frac{\ell}{2}}\, T_{00}(x) dx\,,
\end{equation}

where $T_{00} = \frac{1}{2\pi} (T+\tilde T)$ is the energy density  <cit.>.
Under the modular flow generated by $K_{AB}$, the evolution of $e^{i\mu Q_{BC}}$ can be computed by replacing the defect operator by two twist operators $V_{\pm\mu}$ at the end points, i.e.

\begin{equation}
e^{i\mu Q_{BC}} \mapsto \delta^{(k_L+k_R)(\frac{\mu}{2\pi})^2} V_\mu(x) V_{-\mu}(y)\,,
\end{equation}

where $\delta$ is a UV length scale, $V_{\pm\mu}$ are primary operators with the conformal dimension $h = \frac{k_L}{2} \left(\frac{\mu}{2\pi} \right)^2$, $\tilde h = \frac{k_R}{2} \left(\frac{\mu}{2\pi} \right)^2$ and are not mutually local with the charged operators <cit.>.
The linear response is given by

\begin{equation}
\begin{aligned}
\frac{d}{ds} \log & Z_{BC}(\mu,s) \Big|_{s=0} = -\mu^2\frac{k_L-k_R}{4\pi} \left( 1 - 2 \frac{\sin \frac{x}{2} \sin \frac{\ell-y}{2}}{\sin \frac{\ell}{2} \sin \frac{x-y}{2}} \right) \,.
\end{aligned}
\end{equation}

where the ratio of the sines is the cross ratio of the four points drawn on the circle.
The appearance of $k_L-k_R$ is due to the fact that the chiral and anti-chiral modes move in the opposite direction under the modular flow.
Now, we consider the limit $x\rightarrow 2\pi$ where $ABC$ occupies the entire circle and their union is in the pure state $\rho_{ABC} = \ket{0}\bra{0}$.
On the one hand, the above result reduces to

\begin{equation}
\label{eq:cft result}
\frac{d}{ds} \log Z_{BC}(\mu;s) \big|_{s=0} = -\frac{k_L-k_R}{4\pi} \mu^2\,.
\end{equation}

On the other hand, we have $K_{AB}\ket{\psi} = K_C\ket{\psi}$. Since $K_{C}$ preserves the total charge in $BC$, $Z_{BC}(\mu;s)$ should not have any time-dependence, which implies $k_L=k_R$ for a well-defined quantum state. An intriguing fact is that eq:2+1d result differs from eq:cft result by a factor of $2$, which comes from summing over contribution from the two triple-contact point.

§ MORE ON THE FREE-FERMION NUMERICS

In this appendix, we provide more details on the free-fermion numerics.
In addition to the definition of the $\pi$-flux model, we present more numerical results.
We also show results for a continuum model, the integer quantum Hall system at $\nu=1$, which serves as a complement to numerics on the lattice.

§.§ $\pi$-flux model

Definition of the model:  The $\pi$-flux model is defined on a square lattice as follows. There are two sublattices A and B, which are depicted as the red and blue dots in the figure below.
Each unit cell contains an $A$ and $B$ sublattice site, as is denoted by the circle in the figure.
We use $r$ to label the unit cell, $c_{r,A}$ and $c_{r,B}$ are the fermion annihilation operators.
The Hamiltonian consists of two parts

\begin{equation}
H = H_1 + H_2\,.
\end{equation}

The first part contains the nearest-neighbor hopping terms that carry appropriate phases such that fermions can pick up a $\pi$ flux after hopping around each plaquette. Our convention for the phases is analogous to the Landau gauge. As is depicted in the figure, the solid links represent hopping with a positive amplitude $t_1>0$ and the dashed links represent hopping with a negative amplitude $-t_1$, i.e.,

\begin{equation}
\begin{tikzpicture}[scale=0.6, baseline={(current bounding box.center)}]
\foreach \x in {0,1,...,5}{
\draw (\x,-0.3) -- (\x,3.3);
\foreach \y in {0,2}{
\draw (-0.3,\y) -- (5.3,\y);
\foreach \y in {1,3}{
\draw[dashed] (-0.3,\y) -- (5.3,\y);
\foreach \y in {0,2}{
\foreach \x in {0,1,...,5}{
\filldraw[red!65] (\x,\y) circle (0.05);
\foreach \y in {1,3}{
\foreach \x in {0,1,...,5}{
\filldraw[blue!65] (\x,\y) circle (0.05);
\foreach \y in {0,1,2}{
\foreach \x in {0,1,...,4}{
\node at (\x+0.5,\y+0.5) {\scriptsize $\pi$};
\draw (1.2,0.5) arc (0:360:0.2 and 0.8);
\node[red!65,left] at (-0.2,0) {\scriptsize$A$};
\node[blue!65,left] at (-0.2,1) {\scriptsize$B$};
\node[red!65,left] at (-0.2,2) {\scriptsize$A$};
\node[blue!65,left] at (-0.2,3) {\scriptsize$B$};
\end{tikzpicture}\qquad
\begin{aligned}
H_1 = t_1 \sum_{r} & c_{r,A}^\dag c_{r+\hat{x},A} - c_{r,B}^\dag c_{r+\hat{x},B} \\
& + c_{r,A}^\dag c_{r,B} + c_{r,B}^\dag c_{r-\hat{y},A} + h.c.
\end{aligned}
\end{equation}

where $\hat{x}$ and $\hat{y}$ denote the horizontal and vertical directions.
The second part contains a special next-nearest-neighbor hopping term with amplitude $t_2$, which is shown by the arrow below

\begin{equation}
\begin{tikzpicture}[scale=0.6, baseline={(current bounding box.center)}]
\foreach \x in {0,1,...,5}{
\draw (\x,-0.3) -- (\x,3.3);
\foreach \y in {0,2}{
\draw (-0.3,\y) -- (5.3,\y);
\foreach \y in {1,3}{
\draw[dashed] (-0.3,\y) -- (5.3,\y);
\foreach \y in {0,2}{
\foreach \x in {0,1,...,5}{
\filldraw[red!65] (\x,\y) circle (0.05);
\foreach \y in {1,3}{
\foreach \x in {0,1,...,5}{
\filldraw[blue!65] (\x,\y) circle (0.05);
\foreach \y in {0,2}{
\foreach \x in {0,1,...,4}{
\draw[mid arrow,black!70] (\x,\y) -- (\x+1,\y+1);

\node[red!65,left] at (-0.2,0) {\scriptsize$A$};
\node[blue!65,left] at (-0.2,1) {\scriptsize$B$};
\node[red!65,left] at (-0.2,2) {\scriptsize$A$};
\node[blue!65,left] at (-0.2,3) {\scriptsize$B$};
\end{tikzpicture}\qquad
H_2 = \sum_{r} t_2\, c_{r+\hat{x},B}^\dag c_{r,A} + h.c.\,.
\end{equation}

The system is gappless with two Dirac cones if we turn off the next-nearest neighbor hopping, i.e., $t_2=0$. The system is generally gapped when $t_2\neq 0$.

Without loss of generality, we fix $t_1=1$.
When $t_2$ is purely real, the system is time-reversal symmetric and cannot have a nonzero Hall conductance.
When $t_2$ is purely imaginary, the system breaks the time-reversal symmetry and can become non-trivial. Let $t_2=i|t_2|$, then one can show that the system is a Chern insulator with a Chern number $C=1$ for $0<|t_2|<2$ and a trivial insulator for $|t_2|>2$. In the following, we choose a path in the 2-dimensional complex $t_2$ plane, that is parameterized by $\eta$ such that we have $t_2=|\eta|$ being purely real for $\eta<0$ and $t_2=i|\eta|$ purely imaginary for $\eta>0$. This explains the phase diagram quoted in the main text.

Modular flow of the defect operator
We derive our conjecture by studying the linear response of the defect operator $e^{i\mu Q_{BC}}$ under the modular flow generated by $K_{AB}$.
We did the analysis in the topological limit, i.e., the state is well-captured by a topological quantum field theory, and found that only the $\mu^2$ term appears in the line tension and final expression.
Here, we provide numerical result and provide an independent support to this $\mu^2$ dependence.

For brevity, we denote the linear response of the defect operator by

\begin{equation}
F_{ABC}(\mu) := \frac{d}{ds} \ln Z_{BC}(\mu;s) \big|_{s=0} = \frac{i\braket{[K_{AB}, e^{i\mu Q_{BC}}]}_\psi}{\braket{e^{i\mu Q_{BC}}}_\psi}
\end{equation}

and numerically study its dependence on $\mu$. The results are shown in fig:modular flow defect operator.
fig:modular flow defect operator (a) show that the real part $\Re F_{ABC}(\mu)$ is almost quadratic in $\mu$.
Because the line tension $\alpha$ increases with $\mu$, the absolute value of $\braket{e^{i\mu Q_{BC}}}_\psi$ becomes extremely small near $\mu=\pi$ and causes numerical instability. This may explain the wiggles near $\mu=\pi$.
We then study the discrepancy between $\Re F_{ABC}(\mu)/\mu^2$ and the Hall conductance and found that it decays exponentially with respect to the subsystem size, which is shown in fig:modular flow defect operator (b).
In the lattice system, $F_{ABC}(\mu)$ also acquires a small imaginary part, which has a cubic dependence in $\mu$ when $\mu$ is small and the result is shown in fig:modular flow defect operator (c). This can be understood via a Taylor expansion of $F_{ABC}(\mu)$ and one can find that the non-vanishing leading term of the imaginary part is cubic in $\mu$.

Linear response of the defect operator under the modular flow. (a) The real part of $F_{ABC}(\mu)$ is almost quadratic in $\mu$ (b) Finite size scaling of the discrepancy $\delta = \frac{1}{2\pi} - \frac{\Re F_{ABC}(\mu)}{\mu^2}$ (c) The imaginary part of $F_{ABC}(\mu)$ is cubic in $\mu$ for small $\mu$. Figure (a) and (c) uses $L_x=L_y=26$ and $r=7$ as the radius of the subsystem. Figure (b) uses $L_x=L_y=30$, $t_2 = i$ and $\mu=0.1$. Choosing larger values of $\mu$ gives the same result.

Incomplete disk
The definition of $\Sigma(\psi;A,B,C)$ requires that every three subregions meet once and only once at a point. It is crucial to use this topology in order to get $\sigma_{xy}$, which is manifest from the quasi-particle picture. Here, we present numerical results on what happens when this condition is broken.

An example is shown in fig:incomplete disk (a1), where region $A,B$ and $C$ are represented by red, green and blue dots, and region $D$ is the rest part.
We remove part of the region $C$ such that it together with $A$ and $B$ only forms an incomplete disk, which replaces the $ACD$ triple-contact point with an $ABCD$ quadruple-contact point.
Let $\theta_0=2\pi/3$ be the angle of subregion $A$ and $B$, $\Delta \theta$ be the angle of the removed part. It is found that $\Sigma(\psi;A,B,C)$ decreases as $\Delta\theta$ increases, which is shown in fig:incomplete disk (b). This behavior is qualitatively consistent with the quasi-particle picture. Namely, near the quadruple contact point, quasi-particles that have correlation between $B$ and $D$ will make opposite contribution compared with ones that correlates $B$ and $C$.

To further demonstrate the importance of having only triple-contact point, we start with an incomplete disk and gradually add lattice sites to region $A$, as is shown in fig:incomplete disk (a2), where the added sites are near the center of the figure. Once the size of the added blob is comparable with the correlation length, it is expected that $\Sigma$ should return to the quantized value, $1$ in this case. This intuition is confirmed by the result shown in fig:incomplete disk (c).

Evaluate $\Sigma(\psi;A,B,C)$ for an incomplete disk. (a1) We remove part of region $C$ to have an incomplete disk, the angle of the remove part is $\Delta\theta$. (a2) We add lattice sites that are near the center of the system to region $A$ to remove the quadruple-contact point. (b) $\Sigma(\psi;A,B,C)$ decreases with the angle $\Delta\theta$. (c) $\Sigma(\psi;A,B,C)$ quickly returns to the quantized value as the size of the added region increases. We choose the total system size to be $L_x=L_y=40$ and the radius of the disk $ABC$ to be $R=14$.

§.§ Integer quantum Hall system

In this section, we discuss the integer quantum Hall (IQH) system.
For simplicity, we will focus on the lowest Landau level and formulate the problem in the continuum as a complementary to the lattice-model discussion.

The system is defined on a finite disk of radius $R$ with a uniform magnetic field $B$ perpendicular to the plane. We choose the symmetric Coulomb gauge $A_x = -By/2$, $A_y = Bx/2$ and write the single-particle Hamiltonian as ($e>0$)

\begin{equation}
H = \frac{1}{2m} (p_x-eA_x)^2 + \frac{1}{2m} (p_y-eA_y)^2\,.
\end{equation}

The system has an intrinsic length scale $\ell_B = \sqrt{1/eB}$, which describes the typical radius of the semi-classical cyclotron motion.
For simplicity, we will set $\ell_B=1$ in the following discussion.

We will focus on the lowest Landau level (LLL), which has $R^2/2$ degenerate states. For $B>0$, these states are given by the following wavefunctions

\begin{equation}
\label{eq:LLL wavefunction}
\varphi_{n}(z) = \frac{1}{\sqrt{2\pi 2^n n!}} z^n e^{-\frac{|z|^2}{4}},\quad n=0,\cdots,\frac{R^2}{2}-1\,,
\end{equation}

where $z=x+iy$ is the holomorphic coordinate. Replacing $z$ with $\bar z = x-iy$ gives the wavefunction for $B<0$.
Strictly speaking, (<ref>) are solutions on the infinite plane and cannot be applied to a finite disk. We will choose $R\gg 1$ such that (<ref>) is a good approximation to the actual solution away from the boundary.

The results are shown in fig:IQH numerics.
We first consider the case where $ABC$ forms a complete disk, as depicted in (a1).
When the radius of the disk becomes much larger than the magnetic length, $2\pi \Sigma(\psi;A,B,C)$ quickly converges to $1$, as shown in (b).
We then consider the case where $ABC$ forms an incomplete disk, as depicted in (a2), i.e., the case with a quadruple-contact point.
As the angle of the removed region increases, $\Sigma(\psi;A,B,C)$ also decreases, which is qualitative the same as what we have observed for the $\pi$-flux model.
One quantitative difference is that $\Sigma(\psi;A,B,C)$ decreases in a smooth way for the IQH system while quite discontinuously for the $\pi$-flux model. Such a difference may be attributed to using a smooth wavefunction in the continuum verses discrete lattice model.

Evaluate $\Sigma(\psi;A,B,C)$ for the integer quantum Hall system. (a1) $ABC$ form a complete disk. (a2) $ABC$ form an incomplete disk. (b) $\Sigma(\psi;A,B,C)$ as a function of the radius of the disk. (c) $\Sigma(\psi;A,B,C)$ as a function of the angle of the remove region.


….

# X-rays as dominant excitation mechanism of [Fe ii] and $\rm H_{2}$ emission
lines in active galaxies

Oli L. Dors Jr.1, Rogemar A. Riffel2, Mónica V. Cardaci3,4,5, Guillermo F.
Hägele3,4,5, Ângela C. Krabbe1, Enrique Pérez-Montero6, Irapuan Rodrigues1
1 Universidade do Vale do Paraíba, Av. Shishima Hifumi, 2911, Cep 12244-000,
São José dos Campos, SP, Brazil
2 Universidade Federal de Santa Maria, Av. Roraima, 1000, Cep 97105-900, Santa
Maria, Brazil
3 Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET),
Argentina.
4Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de la La
Plata, Paseo del Bosque s/n, 1900 La Plata, Argentina.
5 Departamento de Física Teórica, C-XI, Universidad Autónoma de Madrid, 28049
Madrid, Spain.
6 Instituto de Astrofísica de Andalucía (CSIC), PO Box 3004, 18080 Granada,
Spain<EMAIL_ADDRESS>
###### Abstract

We investigate the excitation mechanisms of near-infrared [Fe ii] and $\rm
H_{2}$ emission lines observed in Active Galactic Nuclei (AGNs). We built a
photoionization model grid considering a two-component continuum, one accounts
for the Big Bump component peaking at $\rm 1Ryd$ and another represents the
X-ray source that dominates the continuum emission at high energies.
Photoionization models considering as ionizing source a spectral energy
distribution obtained from photometric data of the Sy 2 Mrk 1066 taken from
the literature were considered. Results of these models were compared with a
large sample of observational long-slit and Integral field Unit (IFU)
spectroscopy data of the nuclear region for a sample of active objects. We
found that the correlation between the observational [Fe
ii]$\lambda$1.2570$\,\mu$m/Pa$\beta$ vs.
H${}_{2}\lambda$2.1218$\,\mu$m/Br$\gamma$ is well reproduced by our models as
well as the relationships that involve the H2 emission line ratios observed in
the spectroscopic data. We conclude that the heating by X-rays produced by
active nuclei can be considered a common and very important mechanism of
excitation of [Fe ii] and $\rm H_{2}$.

###### keywords:

galaxies: Seyfert – galaxies: ISM – infrared: galaxies

††pagerange: X-rays as dominant excitation mechanism of [Fe ii] and $\rm
H_{2}$ emission lines in active galaxies–References

##  1 Introduction

The excitation of the Narrow Line Region of Seyfert (Sy) galaxies can reveal
how radiation and mass outflows from the nucleus interact with circumnuclear
gas. In particular, near-infrared (hereafter near-IR) observations are a
powerful tool to investigate this issue, because the obscuration –which can
affect the optical morphology of the emitting gas region– is less important at
these wavelengths (Mulchaey et al., 1996; Ferruit et al., 2000). Relevant
emission-lines in the near-IR include [Fe ii] $\lambda\,1.2570\,\mu$m and
$\lambda\,1.6440\,\mu$m, H i lines such as Pa$\beta$, and H2 at
$\lambda\,1.9576\,\mu$m, $\lambda\,2.1218\,\mu$m, and $\lambda\,2.3085\,\mu$m,
which can be used to map the gas kinematics and excitation (e.g. Riffel &
Storchi-Bergmann, 2011b; Riffel et al., 2010). Nevertheless, the dominant
excitation mechanisms of the [Fe ii] and H2 emission lines in the central
regions of active galaxies are still unclear and have been the subject of
several recent studies (e.g. Riffel et al., 2010, 2008, 2006; Storchi-Bergmann
et al., 2009; Hicks et al., 2009; Müller Sánchez et al., 2009; Ramos Almeida
et al., 2009; Rodríguez-Ardila et al., 2005; Rodríguez-Ardila et al., 2004;
Davies et al., 2007).

The H2 can be excited by two mechanisms: (i) fluorescent excitation through
absorption of soft-UV photons (912–1108 Å) in the Lyman and Werner bands,
existing both in star-forming regions and surrounding the Active Galactic
Nuclei (AGNs) (Black & van Dishoeck, 1987) and (ii) collisional excitation due
to the heating of the gas by shocks, the interaction of a radio jet with the
interstellar medium (Hollenbach & McKee, 1989), or by X-ray photons from the
central AGN (Maloney et al., 1996). Several studies based on intensity-line
ratios (Riffel et al., 2010; Storchi-Bergmann et al., 2009; Rodríguez-Ardila
et al., 2005; Rodríguez-Ardila et al., 2004) have shown that collisional
excitation processes dominate the H2 emission surrounding AGNs. However, which
is the dominant process is an open question. Veilleux et al. (1997), using J
and K-band spectra of a sample of 33 Sy 2 galaxies, found that shocks
associated with nuclear outflows are a likely source of both [Fe ii] and $\rm
H_{2}$ emission rather than circumnuclear starbursts, as suggested by Quillen
et al. (1999).

For the [Fe ii] emission, the [Fe ii] $\lambda\,1.2570\,\mu$m/Pa$\beta$ line
ratio is generally used to investigate the main mechanism of excitation. The
value of this line ratio is controlled by the quotient of the volumes of
partially and fully ionized gas regions, with [Fe ii] emission being excited
in the partially ionized gas (Mouri et al., 1990; Mouri et al., 1993;
Rodríguez-Ardila et al., 2005; Riffel et al., 2010, 2008, 2006; Storchi-
Bergmann et al., 2009). Such zones in AGNs are created by X-ray emission (e.g.
Simpson et al., 1996) and/or shock heating of the gas by mass outflows from
the nuclei which interact with the ambient clouds (e.g. Forbes & Ward, 1993).
This problem was addressed by Mouri et al. (2000), who compared the values of
the line ratios [Fe ii] $\lambda\,1.2570\,\mu$m/Pa$\beta$ and [O i]
$\lambda\,6300$ Å/H$\beta$ predicted by models, considering photoionization
and shock heating, with those observed in a sample of AGNs and Starburst
galaxies. These authors pointed out that in AGNs, X-ray heating is the most
important [Fe ii] excitation mechanism. However, Rodríguez-Ardila et al.
(2004), using near-IR spectroscopy of a sample of galaxies obtained with the
Infrared Telescope Facility, found that X-ray excitation is enough to explain
the H2 emission and part of the [Fe ii] emission observed in Sy 1 galaxies,
but fails to explain the emission of these elements in Sy 2. For these
objects, a combination of shocks and circumnuclear star-formation is required
to explain these emissions. Moreover, it is not clear whether the [Fe ii] and
H2 are excited by the same mechanism. Rodríguez-Ardila et al. (2004) found a
correlation between the [Fe ii] $\lambda\,1.2570\,\mu$m/Pa$\beta$ and the H2
$\lambda\,2.1218\,\mu$m/Br$\gamma$ ratios, indicating that both sets of lines
may be originated by a single dominant mechanism. However, high spatial
resolution spectroscopy data from Integral Field Unit (IFU) of active galaxies
indicate that the H2 and the [Fe ii] emitting gas have distinct flux
distributions and kinematics, with the former being considered a tracer of the
feeding of the AGN and the latter a tracer of its feedback (Riffel et al.,
2010, 2009, 2008; Storchi-Bergmann et al., 2009; Hicks et al., 2009; Müller
Sánchez et al., 2009). This result indicates that the lines of these elements
are formed in distinct regions. Although several works have investigated the
excitation origins of H2 and the [Fe ii], it is still unknown whether a common
mechanism can excite these elements. Fortunately, a large number of near-IR
data of AGNs are currently available in the literature, which enables an
extensive comparison with models yielding a more reliable conclusion about the
likely dominant excitation mechanism of these emission lines.

In this paper, we combined near-IR data of Sy galaxies obtained with IFU and
long-slit spectroscopy with photoionization models to investigate the origin
of the H2 and [Fe ii]. In Section 2, we describe the observational data used
in the analysis. The modelling procedures are presented in Sect. 3. In Sect.
4, the diagnostic diagrams used to compare the observational data with our
model predictions are described. Results and discussion are presented in
Sects. 5 and 6, respectively. A conclusion of the outcome is given in Sect. 7.

##  2 Observational data

We compiled from the literature observational data of the nuclear region of
active galaxies in the near-IR and optical spectral range obtained with long-
slit and IFU spectroscopy. The selection criterion was the presence of bright
infrared emission lines in their spectra. These data are described below.

###  2.1 Long-slit data

Near-IR emission line intensity ratios of 35 active galaxies were obtained
from Rodríguez-Ardila et al. (2004), Reunanen et al. (2002), Knop et al.
(2001), and Riffel et al. (2006). This sample comprises long-slit data of 13
Sy 1 and 21 Sy 2 galaxies, along with 1 Quasar. The intensities of the near-IR
[Fe ii] and H2 emission lines observed in these objects were compared with our
photoionization models. We also used the [O iii] $\lambda$ 5007 Å/H$\beta$ and
[O i] $\lambda$ 6300 Å/H$\alpha$ line intensity ratios of about 600 000
emission-line galaxies listed in the MPA/JHU Data catalogue of the Sloan
Digital Sky Survey DR7 release (available at http://www.mpa-
garching.mpg.de/SDSS/DR7/).

###  2.2 IFU data

For this study we selected two Sy 1 galaxies, Mrk 1157 and NGC 4151, and two
Sy 2 galaxies (ESO 428-G14 and Mrk 1066). All of them were previously observed
by our group using the IFU spectrographs of the Gemini telescopes. We selected
these objects because both J- and K-band spectroscopic data are available. The
observations of Mrk 1066, Mrk 1157, and NGC 4151 were performed using the
Near-IR Integral field Spectrograph (NIFS; McGregor et al., 2003) on Gemini
North, while ESO 428-G14 was observed with the Gemini Near Infra-Red
Spectrograph (GNIRS; Elias et al., 1998) on Gemini South.

##  3 Photoionization model

To analyse the [Fe ii] and $\rm H_{2}$ excitation mechanisms, we built a grid
of models using the photoionization code Cloudy/08 (Ferland et al., 1998), and
then we compared the line intensity ratios predicted by them with those
observed. The spectral energy distribution (SED) of the ionizing source used
as input for the Cloudy code was a two-component continuum ranging from $\sim
10^{15}$ Hz to $\sim 10^{21}$ Hz. The shape of this SED is similar to the one
observed in typical AGNs for that range. The first is the Big Bump component
peaking at $\rm 1\>Ryd$ with a high-energy and an infrared exponential cutoff
and the second one represents the X-ray source that dominates at high energies
and is characterized by a power law with an index $\alpha_{x}=-1$. Its
normalization was computed to produce the required value of the optical to
X-ray spectral index $\alpha_{ox}$. This index describes the continuum between
2 keV and 2500 Å (Zamorani et al., 1981). We assumed the default value of the
Cloudy code $\alpha_{ox}=-1.4$, because that is about the average of the
observed values, which are between -1.0 and -2.0, for the entire range of
observed luminosities of AGNs (Miller et al., 2011; Zamorani et al., 1981).

The cosmic ray emission was considered in the models as a second ionizing
source. Cosmic rays heat the ionized gas and produce secondary ionizations in
the neutral gas, which mostly increase the intensities of the $\rm H_{2}$
emission lines. We assumed a value of the $\rm H_{2}$ ionization rate of
$10^{-15}\>\rm s^{-1}$, which is about the same rate found by McCall et al.
(2003) for a galactic line of sight. It is worth noting that the value of the
cosmic ray rate must be estimated object by object. For example, Suchkov et
al. (1993) found for M 82 a cosmic ray rate several times larger than the one
in the Milk Way. Gamma ray observations of the starburst NGC 253 by Acero et
al. (2009) indicate a cosmic ray rate three orders of magnitude larger than
that for the Milky Way. Also, molecular data of star forming galaxies, such as
Arp 220, show evidence for extremely high cosmic ray rates yielded by the UV
emission from supernova remnants (Meijerink et al., 2011).

We computed a sequence of models assuming an electron density $N_{\rm
e}=10^{4}\rm\>cm^{-3}$, ionization parameter $U$ in the range $-4.0\leq\log
U\leq-1.0$ defined as $U=Q_{ion}/4\pi R^{2}_{\rm S}nc$, where $Q_{ion}$ is the
number of hydrogen ionizing photons emitted per second by the ionizing source,
$R_{\rm S}$ is the Strömgren radius (in cm), $n$ is the particle density (in
$\rm cm^{-3}$), and $c$ is the speed of light. The chosen range of these
values for $U$ is typical of narrow-line regions of Sy galaxies (e.g. Ferland
& Netzer 1983). The H2 emission lines are very dependent on the electron
density value assumed in the models. For example, when $N_{\rm e}$ varies from
$10^{3}$ to $10^{5}\>\rm cm^{-3}$, the logarithm of the H${}_{2}\,\lambda$
2.1218$\,\mu$m/Br$\gamma$ emission line intensity ratio span about 2.6 dex.
The value $N_{\rm e}=10^{4}\rm cm^{-3}$, assumed in our models, is a mean
value from those considered by Mouri et al. (2000). We considered in our
models three values of 12+log(O/H)= 8.38, 8.69, and 9.00, which correspond to
values of the metallicity 0.5, 1, and 2 times the solar value published by
Allende Prieto et al. (2001). The abundances of other metals in the nebula
were scaled linearly to the solar metal composition through the comparison of
the oxygen abundances, with the exception of the N and Fe abundances. The
nitrogen abundance was taken from the relation log(N/O)=log(0.034+120O/H) of
Vila- Costas & Edmunds (1993). The Fe/O abundance ratio has a large scatter
for a fixed O/H value (Izotov et al., 2006) and its value is uncertain because
the Fe and O abundances in grains are poorly known (Peimbert & Peimbert,
2010). Thus, we varied the Fe/O abundance ratio by about 0.7 dex on each
metallicity.

Table 1: Fe/O and O/H gas phase abundances assumed in the models. Metallicity (Z/Z⊙) | 12+log(O/H) |  | log(Fe/O)
---|---|---|---
2 | 9.0 |  | -1.47(a1) | -1.94(a2) | -2.24(a3)
1 | 8.69 |  | -1.77(b1) | -2.24(b2) | -2.77(b3)
0.5 | 8.38 |  | -2.15(c1) | -2.54(c2) | -2.76(c3)

The presence of internal dust was considered and the grain abundances (van
Hoof et al. 2001) were also linearly scaled with the oxygen abundance. To take
into account the depletion of refractory elements onto dust grains, the
abundances of Mg, Al, Ca, Ni, and Na were reduced by a factor of 10, and Si by
a factor of 2 relative to adopted abundances in each model in accordance with
Garnett et al. (1995). In Table 1, the O/H and Fe/O abundance values of the
gas phase assumed in the models are shown. The model of the $\rm H_{2}$
molecule described by Shaw et al. (2005) and the model of the $\rm Fe^{+}$ ion
described by Verner et al. (1999), which consider 371 energy levels, were
assumed in our computations. The outer radius of the modelled nebula is that
where the temperature falls below 1000 K.

##  4 Diagnostic Diagrams

We used four diagnostic diagrams containing predicted and observed emission
line ratios of the [Fe ii] , $\rm H_{2}$, [O iii], and [O i] which are
described below.

  * •

[Fe ii] $\lambda$ 1.2570$\,\mu$m/Pa$\beta$ vs. H${}_{2}\,\lambda$
2.1218$\,\mu$m/Br$\gamma$ (Fig. 1) — Diagnostic diagram suggested by Larkin et
al. (1998) and Rodríguez-Ardila et al. (2004) to separate galaxies according
to their level of nuclear activity. Recently, Riffel et al. (2010) constructed
this diagram with spatially resolved IFU data of an AGN. Typical values for
the nucleus of Sy galaxies are $0.6\,\lesssim$ [Fe ii] $\lambda$
1.2570$\,\mu$m/Pa$\beta$ $\lesssim\,2.0$ and $0.6\,\lesssim$
H${}_{2}\,\lambda$ 2.1218$\,\mu$m/Br$\gamma$ $\lesssim\,2.0$ (Rodríguez-Ardila
et al., 2005). This [Fe ii]/Pa$\beta$ is very dependent on the Fe/O abundance
while the H2 emission lines are dependent on the ionization parameter.

  * •

H${}_{2}\lambda$ 1.957$\,\mu$m/$\lambda$ 2.121$\,\mu$m and H${}_{2}\lambda$
2.033$\,\mu$m/$\lambda$ 2.223$\,\mu$m vs. H${}_{2}\lambda$
2.247$\,\mu$m/$\lambda$ 2.121$\,\mu$m (Fig. 2) — Mouri (1994) proposed these
diagrams to separate gas emission yielded by shocks from emission caused by
fluorescence. The drawback in using the H${}_{2}\lambda$
1.957$\,\mu$m/$\lambda$ 2.121$\,\mu$m ratio is that the H${}_{2}\lambda$
1.957$\,\mu$m may be affected by telluric bands of $\rm H_{2}O$ and $\rm
CO_{2}$, or blended with the [Si iv] $\lambda$ 1.963$\,\mu$m emission line
(Rodríguez-Ardila et al., 2005).

  * •

[O iii] $\lambda\,5007$Å/H$\beta$ vs. [O i] $\lambda$ 6300Å/H$\alpha$ (Fig. 3)
— This diagram was suggested by Baldwin et al. (1981) to separate objects
according to their primary excitation mechanisms, i.e. (a) photoionization by
stars, (b) photoionization by a power law continuum source or (c) shock
heating. In particular, the [O i] $\lambda$ 6300 Å/H$\alpha$ line ratio is
greatly increased by the presence of shocking gas, even when it has low
velocities (e.g. Allen et al. 2008).

Figure 1: Diagnostic diagram showing the observational data taken from the
literature (see Sect. 2) and results from the grid of photoionization models
(see Sect.3). Solid lines connect curves of iso-$Z$, while dotted lines
connect curves of iso-$U$. The values of $\log U$ and $Z$ are indicated. The
three different lines for each $Z$ correspond to the different assumed values
of the Fe/O as indicated by the labels (see Table 1). Circles, squares, and
star represent Sy1, Sy2, and quasar data, respectively. The typical error bar
(not shown) of the emission line ratios is about 10 %.  Figure 2: As in Fig. 1
for $\rm H_{2}$ emission lines. The arrow indicates the direction in which the
ionization parameter increases. Circles, squares, and star represent Sy1, Sy2,
and quasar data, respectively. The hatched area represents the region occupied
by shock model results from Hollenbach & McKee (1989).

##  5 Results

###  5.1 Integrated spectra

In Figs. 1 and 2 we show the first three diagnostic diagrams described above
containing the results of our grid of photoionization models and the data
sample. Sy 1, Sy 2 and quasar are represented by different symbols. For the
IFU data, the emission line ratios represented in these Figs. were estimated
by integrating the spaxels inside a central aperture of 0.5″$\times$0.5″ for
each galaxy, with exception of ESO428-G14 for which an aperture of
0.75″$\times$0.75″ was considered. These values are presented in Table 2.

In Fig. 1, we can see that almost all the observational ratios are within the
parameter space defined by our grid of photoionization models. A lower
metallicity than those assumed in our models is required to reproduce the data
of the galaxies out of the grid. The observed correlation between [Fe ii]
$\lambda$ 1.2570$\,\mu$m/Pa$\beta$ and H2 $\lambda$ 2.1218$\,\mu$m/Br$\gamma$
is explained by an increase in metallicity and ionization parameter.
Noteworthy that the parameter space defined by the models built using
$Z\,=\,0.5\,Z_{\odot}$ is almost completely contained in the one defined by
the models built using the solar metallicity.

Table 2: Integrated line ratio intensities of IFU data Object | [Fe ii] $\lambda$ 1.2570$\,\mu$m/Pa$\beta$ | H${}_{2}\,\lambda$ 2.1218$\,\mu$m/Br$\gamma$
---|---|---
ESO428-G14 | 0.75 | 1.10
Mrk 1066 | 0.52 | 0.96
Mrk1 157 | 0.73 | 2.24
NGC 4151 | 0.45 | 0.26

In the case of the diagnostic diagrams that only involve H2 line ratios (Fig.
2), the photoionization models are slightly dependent on the assumed
metallicities, covering almost the same parameter space, and strongly
dependent on variations in the ionization parameter. Taking into account the
observational error bars, our models are in good agreement with the observed
H${}_{2}\lambda$ 2.033$\,\mu$m/$\lambda$ 2.223$\,\mu$m and H${}_{2}\lambda$
2.247$\,\mu$m/$\lambda$ 2.121$\,\mu$m ratios (upper panel). On the other hand,
in the lower panel of Fig. 2 can be noticed larger dispersion of the
observational data which is not well reproduce by the models. This dispersion
could be the result of a contamination of the measurements of the
H${}_{2}\,\lambda$ 1.957 $\mu$m emission line intensities due to a blend with
the [Si iv] $\lambda$ 1.963$\,\mu$m line (as explained above). Therefore, the
predicted H${}_{2}\lambda$1.957$\,\mu$m intensities are somewhat lower than
the observed ones. In Fig. 2, we also show the area occupied by the
theoretical intensities of the line ratio H${}_{2}\lambda$
2.247$\,\mu$m/$\lambda$ 2.121$\,\mu$m from shock models performed by
Hollenbach & McKee (1989). These authors computed emission-line spectra of
steady interstellar shocks in molecular gas considering velocities from 30 to
150 km/s and particle densities of $10^{3}-10^{6}\>\rm cm^{-3}$. We can see
that most of the objects of our sample appear to have the X-rays as main
ionizing source while for the remaining ones a composite ionization by X-rays
and shock can be considered.

Fig. 3 shows the [O iii] $\lambda$ 5007 Å/H$\beta$ vs. [O i $\lambda$ 6300
Å/H$\alpha$ diagnostic diagram. In this Fig. we can see that the observational
data of AGNs are well describe by our models. If our models use the lower
values of the ionization parameter ($\log\,U\,<\,-3.5$; these models are not
shown), we can extend the parameter space to include the objects that have
values of the logarithm of the [O iii] $\lambda$ 5007 Å/H$\beta$ ratio lower
than zero. As in the case of the [Fe ii] $\lambda$ 1.2570$\,\mu$m/Pa$\beta$
and H2 $\lambda$ 2.1218$\,\mu$m/Br$\gamma$ diagnostic diagram (Fig. 1), the
parameter space of the models with solar metallicity almost contains that of
the $Z\,=\,0.5\,Z_{\odot}$ models.

Figure 3: [O iii $\lambda$ 5007 Å/H$\beta$ vs. [O i $\lambda$ 6300 Å/H$\alpha$
diagnostic diagram. The yellow line separates objects ionized by massive stars
from those containing active nucleus (Kewley et al., 2001). Blue, green and
red solid lines are as in Fig. 1. Points represent emission-line galaxies
listed in the MPA/JHU Data catalogue of the Sloan Digital Sky Survey DR7
release (see Sect. 2).

###  5.2 IFU data

We plot the [Fe ii] $\lambda$ 1.2570$\,\mu$m/Pa$\beta$ and H2 $\lambda$
2.1218$\,\mu$m/Br$\gamma$ diagnostic diagram for each spaxel of our four
objects with our model results (see upper panels of Figs. 4 and 5). In these
Figs., the spaxel data are separated by their ionization mechanism according
to the place in the diagnostic diagram, with different colours for each
mechanism. The different ionization mechanism zones are delimited in the Figs.
by dashed-lines, following the work of Rodríguez-Ardila et al. (2004). The
spaxels showing typical values of starbursts, Seyferts, and low-ionization
nuclear emission-line regions (LINERs) are represented by green open circles,
black filled circles, and red crosses, respectively. With the same colour
code, we show the spatial position of each spaxel in the IFU field of view
(see lower panels of Figs. 4 and 5). Our models completely represent the
region occupied by Seyfert and LINERs data.

Figure 4: Top panels: [Fe ii] $\lambda\,1.2570\,\mu$m/Pa$\beta$ vs. H2
$\lambda\,2.1212\,\mu$m/Br$\gamma$ line-ratio diagnostic diagram for Mrk 1157
(left) and ESO 428-G14 (right). The dashed lines delimit regions with ratios
typical of Starbursts (green open circles), Seyferts (black filled circles)
and LINERs (red crosses). Blue, green and red solid lines are as in Fig. 1.
Bottom panels: spatial position of each spaxel in the IFU field of view from
the diagnostic diagram.

Figure 5: As in Fig. 4 but for Mrk 1066 (left) and NGC 4151 (right).

##  6 Discussion

The excitation mechanism of the near-IR emission lines of the [Fe ii] and H2
in active galaxies have been the subject of several works. For example, Mouri
et al. (2000) compared results of models considering photoionization by X-rays
and shock heating with observational data of AGNs and starburst galaxies.
These authors built their models considering large ranges in shock velocities,
gas density, metallicity, and different ionizing continua. Mouri and
collaborators showed that the [Fe ii] emission is enhanced when a partially
ionized zone is produced by photoionization by X-rays (described by a power-
law) and shock heating. These two processes can be discriminated by the
electron temperature of the [Fe ii] region: 8000 K in heating by X-rays and
6000 K in shock heating. Comparing the electron temperature of the [Fe ii]
region estimated by Thompson (1995) for NGC 4151 ($8000\><\>T_{\rm
e}\><\>12000$ K) with their models, Mouri et al. (2000) showed that, at least
for this galaxy, it indicates that X-rays are the more important mechanism to
yield the [Fe ii] flux. These authors arrived to the same conclusion using the
[O i] $\lambda$ 6300 Å/H$\alpha$ vs. [Fe ii]$\lambda$1.2570$\,\mu$m]/Pa$\beta$
diagnostic diagram. A similar result was also obtained by Jackson & Beswick
(2007) by analysing J-band spectra of three Sy 2 galaxies.

For our models, we assumed an incident continuum whose shape is given by two
components, a Big Bump and an X-ray power law, varying the Fe/O abundance.
With these models that do not consider shock heating, we are able to explain
the observational data. Nevertheless, we do not exclude some contribution by
shock heating to the [Fe ii] emission. Comparing our models with the SDSS DR7
emission-line galaxies (Fig. 3), we are able to describe the [O i]/H$\alpha$
line ratio observed in AGNs. This diagram cannot be used to discriminate [Fe
ii] excitation mechanism, nevertheless, we must take into account that the [O
i]/H$\alpha$ line ratio is shock sensitive. Hence, although shock contribution
in the ionization of Fe cannot be ruled out, models considering a continuum
described by a Big Bump and an X-ray power law as the ionization source can
also reproduce the [Fe ii] emission lines as well as the behaviour of shock
sensitive emission lines such as [O i] $\lambda$ 6300 Å. Analysing IFU
observations of the Sy galaxy NGC 4151, Turner et al. (2002) found that the
[Fe ii] emission mainly arises in the visible narrow-line region in which the
dominant excitation mechanism is the photoionization by collimated X-ray
emission from the nucleus. Oliva et al. (2001) pointed out that in regions
where shocks are the dominant mechanism the iron-based grains are destroyed
but the phosphorus is not, yielding a larger [Fe ii] $\lambda\,1.2570$
$\mu$m/[P ii] $\lambda\,1.188$ $\mu$m line-ratio intensity than that observed
in the regions dominated by X-ray. In order to verify this, in Fig. 6 we show
a histogram containing this observed line intensity ratio for 17 Seyfert
galaxies, 5 Sy 1 and 12 Sy 2, taken from Jackson & Beswick (2007), Riffel et
al. (2006), and Oliva et al. (2001). It can be seen that the [Fe ii]/[P ii]
for Sy galaxies ranges from 1.5 to 6 (with a mean value of 2.7). The mean
value of this ratio is about 20 for SNRs, which indicate that the emitting gas
has recently passed through a fast shock (Oliva et al., 2001). Therefore,
these results confirm that shocks have little influence on the [Fe ii]
emission.

Figure 6:  Histogram showing the [Fe ii]/[P ii] emission line intensity ratios
of a sample of objects collected from the literature.

Regarding the $\rm H_{2}$, this molecule can be excited via three distinct
mechanisms: (1) UV fluorescence, where photons with $\lambda\,>\,912$ Å are
absorbed by the $\rm H_{2}$ molecule and then re-emitted, resulting in the
population of various vibro-rotational levels, (2) shocks, where high-velocity
gas motions heat and accelerate this molecule; and (3) X-ray illumination,
where hard X-ray photons penetrate deep into molecular clouds, heating large
amounts of molecular gas resulting in the $\rm H_{2}$ emission (see Rodríguez-
Ardila et al., 2004, and references therein). Rodríguez-Ardila and
collaborators used the diagrams shown in Fig. 2 to compare observational data
of 22 objects with models considering a thermal emission, a non-thermal UV
excitation, a thermal UV excitation, and a mixture of thermal and low-density
fluorescence. These authors found that for 4 objects the excitation mechanism
is clearly thermal, while for the remaining objects a mixing with a non-
thermal process cannot be discarded, even though the results point out to a
dominant thermal mechanism.

To analyse the relative weight of the X-ray emission with respect to the other
model components (mainly with fluorescence and UV photons), not only for the
$\rm H_{2}$ emission but also for the [Fe ii], we made models fixing all
parameters with the exception of the $\alpha_{ox}$ value (see Fig. 7), which
is related to the X-ray power law normalization (see Section §3). We assumed
$Z$ = $Z$⊙ and log $U$ = $-$2.5 since the models built using the solar
metallicity and this value of the ionization parameter cover almost all the
parameter space occupied by the observational data (see Fig. 1). Taking into
account the $\alpha_{ox}$ definition (Tananbaum et al., 1979), which fixes the
Big Bump parameters, a decrement of the $\alpha_{ox}$ value implies that the
amount of the X-rays emitted by the source decreases. In Fig. 7 we can see
that our models with $\alpha_{ox}=-1.4$ reproduce well the observational AGN
data. Nevertheless, when we use lower values of this parameter, the ratios
predicted by the models go out of the region typically occupied by the AGNs
(Rodríguez-Ardila et al., 2004). Therefore, our models favour the scenario
suggested by Maloney et al. (1996), where the $\rm H_{2}$ molecule emission is
mainly governed by photons emitted at X-ray wavelengths from the central AGN.
This also can be inferred from the dependence of the $\rm H_{2}$ emission
lines on the ionization parameter $U$.

To verify if shock models can fit the observational data, we compared shock
model results by Hollenbach & McKee (1989) with our sample (see Fig. 2). Only
few observational points are located in the area occupied by these shock
models and, even in these cases, models considering X-rays also describe the
data.

On the other hand, varying in our models the $\rm H_{2}$ ionization rate by
cosmic rays by a factor of 200, we found that the H${}_{2}\,\lambda$
2.1218$\,\mu$m/Br$\gamma$ line ratio only increases by about 0.15 dex, which
shows that the additional ionization by cosmic rays has little influence on
the $\rm H_{2}$ emission lines.

Figure 7: Model results using solar metallicity, Ne = $10^{4}$ cm-3, log U =
-2.5 and varying only the $\alpha_{ox}$ parameter to see the influence of the
X-rays on the [Fe ii] $\lambda\,1.2570\,\mu$m/Pa$\beta$ and H2
$\lambda\,2.1218\,\mu$m/Br$\gamma$ emission line ratios. To delimit the region
occupied by the AGNs we follow Rodríguez-Ardila et al. (2004). Figure 8:
Spectral energy distribution at the Schwarzschild radio (10-5 pc) of the Sy 1
galaxy Mrk 1066 used as the photoionization source for some models of this
galaxy. We assumed a galaxy distance of 50 Mpc (Mould et al., 2000). The
photometric data were taken from Dressel & Condon (1978), Moshir et al.
(1990), de Vaucouleurs et al. (1991), Becker et al. (1991), Douglas et al.
(1996), Condon et al. (2002), Skrutskie et al. (2003), Braatz et al. (2004),
Guainazzi et al. (2005), Muñoz-Marín et al. (2007), and Cardamone et al.
(2007).  Figure 9: Comparison between the grid model results shown in Fig. 1
(solid lines) and the models built considering the semi-empirical SED of Mrk
1066 shown in Fig. 8 (dashed lines).

A simple scenario where both [Fe ii] and $\rm H_{2}$ emissions are mainly due
to the X-ray continuum coming from the active nucleus has also been proposed
by other authors. For example, Blietz et al. (1994) and Knop et al. (2001)
showed that X-rays from the nucleus can heat the gas located in the narrow
line region driving the [Fe ii] and H 2 emission. Because 98 % of the iron is
tied up in dust grains, this process must free the iron through dust
destruction and yet not destroy the $\rm H_{2}$ molecules (Rodríguez-Ardila et
al., 2004). These authors computed the emergent [Fe ii] $\lambda$ 1.2570
$\,\mu$m and H${}_{2}\,\lambda$ 2.1218$\,\mu$m flux using the X-ray models by
Maloney et al. (1996) and compared their predictions with observational data
of seven objects. They found that X-ray heating can only explain a fraction of
the [Fe ii] and H2 emission, and they stated that the discrepancy found can be
alleviated if the emitting gas is located closer than the distance adopted in
their models. The X-ray data, provided by the XMM-Newton and Chandra space
telescopes, and their detailed analysis (see e.g. Piconcelli et al., 2005;
Longinotti et al., 2007; Bianchi et al., 2009; Krongold et al., 2009; Cardaci
et al., 2011; Corral et al., 2011, and references therein), provide
information about the continuum shape and the particular spectral features of
the AGNs in this wavelength range. For Mrk 1066, we compared the results
obtained using this simple scenario that only involves a continuum modelled by
a Big Bump and an X-ray power law with those obtained using its intrinsic SED.
We built the observational SED taking the photometric data from the NASA/IPAC
Extragalactic Database (NED), following Cardaci et al. (2009). To enhance the
number of points of the SED as needed by Cloudy, we performed a linear
interpolation among the semi-empirical points (see Fig. 8). We built a new
grid of photoionization models under the same assumptions of abundances,
ionization parameters and density, but only for one value of the Fe/O ratio
for each metallicity. In Fig. 9 the predictions of our models using the SED of
Mrk 1066 and the model results presented in Fig. 1 assuming the same Fe/O
abundance as the Mrk 1066 models are shown. The model results derived using
the two different ionizing sources are mostly in agreement.

The semi-empirical SED of Mrk 1066 includes not only the range covered by the
Cloudy model but also the radio and IR wavelengths. Hence, the agreement
between solid and dashed lines in Fig. 9 only indicates that the assumed
multicomponent model is a good representation of the AGN continuum when
studying the [Fe ii] and H2 emission.

Recent resolved integral field spectroscopy of the central region of active
galaxies shows that the ionized (in particular the [Fe ii] emitting gas) and
the molecular (traced by the H2 emission) gas have distinct flux distributions
and kinematics. The molecular component is more restricted to the plane of the
galaxies and the ionized one extends to high latitudes above it, which is in
most cases co-spatially with the radio jet (e.g., Riffel et al., 2006, 2008,
2009, 2010; Riffel & Storchi-Bergmann, 2011a, b; Storchi-Bergmann et al.,
2009; Storchi-Bergmann et al., 2010). Usually the [Fe ii] has an enhancement
in flux and velocity dispersion in regions surrounding the radio structure,
suggesting that the radio jet plays an important role in the [Fe ii] emission.
Our models are able to reproduce the [Fe ii] emission of active galaxies
without considering shock excitation by the radio jet. Thus, the enhancement
in the [Fe ii] flux in the vicinity of radio structures can be interpreted as
being due to an enhancement in the gas density, caused by the interaction of
the radio jet with the emitting gas, and mainly excited by X-rays from the
central engine.

The main exciting mechanism of infrared emission lines of ESO 428-G14, Mrk
1157, Mrk 1066 and NGC 4151 was discussed by Riffel et al. (2006, 2008, 2010)
and Storchi-Bergmann et al. (2009), respectively. However, these authors did
not reach conclusive results. For example, Riffel et al. (2006) suggested that
the [Fe ii] excitation in ESO 428-G14 is mainly due to shocks. Nevertheless,
the detailed analysis performed in the present work confronting our models
with the IFU data shows that X-rays are a more reliable dominant excitation
mechanism even in the case of ESO 428-G14.

##  7 Conclusions

In this work we show that a photoionization model grid built by adopting a
continuum source characterized by two components, one accounting for the Big
Bump component peaking at $\rm 1\>Ryd$ and the other describing the X-rays
emission, is able to reproduce the [Fe ii] and H2 infrared emission lines of a
sample of AGNs. Testing the influence of the X-rays on the intensity of these
emission lines, we found that a decrement in the X-ray content of the
continuum source translates into a weakening of these lines, and the models
are no longer compatible with the observations. This implies that the heating
by the X-ray emission from the active nuclei can be considered as the most
important mechanism of excitation for the IR emission lines of these elements.

## Acknowledgments

We are grateful to the referee, Neal Jackson, for a thorough reading of the
manuscript and for suggestions that greatly improved its clarity.

Based on observations obtained at the Gemini Observatory, which is operated by
the Association of Universities for Research in Astronomy, Inc., under a
cooperative agreement with the NSF on behalf of the Gemini partnership: the
National Science Foundation (United States), the Science and Technology
Facilities Council (United Kingdom), the National Research Council (Canada),
CONICYT (Chile), the Australian Research Council (Australia), Ministério da
Ciência e Tecnologia (Brazil) and south-eastCYT (Argentina). This research has
made use of the NASA/IPAC Extragalactic Database (NED) which is operated by
the Jet Propulsion Laboratory, California Institute of Technology, under
contract with the National Aeronautics and Space Administration. OLD and ACK
are grateful to the FAPESP for support under grant 2009/14787-7 and
2010/01490-3. MC and GH are grateful to the Spanish Ministerio de Ciencia e
Innovación for support under grant AYA2010-21887-C04-03, and the Comunidad de
Madrid under grant S2009/ESP-1496 (ASTROMADRID). EPM is grateful to the
Spanish Ministerio de Ciencia e Innovación for support under grant
AYA2010-21887-C04-02, and the Junta de Andalucía under grant TIC114.

## References

  * Acero et al. (2009) Acero F., Aharonian F., Akhperjanian A. G., Anton G., Barres de Almeida U., Bazer-Bachi A. R., Becherini Y., Behera B., Bernlöhr K., Bochow A., and 160 coauthors, 2009, Science, 326, 1080
  * Allen et al. (2008) Allen M. G., Groves B. A., Dopita M. A., Sutherland R. S., Kewley L. J., 2008, Astrophys. J., Suppl. Ser., 178, 20
  * Allende Prieto et al. (2001) Allende Prieto C., Lambert D. L., Asplund M., 2001, Astrophys. J. Letters, 556, L63
  * Baldwin et al. (1981) Baldwin J. A., Phillips M. M., Terlevich R., 1981, Publ. Astron. Soc. Pac., 93, 5
  * Becker et al. (1991) Becker R. H., White R. L., Edwards A. L., 1991, Astrophys. J., Suppl. Ser., 75, 1
  * Bianchi et al. (2009) Bianchi S., Guainazzi M., Matt G., Fonseca Bonilla N., Ponti G., 2009, Astron. Astrophys., 495, 421
  * Black & van Dishoeck (1987) Black J. H., van Dishoeck E. F., 1987, Astrophys. J., 322, 412
  * Blietz et al. (1994) Blietz M., Cameron M., Drapatz S., Genzel R., Krabbe A., van der Werf P., Sternberg A., Ward M., 1994, Astrophys. J., 421, 92
  * Braatz et al. (2004) Braatz J. A., Henkel C., Greenhill L. J., Moran J. M., Wilson A. S., 2004, Astrophys. J. Letters, 617, L29
  * Cardaci et al. (2011) Cardaci M. V., Santos-Lleó M., Hägele G. F., Krongold Y., Díaz A. I., Rodríguez-Pascual P., 2011, Astron. Astrophys., 530, A125+
  * Cardaci et al. (2009) Cardaci M. V., Santos-Lleó M., Krongold Y., Hägele G. F., Díaz A. I., Rodríguez-Pascual P., 2009, Astron. Astrophys., 505, 541
  * Cardamone et al. (2007) Cardamone C. N., Moran E. C., Kay L. E., 2007, Astron. J., 134, 1263
  * Condon et al. (2002) Condon J. J., Cotton W. D., Broderick J. J., 2002, Astron. J., 124, 675
  * Corral et al. (2011) Corral A., Della Ceca R., Caccianiga A., Severgnini P., Brunner H., Carrera F. J., Page M. J., Schwope A. D., 2011, Astron. Astrophys., 530, A42+
  * Davies et al. (2007) Davies R. I., Müller Sánchez F., Genzel R., Tacconi L. J., Hicks E. K. S., Friedrich S., Sternberg A., 2007, Astrophys. J., 671, 1388
  * de Vaucouleurs et al. (1991) de Vaucouleurs G., de Vaucouleurs A., Corwin Jr. H. G., Buta R. J., Paturel G., Fouque P., 1991, Third Reference Catalogue of Bright Galaxies. Volume 1-3, XII, 2069 pp. 7 figs.. Springer-Verlag Berlin Heidelberg New York
  * Douglas et al. (1996) Douglas J. N., Bash F. N., Bozyan F. A., Torrence G. W., Wolfe C., 1996, Astron. J., 111, 1945
  * Dressel & Condon (1978) Dressel L. L., Condon J. J., 1978, Astrophys. J., Suppl. Ser., 36, 53
  * Elias et al. (1998) Elias J. H., Vukobratovich D., Andrew J. R., Cho M. K., Cuberly R. W., Don K., Gerzoff A., Harmer C. F., Harris D., Heynssens J. B., Hicks J., Kovacs A., Li C., Liang M., Moon I. K., Pearson E. T., Plum G., Roddier N. A., Tvedt J., Wolff R. J., Wong W.-Y., 1998, in A. M. Fowler ed., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 3354 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Design of the Gemini near-infrared spectrometer. pp 555–565
  * Ferland et al. (1998) Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, Publ. Astron. Soc. Pac., 110, 761
  * Ferland & Netzer (1983) Ferland G. J., Netzer H., 1983, Astrophys. J., 264, 105
  * Ferruit et al. (2000) Ferruit P., Wilson A. S., Mulchaey J., 2000, Astrophys. J., Suppl. Ser., 128, 139
  * Forbes & Ward (1993) Forbes D. A., Ward M. J., 1993, Astrophys. J., 416, 150
  * Garnett et al. (1995) Garnett D. R., Dufour R. J., Peimbert M., Torres-Peimbert S., Shields G. A., Skillman E. D., Terlevich E., Terlevich R. J., 1995, Astrophys. J. Letters, 449, L77+
  * Guainazzi et al. (2005) Guainazzi M., Matt G., Perola G. C., 2005, Astron. Astrophys., 444, 119
  * Hicks et al. (2009) Hicks E. K. S., Davies R. I., Malkan M. A., Genzel R., Tacconi L. J., Müller Sánchez F., Sternberg A., 2009, Astrophys. J., 696, 448
  * Hollenbach & McKee (1989) Hollenbach D., McKee C. F., 1989, Astrophys. J., 342, 306
  * Izotov et al. (2006) Izotov Y. I., Stasińska G., Meynet G., Guseva N. G., Thuan T. X., 2006, Astron. Astrophys., 448, 955
  * Jackson & Beswick (2007) Jackson N., Beswick R. J., 2007, Mon. Not. R. Astron. Soc., 376, 719
  * Kewley et al. (2001) Kewley L. J., Dopita M. A., Sutherland R. S., Heisler C. A., Trevena J., 2001, Astrophys. J., 556, 121
  * Knop et al. (2001) Knop R. A., Armus L., Matthews K., Murphy T. W., Soifer B. T., 2001, Astron. J., 122, 764
  * Krongold et al. (2009) Krongold Y., Jiménez-Bailón E., Santos-Lleo M., Nicastro F., Elvis M., Brickhouse N., Andrade-Velazquez M., Binette L., Mathur S., 2009, Astrophys. J., 690, 773
  * Larkin et al. (1998) Larkin J. E., Armus L., Knop R. A., Soifer B. T., Matthews K., 1998, Astrophys. J., Suppl. Ser., 114, 59
  * Longinotti et al. (2007) Longinotti A. L., Sim S. A., Nandra K., Cappi M., 2007, Mon. Not. R. Astron. Soc., 374, 237
  * Maloney et al. (1996) Maloney P. R., Hollenbach D. J., Tielens A. G. G. M., 1996, Astrophys. J., 466, 561
  * McCall et al. (2003) McCall B. J., Huneycutt A. J., Saykally R. J., Geballe T. R., Djuric N., Dunn G. H., Semaniak J., Novotny O., Al-Khalili A., Ehlerding A., Hellberg F., Kalhori S., Neau A., Thomas R., Österdahl F., Larsson M., 2003, Nature, 422, 500
  * McGregor et al. (2003) McGregor P. J., Hart J., Conroy P. G., Pfitzner M. L., Bloxham G. J., Jones D. J., Downing M. D., Dawson M., Young P., Jarnyk M., Van Harmelen J., 2003, in M. Iye & A. F. M. Moorwood ed., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 4841 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Gemini near-infrared integral field spectrograph (NIFS). pp 1581–1591
  * Meijerink et al. (2011) Meijerink R., Spaans M., Loenen A. F., van der Werf P. P., 2011, Astron. Astrophys., 525, A119+
  * Miller et al. (2011) Miller B. P., Brandt W. N., Schneider D. P., Gibson R. R., Steffen A. T., Wu J., 2011, Astrophys. J., 726, 20
  * Moshir et al. (1990) Moshir M., Kopan G., Conrow T., McCallon H., Hacking P., Gregorich D., Rohrbach G., Melnyk M., Rice W., Fullmer L., White J., Chester T., 1990, in Bulletin of the American Astronomical Society Vol. 22 of Bulletin of the American Astronomical Society, The IRAS Faint Source Catalog, Version 2. pp 1325–+
  * Mould et al. (2000) Mould J. R., Huchra J. P., Freedman W. L., Kennicutt Jr. R. C., Ferrarese L., Ford H. C., Gibson B. K., Graham J. A., Hughes S. M. G., Illingworth G. D., Kelson D. D., Macri L. M., Madore B. F., Sakai S., Sebo K. M., Silbermann N. A., Stetson P. B., 2000, Astrophys. J., 529, 786
  * Mouri (1994) Mouri H., 1994, Astrophys. J., 427, 777
  * Mouri et al. (1993) Mouri H., Kawara K., Taniguchi Y., 1993, Astrophys. J., 406, 52
  * Mouri et al. (2000) Mouri H., Kawara K., Taniguchi Y., 2000, Astrophys. J., 528, 186
  * Mouri et al. (1990) Mouri H., Nishida M., Taniguchi Y., Kawara K., 1990, Astrophys. J., 360, 55
  * Muñoz-Marín et al. (2007) Muñoz-Marín V. M., González Delgado R. M., Schmitt H. R., Cid Fernandes R., Pérez E., Storchi-Bergmann T., Heckman T., Leitherer C., 2007, Astron. J., 134, 648
  * Mulchaey et al. (1996) Mulchaey J. S., Wilson A. S., Tsvetanov Z., 1996, Astrophys. J., 467, 197
  * Müller Sánchez et al. (2009) Müller Sánchez F., Davies R. I., Genzel R., Tacconi L. J., Eisenhauer F., Hicks E. K. S., Friedrich S., Sternberg A., 2009, Astrophys. J., 691, 749
  * Oliva et al. (2001) Oliva E., Marconi A., Maiolino R., Testi L., Mannucci F., Ghinassi F., Licandro J., Origlia L., Baffa C., Checcucci A., and 12 coauthors, 2001, Astron. Astrophys., 369, L5
  * Peimbert & Peimbert (2010) Peimbert A., Peimbert M., 2010, Astrophys. J., 724, 791
  * Piconcelli et al. (2005) Piconcelli E., Jimenez-Bailón E., Guainazzi M., Schartel N., Rodríguez-Pascual P. M., Santos-Lleó M., 2005, Astron. Astrophys., 432, 15
  * Quillen et al. (1999) Quillen A. C., Alonso-Herrero A., Rieke M. J., Rieke G. H., Ruiz M., Kulkarni V., 1999, Astrophys. J., 527, 696
  * Ramos Almeida et al. (2009) Ramos Almeida C., Pérez García A. M., Acosta-Pulido J. A., 2009, Astrophys. J., 694, 1379
  * Reunanen et al. (2002) Reunanen J., Kotilainen J. K., Prieto M. A., 2002, Mon. Not. R. Astron. Soc., 331, 154
  * Riffel et al. (2006) Riffel R., Rodríguez-Ardila A., Pastoriza M. G., 2006, Astron. Astrophys., 457, 61
  * Riffel & Storchi-Bergmann (2011a) Riffel R. A., Storchi-Bergmann T., 2011a, Mon. Not. R. Astron. Soc., 411, 469
  * Riffel & Storchi-Bergmann (2011b) Riffel R. A., Storchi-Bergmann T., 2011b, Mon. Not. R. Astron. Soc., pp 1320–+
  * Riffel et al. (2009) Riffel R. A., Storchi-Bergmann T., Dors O. L., Winge C., 2009, Mon. Not. R. Astron. Soc., 393, 783
  * Riffel et al. (2010) Riffel R. A., Storchi-Bergmann T., Nagar N. M., 2010, Mon. Not. R. Astron. Soc., 404, 166
  * Riffel et al. (2006) Riffel R. A., Storchi-Bergmann T., Winge C., Barbosa F. K. B., 2006, Mon. Not. R. Astron. Soc., 373, 2
  * Riffel et al. (2008) Riffel R. A., Storchi-Bergmann T., Winge C., McGregor P. J., Beck T., Schmitt H., 2008, Mon. Not. R. Astron. Soc., 385, 1129
  * Rodríguez-Ardila et al. (2004) Rodríguez-Ardila A., Pastoriza M. G., Viegas S., Sigut T. A. A., Pradhan A. K., 2004, Astron. Astrophys., 425, 457
  * Rodríguez-Ardila et al. (2005) Rodríguez-Ardila A., Riffel R., Pastoriza M. G., 2005, Mon. Not. R. Astron. Soc., 364, 1041
  * Shaw et al. (2005) Shaw G., Ferland G. J., Abel N. P., Stancil P. C., van Hoof P. A. M., 2005, Astrophys. J., 624, 794
  * Simpson et al. (1996) Simpson C., Forbes D. A., Baker A. C., Ward M. J., 1996, Mon. Not. R. Astron. Soc., 283, 777
  * Skrutskie et al. (2003) Skrutskie M. F., Cutri R. M., Stiening R., Weinberg M. D., Schneider S., Carpenter J. M., Beichman C., Capps R., Chester T., Elias J., and 21 coauthors, 2003, VizieR Online Data Catalog, 7233, 0
  * Storchi-Bergmann et al. (2010) Storchi-Bergmann T., Lopes R. D. S., McGregor P. J., Riffel R. A., Beck T., Martini P., 2010, Mon. Not. R. Astron. Soc., 402, 819
  * Storchi-Bergmann et al. (2009) Storchi-Bergmann T., McGregor P. J., Riffel R. A., Simões Lopes R., Beck T., Dopita M., 2009, Mon. Not. R. Astron. Soc., 394, 1148
  * Suchkov et al. (1993) Suchkov A., Allen R. J., Heckman T. M., 1993, Astrophys. J., 413, 542
  * Tananbaum et al. (1979) Tananbaum H., Avni Y., Branduardi G., Elvis M., Fabbiano G., Feigelson E., Giacconi R., Henry J. P., Pye J. P., Soltan A., Zamorani G., 1979, Astrophys. J. Letters, 234, L9
  * Thompson (1995) Thompson Jr. R. J., 1995, Astrophys. J., 454, 660
  * Turner et al. (2002) Turner J. E. H., Allington-Smith J., Chapman S., Content R., Done C., Haynes R., Lee D., Morris S., 2002, Mon. Not. R. Astron. Soc., 331, 284
  * Veilleux et al. (1997) Veilleux S., Goodrich R. W., Hill G. J., 1997, Astrophys. J., 477, 631
  * Verner et al. (1999) Verner E. M., Verner D. A., Korista K. T., Ferguson J. W., Hamann F., Ferland G. J., 1999, Astrophys. J., Suppl. Ser., 120, 101
  * Zamorani et al. (1981) Zamorani G., Henry J. P., Maccacaro T., Tananbaum H., Soltan A., Avni Y., Liebert J., Stocke J., Strittmatter P. A., Weymann R. J., Smith M. G., Condon J. J., 1981, Astrophys. J., 245, 357