A&A 392, 827-839 (2002)
DOI: 10.1051/0004-6361:20020951

Radiative transfer with finite elements

II. Ly$\alpha$ line transfer in moving media

E. Meinköhn^1,2 - S. Richling¹

1 - Institut für Theoretische Astrophysik, Universität Heidelberg, Tiergartenstr. 15, 69121 Heidelberg
2 - Institut für Angewandte Mathematik, Universität Heidelberg, Im Neuenheimer Feld 294, 69120 Heidelberg

Received 8 May 2002 / Accepted 24 June 2002

Abstract
A finite element method for solving the resonance line transfer problem in moving media is presented. The algorithm works in three spatial dimensions on unstructured grids which are adaptively refined by means of an a posteriori error indicator. Frequency discretization is implemented via a first-order Euler scheme. We discuss the resulting matrix structure for coherent isotropic scattering and complete redistribution. The solution is performed using an iterative procedure, where monochromatic radiative transfer problems are successively solved. The present implementation is applicable for arbitrary model configurations with an optical depth up to 10^3-4. Results of Ly$\alpha$ line transfer calculations for a spherically symmetric model, a disk-like configuration, and a halo containing three source regions are discussed. We find the characteristic double-peaked Ly$\alpha$ line profile for all models with an optical depth $\ga$1. In general, the blue peak of the profile is enhanced for models with infall motion and the red peak for models with outflow motion. Both velocity fields produce a triangular shape in the two-dimensional Ly$\alpha$ spectra, whereas rotation creates a shear pattern. Frequency-resolved Ly$\alpha$ images may help to find the number and position of multiple Ly$\alpha$ sources located in a single halo. A qualitative comparison with observations of extended Ly$\alpha$ halos associated with high redshift galaxies shows that even models with lower hydrogen column densities than required from profile fitting yield results which reproduce many features in the observed line profiles and two-dimensional spectra.

Key words: line: formation - radiative transfer - scattering - methods: numerical - galaxies: high-redshift

1 Introduction

Hydrogen Ly$\alpha$ as a prominent emission line of high redshift galaxies is important for the understanding of galaxy formation and evolution in the early universe. Recent narrow-band imaging and spectroscopic surveys used the Ly$\alpha$ line to identify galaxies at very high redshift (e.g. Hu et al. 1998; Kudritzki et al. 2000; Rhoads et al. 2000; Fynbo et al. 2000, 2001).

But beside from being a good redshift indicator, Ly$\alpha$ emission also bears information on the distribution and kinematics of the interstellar gas as well as the nature of the energy source. Ly$\alpha$ observations of high redshift radio galaxies (e.g. Hippelein & Meisenheimer 1993; van Ojik et al. 1996, 1997; Villar-Martín et al. 1999) reveal extended Ly$\alpha$ halos with sizes >100 kpc which are aligned with the radio jet. The line profiles are often double-peaked and the two-dimensional spectra point to complex kinematics involving velocities >1000 km $\;\mbox{s}^{-1}$.

The interpretation of Ly$\alpha$ observations is difficult, because high-redshift radio galaxies tend to be in the center of proto clusters, where the radio jet interacts with a clumpy environment influenced by merging processes (e.g. Bicknell et al. 2000; Kurk et al. 2001). Actually, the three-dimensional structure of the objects and the fact that Ly$\alpha$ is a resonance line require detailed radiative transfer modeling.

The transfer of resonance line photons is profoundly determined by scattering in space and frequency. Analytical (Neufeld 1990) as well as early numerical methods (Adams 1972; Hummer & Kunasz 1980) were restricted to one-dimensional, static media. Only recently, codes based on the Monte Carlo method were developed which are capable to investigate the more general case of a multi-dimensional medium (Ahn et al. 2001, 2002).

In this paper, we introduce a finite element method for the solution of the resonance line transfer problem in three-dimensional, moving media and present the results of some simple, slightly optically thick model configurations. The basic, monochromatic code was originally developed by Kanschat (1996) and is described in Richling et al. (2001, hereafter Paper I). The three-dimensional method is particularly useful for scattering dominated problems in inhomogeneous media. Steep gradients are resolved by means of an adaptively refined grid. Here, we only specify the extension from the monochromatic to the polychromatic problem including the implementation of the Doppler-effect and complete redistribution.

In Sect. 2, we review the equations for resonance line transfer in moving media. In Sect. 3, we describe some details regarding frequency discretization and the form of the resulting matrices for coherent scattering and complete redistribution and give a short outline of the complete finite element algorithm. Then, the results of a spherically symmetric model (Sect. 4), a disk-like configuration (Sect. 5) and a model with three separate source regions (Sect. 6) are presented. A summary is given in Sect. 7.

2 Line transfer in moving media

We calculate the frequency-dependent radiation field in moving media by solving the non-relativistic radiative transfer equation in the comoving frame for a three-dimensional domain $\Omega$ which in operator form simply reads

$$\bigl({\cal T}+{\cal D}+{\cal S}+\chi(\vec{x},\nu)\bigr){\cal I}(\vec{x},\vec{n},\nu)=f(\vec{x},\nu).$$ (1)

The transfer operator ${\cal T}$, the "Doppler'' operator ${\cal D}$, and the scattering operator ${\cal S}$ are defined by

$${{\cal T}{\cal I}(\vec{x},\vec{n},\nu)=\vec{n}\cdot\nabla_x {\cal I}(\vec{x},\vec{n},\nu),}$$

$${{\cal D}{\cal I}(\vec{x},\vec{n},\nu)= w(\vec{x},\vec{n})~\nu\frac{\partial}{\partial\nu}{\cal I}(\vec{x},\vec{n},\nu),}$$

$${{\cal S}{\cal I}(\vec{x},\vec{n},\nu)= -\frac{\sigma(\vec{x})}{4... ...\cal I}(\vec{x},\hat{\vec{n}},\hat{\nu})~{\rm d}\hat{\omega}~{\rm d}\hat{\nu}.}$$

Considering three dimensions in space, the relativistic invariant specific intensity ${\cal I}$ is six-dimensional and depends on the space variable $\vec{x}$, the direction $\vec{n}$ (pointing to the solid angle ${\rm d}\omega$ of the unit sphere S²), and the frequency $\nu$.

The extinction coefficient $\chi(\vec{x},\nu)=\kappa(\vec{x},\nu)+\sigma(\vec{x},\nu)$ is the sum of the absorption coefficient $\kappa(\vec{x},\nu)=\kappa(\vec{x})~\phi(\nu)$ and the scattering coefficient $\sigma(\vec{x},\nu)=\sigma(\vec{x})~\phi(\nu)$. The frequency-dependence is given by a normalized profile function $\phi$. Usually, we adopt a Doppler profile

$$\phi(\nu)=\frac{1}{\sqrt{\pi}~\Delta\nu_{\rm D}} \exp{\left[-\left(\frac{\nu-\nu_0}{\Delta\nu_{\rm D}}\right)^2\right]},$$ (2)

where $\nu_0$ is the frequency of the line center. The Doppler width  $\Delta\nu_{\rm D}$ and the Doppler velocity $v_{\rm D}$ are determined by a thermal velocity $v_{\rm therm}$ as well as a turbulent velocity $v_{\rm turb}$

$$\Delta\nu_{\rm D}=\frac{\nu_0}{c} v_{\rm D}= \frac{\nu_0}{c} \sqrt{v_{\rm therm}^2+v_{\rm turb}^2}.$$ (3)

c is the speed of light. For situations, where $v_{\rm turb}\gg v_{\rm therm}$, the Doppler core is very broad and would dominate the Lorentzian wings of a Voigt profile. Then, the Doppler profile is a reasonable description of a line profile. Throughout this paper we use $v_{\rm D}=10^{-3}c\sim300\;\mbox{km}~\mbox{s}^{-1}$ (see also Table 1).

For the source term

$$f(\vec{x},\nu)=\kappa(\vec{x},\nu)B(T(\vec{x}),\nu)+\epsilon(\vec{x},\nu),$$ (4)

we can consider thermal radiation and non-thermal radiation. In the case of thermal radiation, f is calculated from a temperature distribution $T(\vec{x})$, where $B(T,\nu)$ is the Planck function.

The "Doppler'' operator ${\cal D}$ is responsible for the Doppler shift of the photons. A derivation of the operator for non-relativistic velocities (v/c < 0.1) can be found in Wehrse et al. (2000). In contrast to the full relativistic transfer equation (cf. Mihalas & Weibel-Mihalas 1984), we neglect any terms responsible for aberration or advection effects. The function

$$w(\vec{x},\vec{n})=-\vec{n}\cdot\nabla_x \left(\vec{n}\cdot \frac{\vec{v}(\vec{x})}{c} \right)$$ (5)

is the gradient of the velocity field $\vec{v}(\vec{x})$ in direction $\vec{n}$. Note that the sign of w may change depending on the complexity of the velocity field $\vec{v}$.

The scattering operator ${\cal S}$ depends on the general redistribution function $R(\hat{\vec{n}},\hat{\nu};\vec{n},\nu)$ giving the probability that a photon $(\hat{\vec{n}},\hat{\nu})$ is absorbed and a photon $(\vec{n},\nu)$ is emitted. In the following, we assume isotropic scattering

$${\cal S} {\cal I} = -\frac{\sigma(\vec{x})}{4\pi} \int\limit... ...\hat{\vec{n}},\hat{\nu})~{\rm d}\hat{\omega}~{\rm d}\hat{\nu},$$ (6)

where $R(\hat{\nu},\nu)$ is the angle-averaged redistribution function

$$R(\hat{\nu},\nu)=\frac{1}{(4\pi)^2}\int_{S^2}\int_{S^2} R(\h... ...{n}},\hat{\nu};\vec{n},\nu)~{\rm d}\hat{\omega}~{\rm d}\omega.$$ (7)

The function defined by (7) is normalized such that

$$\int\limits_0^\infty \int\limits_0^\infty R(\hat{\nu},\nu)~{\rm d}\hat{\nu}~{\rm d}\nu = 1.$$ (8)

For the calculations in Sect. 3, we consider two limiting cases: strict coherence and complete redistribution in the comoving frame. In the former case, we have

$$R(\hat{\nu},\nu) = \phi(\hat{\nu}) \delta (\nu-\hat{\nu})$$ (9)

and in the latter

$$R(\hat{\nu},\nu) = \phi(\hat{\nu}) \phi(\nu).$$ (10)

Thus, for coherent isotropic scattering, the scattering term simplifies to

$${\cal S}^{\rm coh}{\cal I}(\vec{x},\vec{n},\nu)=-\frac{\sigma... ...{S^2}{\cal I}(\vec{x},\hat{\vec{n}},\nu)~{\rm d}\hat{\omega}.$$ (11)

In the case of complete redistribution, the photons are scattered isotropically in angle, but are randomly redistributed over the line profile. Then, the scattering term reads

$${\cal S}^{\rm crd}{\cal I}(\vec{x},\vec{n},\nu) = -\frac{\s... ...at{\vec{n}},\hat{\nu})~{\rm d}\hat{\omega}~ {\rm d}\hat{\nu}.$$ (12)

3 The finite element method

3.1 Boundary conditions

For the modeling of prominent resonance lines, in particular Ly$\alpha$, we restrict the frequency discretization to a finite interval $\Lambda:=[\nu_0,\nu_{N+1}]$, where $\nu_0$ and $\nu_{N+1}$are located far out of the line center in the continuum. The function $w(\vec{x},\vec{n})$ from Eq. (5) defines the Doppler shift of the spectrum towards smaller (w<0) or larger (w>0) frequencies in the comoving frame. Therefore, boundary conditions of the form

$${\cal I}(\vec{x},\vec{n},\nu)={\cal I}_{\rm cont}(\vec{x},\vec{n},\nu)$$

are necessary on the upper and lower frequency interval boundary of the whole computational domain

$$\Sigma = \Omega \times {\cal S}^2 \times \{\nu_0,\nu_{N+1}\}.$$

Furthermore, to be able to solve Eq. (1), boundary conditions of the form

$${\cal I}(\vec{x},\vec{n},\nu)={\cal I}_{\rm in}(\vec{x},\vec{n},\nu)$$

must be imposed on the "inflow boundary''

$$\Gamma^- \times \Lambda = \bigl\{ (\vec{x},\vec{n},\nu) \in \Gamma \vert \mbox{~} \vec{n}_\Gamma \cdot \vec{n} < 0 \bigr\},$$

where $\vec{n}_\Gamma$ is the unit vector perpendicular to the boundary surface $\Gamma$ of the spatial domain $\Omega$. The sign of the product $\vec{n}_\Gamma \cdot \vec{n}$ describes the "flow direction'' of the photons across the boundary. If we neglect any continuum emission ( ${\cal I}_{\rm cont}=0$) and assume that there are no light sources outside the modeled domain as in the case of a non-illuminated atmosphere ( ${\cal I}_{\rm in}=0$), the two boundary conditions for the solution of the transfer Eq. (1) in moving media are

$${\cal I}(\vec{x},\vec{n},\nu) 0 \mbox{\qquad on }\; \Sigma ,$$ = (13)

$${\cal I}(\vec{x},\vec{n},\nu) 0 \mbox{\qquad on }\; \Gamma^- \times \Lambda .$$ = (14)

3.2 Discretization for coherents scattering

In order to solve the radiative transfer Eq. (1), we first carry out the frequency discretization including coherent scattering. With the abbreviation

$${\cal A}^{\rm coh}={\cal T}+{\cal S}^{\rm coh}+\chi(\vec{x},\nu)$$ (15)

Eq. (1) can be written as

$${\cal A}^{\rm coh}{\cal I}(\vec{x},\vec{n},\nu)+w(\vec{x},\v... ...l}{\partial \nu}{\cal I}(\vec{x},\vec{n},\nu) =f(\vec{x},\nu).$$ (16)

In Paper I, we described a Galerkin discretization for the monochromatic radiative transfer equation. Considering memory and CPU requirements this approach is by far too "expensive'' for the frequency-dependent transfer problem. Instead, we use an implicit Euler method for N equidistantly distributed frequency points $\nu_i \in \bigl\{ \nu_1,\nu_2,...,\nu_N \bigr\} \subset \Lambda$. Since the function $w(\vec{x},\vec{n})$ may change its sign, this simple difference scheme for the Doppler term reads

$$w(\vec{x},\vec{n})\nu\frac{\partial {\cal I}}{\partial \nu} \... ...{{\cal I}_i - {\cal I}_{i-1}}{\Delta \nu} \mbox{\qquad} (w>0)$$ (17)

and

$$w(\vec{x},\vec{n})\nu\frac{\partial {\cal I}}{\partial \nu} \... ...{\cal I}_{i+1} - {\cal I}_i}{\Delta \nu} \mbox{\qquad} (w<0),$$ (18)

where $\Delta \nu$ is the constant frequency step size. All quantities referring to the discrete frequency point $\nu_i$ are denoted by an index "i''. Employing the Euler method, we get a semi-discrete representation of Eq. (16)

$$\Bigl({\cal A}_i^{\rm coh}+\frac{\vert w\vert\nu_i}{\Delta \... ... & (w > 0) \\ {\cal I}_{i+1} & (w < 0) \end{array} \right. .$$ (19)

The additional term on the left side of Eq. (19) can be interpreted as an artificial opacity, which is advantageous for the solution of the resulting linear system of equations. The additional term on the right side of Eq. (19) is included as an artificial source term. Equation (19) can be written in a compact operator form

$$\tilde{{\cal A}}_i^{\rm coh}{\cal I}_i=\tilde{f_i}.$$ (20)

Given a discretization with L degrees of freedom in $\Omega$, Mdirections on the unit sphere S² and N frequency points, the overall discrete system has the matrix form

$${\bf A}^{\rm coh} \vec{u} = \vec{f} ,$$ (21)

with the vector $\vec{u}$ containing the discrete intensities and the vector $\vec{f}$ the values of the source term for all frequency points $\nu_i$. Both vectors are of length $(L \cdot M \cdot N)$ and ${\bf A}^{\rm coh}$ is a $(L \cdot M \cdot N)\times(L \cdot M \cdot N)$ matrix. The fact that the function $w(\vec{x},\vec{n})$ may change its sign, results in a block-tridiagonal structure of  ${\bf A}^{\rm coh}$ and we get

$$\begin{array}{cccc} \left( \begin{array}{ccccc} {\tilde {\bf A}}_... ... \vec{f}_2 \\ \vdots \\ \vdots \\ \vec{f}_N \end{array}\right) . \end{array}$$

According to the sign of $w(\vec{x},\vec{n})$ the block matrices ${\bf R}_i$ and ${\bf B}_i$ hold entries of $w(\vec{x},\vec{n})\nu_i/\Delta\nu$ causing a redshift and blueshift of the photons in the medium, respectively. Requiring a reasonable resolution, the resulting linear system of equations of the total system is too large to be solved directly. Hence, we are treating N "monochromatic'' radiative transfer problems

$${\tilde {\bf A}}_i^{\rm coh} \vec{u}_i = \tilde{\vec{f}}_i ,$$ (22)

with a slightly modified right hand side

$$\tilde{\vec{f}}_i = \vec{f}_i + {\bf R}_i \vec{u}_{i+1} + {\bf B}_i \vec{u}_{i-1} .$$ (23)

Note that using simple velocity fields (e.g. infall or outflow) the sign of $w(\vec{x},\vec{n})$ is fixed and the matrix ${\bf A}^{\rm coh}$has a block-bidiagonal structure. This is favorable for our solution strategy, since we only need to solve Eq. (22) once for each frequency.

3.3 Discretization for complete redistribution

Considering complete redistribution, we use an implicit Euler scheme to discretize the Doppler term as described above and a simple quadrature method for the frequency integral in the scattering operator ${\cal S}^{\rm crd}$ in Eq. (12). Starting from N equidistantly distributed frequency points $\nu_i \in \bigl\{ \nu_1,\nu_2,...,\nu_N \bigr\} \subset \Lambda$ and N weights q₁,q₂,...,q_N, we define a quadrature method

$$Q(\nu_i) := \sum_{j=1}^N q_j \xi(\nu_j)$$ (24)

for integrals $\int_{\Lambda} \xi(\nu')~{\rm d}\nu'$. In the case of complete redistribution the kernel is

$$\xi(\nu_j) = \frac{\phi(\nu_j)}{4\pi}\int_{S^2}{\cal I}(\vec{x},\hat{\vec{n}},\nu_j) {\rm d}\hat{\omega} .$$ (25)

Separating the terms with the unknown intensities ${\cal I}_i$ from the known quantities ${\cal I}_j$ the discretized scattering integral (12) reads

$$\frac{\sigma_i}{4\pi} \phi_i q_i \int_{S^2} {\cal I}_i{\rm d}... ...^2} {\cal I}(\vec{x},\hat{\vec{n}},\nu_j){\rm d}\hat{\omega} .$$ (26)

Using this quadrature scheme and employing the Euler method for the frequency derivatives we get a semi-discrete formulation of the transfer problem for each frequency point including complete redistribution

 

$$\Bigl({\cal A}^{\rm crd}_i+ \frac{\vert w\vert\nu_i}{\Delta\nu}\B... ... \phi_jq_j \int_{S^2}{\cal I}(\vec{x},\hat{\vec{n}},\nu_j){\rm d}\hat{\omega} ,$$ (27)

where

$${\cal A}^{\rm crd}_i = {\cal T}+\chi_i+\phi_i q_i {\cal S}^{\rm coh} .$$ (28)

The additional terms on the right hand side of Eq. (27) must be interpreted as artificial source terms. Equation (27) can also be written in a compact operator form

$$\tilde{{\cal A}}_i^{\rm crd}{\cal I}_i=\hat{f_i}.$$ (29)

The total discrete system has the matrix form (cf. Eq. (21))

$${\bf A}^{\rm crd} \vec{u} = \vec{f} .$$ (30)

Unfortunately, the global frequency coupling via the scattering integral (12) results in a full block matrix and we get

$$\begin{array}{ccc} {\bf A}^{\rm crd} & = & \left( \begin{array}{c... ...& \ldots& \ldots&{\tilde {\bf A}}_N^{\rm crd} \end{array} \right) . \end{array}$$

${\bf B}_i$ and ${\bf R}_i$ again contain the contribution of the Doppler factor $w(\vec{x},\vec{n})\nu_i/\Delta\nu$, whereas ${\bf Q}_j$ holds the terms from the quadrature scheme. As we already explained in the case of coherent scattering, we do not solve the total system, but N"monochromatic'' radiative transfer problems

$${\tilde {\bf A}}_i^{\rm crd} \vec{u}_i = \hat{\vec{f}}_i ,$$ (31)

with a modified right hand side

$$\hat{\vec{f}}_i = \vec{f}_i + {\bf R}_i \vec{u}_{i+1} + {\bf B}_i \vec{u}_{i-1} + \sum_{j\not=i} {\bf Q}_j \vec{u}_j .$$ (32)

3.4 Full solution algorithm

Equations (20) and (29) are of the same form as the monochromatic radiative transfer equation, ${\cal A}{\cal I}=f$, for which we proposed a solution method based on a finite element technique in Paper I. The finite element method employs unstructured grids which are adaptively refined by means of an a posteriori error estimation strategy. Now, we apply this method to Eq. (20) or Eq. (29). In brief, the full solution algorithm reads:

1.

Start with ${\cal I}=0$ for all frequencies.

2.

Solve Eq. (20) or Eq. (29) for i=1,..,N.

3.

Repeat step 2 until convergence is reached.

4.

Refine grid and repeat step 2 and 3.

We start with a relatively coarse grid, where only the most important structures are pre-refined, and assure that the frequency interval $[\nu_1,\nu_N]$ is wide enough to cover the total line profile. Then, we solve Eq. (20) or Eq. (29) for each frequency several times depending on the choice of the redistribution function and the velocity field. During this fix point iteration, we monitor the changes of the resulting line profile in a particular direction $\vec{n}_{\rm out}$. Not until the line profile remains unchanged, we go over to step 4 and refine the spatial grid. Again, we apply the fixed fraction grid refinement strategy: the cells are ordered according to the size of the local refinement indicator $\eta_K=\max(\eta_K(\nu_i))\vert _{\nu_i}$ and a fixed portion of the cells with largest $\eta_K$ is refined. $\eta_K(\nu_i)$ is an indicator for the error of the solution in cell K at frequency $\nu_i$.

In Paper I (Sect. 2.3), we introduced various possibilities to determine $\eta_K$. If we use the total flux leaving the computational domain in direction $\vec{n}_{\rm out}$ as the quantity whose error functional is used to calculate $\eta_K$ (dual method), our convergence criterion described above is perfectly reasonable. We found that this criterion is also useful for the much simpler method (L² method), where $\eta_K$ is determined from the global error functional. In this case, the rate of convergence of the line profile in the directions $\vec{n} \ne \vec{n}_{\rm out}$ and for the line profile in direction $\vec{n}_{\rm out}$ are very similar. We are aware that this result may not be valid for more complex model configuration.

4 Model configurations

4.1 Halo and sources

We investigate spherically symmetric and elliptical distributions of the extinction coefficient $\chi(\vec{x})=\chi(x,y,z)$ of the form

$$\chi(\vec{x})=\left\{ \begin{array}{ll} \chi_0/(1+\alpha r_{... .../10^3 & \quad\mbox{for}\; r > r_{\rm h} \end{array} \right. ,$$ (33)

where r²=(x/a)²+(y/b)²+(z/c)². Note, that the value of the extinction coefficient is constant in the center within core radius $r_{\rm c}$ and outside halo radius $r_{\rm h}$. If not stated otherwise, $\chi_0$is determined from the line center optical depth

$$\tau=\int_{r_{\rm c}}^{r_{\rm h}}\chi(r)\phi(\nu_0)\;\vec{n}_{\rm thick}~{\rm d}\vec{x}$$ (34)

between $r_{\rm c}$ and $r_{\rm h}$ along the most optically thick direction $\vec{n}_{\rm thick}$. In total, the spatial distribution of $\chi$ is determined by six parameters: the length of the half-axes a, b and c, the radii $r_{\rm c}$ and $r_{\rm h}$, the dimensionless parameter $\alpha$, and the optical depth $\tau$. For $r_{\rm c}$, $r_{\rm h}$ and $\alpha$ we use the values given in Table 1.

Since we are predominantly interested in the transfer of radiation in resonance lines like Ly$\alpha$, we assume $\sigma(\vec{x})=\chi(\vec{x})$ and $\kappa(\vec{x})=0$ for all calculation presented here. This restricts us to the use of purely non-thermal source functions. In particular, we consider one or several spatially confined source regions with radius $r_{\rm s}$ each centered at a position $\vec{x}_i$:

$$f(\vec{x},\nu)=\left\{ \begin{array}{lr} \phi(\nu) & \quad\m... ...\vert\vec{x}-\vec{x}_i\vert > r_{\rm s} \end{array} \right. .$$ (35)

The function $\phi(\nu)$ is the Doppler profile defined in Eq. (2).

4.2 Velocity fields

In general, the finite element code is able to consider arbitrary velocity fields. For velocity fields defined on a discrete grid, for example for velocity fields resulting from hydrodynamical simulations, the velocity gradient in direction $\vec{n}$ must be obtained numerically. Here, we use two simple velocity fields and calculate the function w analytically.

The first velocity field describes a spherically symmetric inflow (v₀<0) or outflow (v₀>0) and is of the form

$$\vec{v}_{\rm io}=v_0 \left(\frac{r_0}{r}\right)^l \frac{\vec{x}}{r},$$ (36)

where $r=\vert\vec{x}\vert$ and v₀ the scalar velocity at radius r₀. The corresponding w function is

$$w(\vec{x},\vec{n})=v_0 \left(\frac{r_0}{r}\right)^l \left(\frac{1}{r}-(l+1) \frac{(\vec{n}\cdot\vec{x})^2}{r^3} \right)\cdot$$ (37)

For the second velocity field, we assume rotation around the z-axis

$$\vec{v}_{\rm rot}=v_0 \left(\frac{R_0}{R}\right)^l R^{-1} \left(\begin{array}{rrr}y\\ -x\\ 0\end{array}\right),$$ (38)

where R²=x²+y² is the distance from the rotational axis and v₀ the scalar velocity at distance R₀. If $\vec{n}=(n_x,n_y,n_z)$, the w function reads

$$w=v_0 \left(\frac{R_0}{R}\right)^l (l+1) \left(\frac{xy (n_y^2-n_x^2) + n_x n_y (x^2-y^2)}{R^3} \right)\cdot$$ (39)

In the special case of a Keplerian velocity field (l=0.5), we obtain

$$w=v_0 \frac{3}{2} R_0^{0.5} R^{-3.5} \left[x y (n_y^2-n_x^2) + n_x n_y (x^2-y^2)\right]$$ (40)

(see also Baschek & Wehrse 1999). Again, the values used for v₀, r₀ and R₀ are given in Table 1.

  
4.3 Normalization

Since we use a normalized form of the equation of transfer, where the computational domain is the unit cube, the results are valid for model configurations with different length scales and for different resonance lines. When $2x_{\max}$ is the size of the unit cube, coordinates in physical units are simply obtained by the transformation $\vec{x}\rightarrow \vec{x}\cdot x_{\rm max}$.

A particular line and the corresponding frequency scale must be defined when transforming to the observers frame. From the solution ${\cal I}(\vec{x},\vec{n},\nu)$ in the comoving frame, we obtain the solution $\tilde{{\cal I}}(\vec{x},\vec{n},\tilde{\nu})$ in the observers frame by performing the transformation

$$\tilde{{\cal I}}(x,\vec{n},\tilde{\nu})={\cal I}(\vec{x},\vec{n},\tilde{\nu}) \left(\frac{\tilde{\nu}}{\nu}\right)^3,$$ (41)

where

$$\tilde{\nu}=\nu\left(1+\vec{n}\cdot \frac{\vec{v}(\vec{x})}{c}\right)\cdot$$ (42)

The relation between optical depth $\tau$ and column density of the scattering media depends on the central line absorption cross section $\chi(\nu_0)$ and therefore on the Doppler velocity v_D. For the Ly$\alpha$ line of neutral hydrogen the central line absorption cross section is

$$\chi_{{\rm Ly}\alpha}=2.0\times10^{-15}\; \left(\frac{300\;\mbox{km}~\mbox{s}^{-1}}{v_{\rm D}}\right) \;\mbox{cm}^2.$$ (43)

Hence, the column density of neutral hydrogen for a halo with given $\tau$ is

$$N_{\rm H}=0.5\times10^{14}\;\tau\; \left(\frac{v_{\rm D}}{300\;\mbox{km}~\mbox{s}^{-1}}\right) \;\mbox{cm}^{-2}.$$ (44)

We would like to emphasize that this value is the column density of the neutral halo alone. Considering the column density of the source region, the total column density is about twice as high and is in the range $N_{\rm H}\sim [10^{13},10^{16}]\;\mbox{cm}^{-2}$ for the investigated configurations with $\tau=[0.1,100]$.

In comparison to column densities deduced from Voigt profile fitting procedures of Ly$\alpha$ profiles of high redshift radio galaxies (van Ojik et al. 1997), our values are at least an order of magnitude too low. With the present implementation of our finite element code, we could calculate systems with column densities up to $N_{\rm H}\sim10^{17}{-}10^{18}\;\mbox{cm}^{-2}$, but only with a large amount of computational time. Still higher column densities are beyond feasibility.

 

 

Table 1: Parameters used for all calculations. Distances are given in units of the computational cube (see Sect. 4.3).

$r_{\rm h}$	$r_{\rm c}$	$\alpha$	$r_{\rm s}$	$v_{\rm D}$	v0	r0	R0
1.0	0.2	103	0.2	10-3 c	-10-3 c	0.2	1.0

Figure 1: Ly$\alpha$ line profiles calculated with the FE code for a spherically symmetric model configuration: a) a static halo with coherent scattering, b) an infalling halo with coherent scattering, c) a static halo with complete redistribution, and d) an infalling halo with complete redistribution. For the static cases a) and c) the line styles refer to calculations with different optical depth $\tau$ as indicated. The small window in a) enlarges the peak of the line. The crosses mark the results of the analytical solution. For the moving halos we show in b) the results for $\tau =1$ (thin lines) and $\tau =10$ (thick lines) and in d) only for $\tau =10$. Here, the line styles refer to the exponent l used for the velocity fields.

Open with DEXTER

Figure 2: Results obtained for a static, disk-like model configuration: a) Emergent line profiles for different viewing angles and optical depths $\tau=\tau(\vec{n}_{\rm thick})$. b) Two-dimensional spectra for an edge-on view (90$^\circ$) and a nearly face-on view (20$^\circ$) and different optical depths. The position and width of the slit is indicated in Fig. 3a. The contours are given for 2.5%, 5%, 7.5%, 10%, 20%, ..., 90% of the maximum value.

Open with DEXTER

	Figure 3: Spatial intensity distribution for a static disk-like model configuration with $\tau(\vec{n}_{\rm thick})=10$ for different viewing angles: a) Frequency integrated intensity distribution. b) Intensity distribution for specific frequencies. The frequency is given in Doppler units relative to the line center $(\nu -\nu _0)/\Delta \nu _{\rm D}$ at the right hand side.
Open with DEXTER

5 Test calculations

We start with a spherically symmetric model configuration (a=b=c=1) and a single source region located at $\vec{x}=0$. Figure 1 shows the results of the finite element code (FE code) for different optical depths, velocity fields and redistribution functions. We used 41 frequencies equally spaced in the interval $(\nu-\nu_0)/\Delta\nu_{\rm D}=[-4,6]$ and 80 directions. Starting with a grid of 4³ cells and a pre-refined source region, we needed 3-5 spatial refinement steps.

The simplest case is a static model with coherent isotropic scattering. Figure 1a displays the emergent line profiles for different $\tau$. As expected, the Doppler profile is preserved and the flux $F_{\nu}$ is independent on $\tau$. The deviation of the numerical results from the analytical solution indicated with crosses is very small. The line profiles for $\tau=0.1$ and $\tau =1$ are identical even in the little window which shows the peak of the line in more detail. For $\tau=100$, the total flux is still conserved better than 99%. This result demonstrates that the frequency-dependent version of our FE code works correctly.

Next, we consider an infalling halo with coherent scattering and show the effects of frequency coupling due to the Doppler term. The emergent line profiles in Fig. 1b are plotted for different exponents l of the velocity field $\vec{v}_{\rm io}$ defined in Eq. (36). The line profiles are redshifted for $\tau =1$(thin lines). Most of the photons directly travel through the halo moving away from the observer. Since the Doppler term is proportional to the gradient of the velocity field, the redshift is larger for a greater exponent l. For $\tau =10$ (thick lines), the line profiles are blueshifted. Before photons escape from the optically thick halo in front of the source, they are back-scattered and blueshifted in the approaching halo behind the source. The blueshift is less pronounced for the accelerated infall with l=2, because the strong gradient of the velocity field leads to a slight redshift in the very inner parts of the halo. In this region, the total optical depth is still small. Further out, where the total optical depth increases, the line profile becomes blueshifted.

Complete redistribution is another method of frequency coupling which leads to a stronger coupling than the Doppler effect (see Sect. 3). The line profiles obtained for a static model with complete redistribution are displayed in Fig. 1c for different $\tau$. With increasing optical depth the photons more and more escape through the line wings. For $\tau\ge 1$, we get a double-peaked line profile with an absorption trough in the line center. The greater $\tau$ the larger the distance between the peaks and the depth of the absorption trough. Since our frequency resolution is too poor for the pointed wings, the flux conservation is only 96% for $\tau=100$.

Figure 1d shows the results for an infalling halo with complete redistribution for $\tau =10$ and different exponents l. For l=0 and l=0.5 the infalling motion of the halo enhances the blue wing of the double-peaked line profile. Equally, an outflowing halo would enhance the red peak. But for l=2, the red peak is slightly enhanced for an infalling halo due to the strong velocity gradient, as explained above. This example affords an insight into the mechanisms of resonance line formation and shows the necessity of a profound multi-dimensional treatment.

6 Applications

All calculations discussed in this section were performed with complete redistribution. We used 49 equidistant frequencies covering the interval $(\nu-\nu_0)/\Delta\nu_{\rm D}=[-6,6]$, 80 directions and started from a grid with 4³ cells and pre-refined source regions, which results in several 10³ cells for the initial mesh. 3-7 refinement steps were performed leading to approximately 10⁵ cells for the finest grid.

6.1 Elliptical halos

Figure 4: Results obtained for a disk-like model configuration with global infall motion (l=0.5): a) Emergent line profiles for different viewing angles and optical depths $\tau=\tau(\vec{n}_{\rm thick})$. b) Two-dimensional spectra for an edge-on view (90$^\circ$) and a nearly face-on view (20$^\circ$) and different optical depths. The position and width of the slit is indicated in Fig. 3a. The contours are given for 2.5%, 5%, 7.5%, 10%, 20%, ..., 90% of the maximum value.

Open with DEXTER

Figure 5: Results obtained for a disk-like model configuration with Keplerian rotation (l=0.5) around the z-axis: a) Emergent line profiles for different viewing angles and optical depths $\tau=\tau(\vec{n}_{\rm thick})$. b) Two-dimensional spectra for an edge-on view (90$^\circ$) and a nearly face-on view (20$^\circ$) and different optical depths. The position and width of the slit is indicated in Fig. 3a. The contours are given for 2.5%, 5%, 7.5%, 10%, 20%, ..., 90% of the maximum value.

Open with DEXTER

	Figure 6: Spatial intensity distribution for a static spherically symmetric model configuration with $\tau =10$ and three source regions for different viewing angles: a) Frequency integrated intensity distribution. b) Intensity distribution for specific frequencies. The frequency is given in Doppler units relative to the line center $(\nu -\nu _0)/\Delta \nu _{\rm D}$ at the right hand side.
Open with DEXTER

Figure 7: Results obtained for a spherically symmetric model configuration with $\tau =10$ and three source regions for the static case and two different velocity fields (inflow and rotation): a) Emergent line profiles for different viewing angles. b) Two-dimensional spectra for a nearly face-on view (70$^\circ$) and a nearly edge-on view (10$^\circ$). The position and width of the slit is indicated in Fig. 6a. The contours are given for 2.5%, 5%, 7.5%, 10%, 20%, ..., 90% of the maximum value.

Open with DEXTER

As a first step towards a full three-dimensional problem without any symmetries, we investigated an axially symmetric, disk-like model configuration (a:b:c=3:3:1) with a single source located at $\vec{x}=0$for several $\tau=\tau(\vec{n}_{\rm thick})$ and for different velocity fields. The most optical thick direction $\vec{n}_{\rm thick}$ is the direction within the equatorial plane of the disk. The direction perpendicular to the equatorial plane is the z-axis which we also call rotational axis even for cases without rotation.

By means of the static case we show which kind of information can be obtained from the results of our FE code. The emergent line profiles provide global information on the underlying system. They are displayed in Fig. 2a for different optical depths and viewing angles. The viewing angle is defined as the angle between the rotational axis and the direction towards the observer. For instance, a viewing angle of 90$^\circ$ means, the disk is "observed'' edge-on. Again, we obtain the characteristical double-peaked line profile for $\tau\ga 1$. The optical thickness decreases for smaller viewing angles. For a viewing angle of $0^\circ$the effective optical depth is only $\tau(\vec{n}_{\rm thick})/3$. Consequently, the distance of the peaks in the line profile and the depth of the absorption feature is growing with increasing viewing angle. Note, that the total relative flux $F(\vec{n}=0^{\circ})/F(\vec{n}=90^{\circ})$ escaping along the z-axis is increasing with increasing $\tau(\vec{n}_{\rm thick})$. Figure 2b will be explained later.

Far more information is contained in the spatial distribution of the line intensity projected into the s-t-plane which is the plane perpendicular to the viewing direction. In Fig. 3a the frequency-integrated intensity distribution is shown for $\tau =10$ and different viewing angles. At 90$^\circ$ the Ly$\alpha$ image shows an elliptical emission line region. The diffuse halo of scattered radiation is most clearly visible. With decreasing viewing angle the diffuse halo becomes less pronounced. For a nearly face-on view (20$^\circ$), the halo is optically thinner and the intensity distribution appears spherically symmetric.

Figure 3b displays the spatial distribution of the Ly$\alpha$ intensity for different viewing angles (columns) at selected frequencies (rows). For the edge-on view, the intensity distribution strongly depends on the frequency. In the outer line wings at $(\nu-\nu_0)/\Delta\nu_{\rm D}=\pm2$ the contours are elliptical and reproduce the disk-like form of the source region. The peaks of the line profile are at $(\nu-\nu_0)/\Delta\nu_{\rm D}\sim\pm1$. Here, the major axis of the elliptical intensity distribution in the inner parts of the model is parallel to the rotational axis, whereas the major axis of the elliptical contour lines in the outer parts remains parallel to the equatorial plane of the disk. In the very center of the line ( $(\nu-\nu_0)/\Delta\nu_{\rm D}=0$), the direct view to the source region is blocked by optically thick material. There appear two knots of Ly$\alpha$ emission directly above and below the disk, which are surrounded by an extended elliptical halo. The Ly$\alpha$ knot below the disk vanishes with decreasing viewing angle. At 20$^\circ$ the intensity distribution is spherically symmetric for all frequencies. The observable emission region is most extended at the frequency of the line center.

Two-dimensional spectra from high-resolution spectroscopy provide frequency-dependent data for only one spatial direction. To be able to compare our results with these observations we calculated two-dimensional spectra using the data within the slits indicated in Fig. 3a for the edge-on and nearly face-on view. The results are shown in Fig. 2b for different optical depths. We find that the form of the two-dimensional spectra is relatively insensible to the width of the slit. For $\tau =1$ a single emission region is visible. But already at 90$^\circ$, the highest contour lines reproduce the faint absorption trough of the corresponding line profile (Fig. 2a). The higher the optical depth the wider the gap and the spatial extent of the two peaks. Note that the spatial extent of the outer contour line only depends on $\tau(\vec{n}_{\rm thick})$ and not on the viewing direction.

What changes when imposing a macroscopic velocity field is shown in the following two figures. First, we consider an infalling velocity field with l=0.5 suitable for a gravitational collapse. Figure 4a displays the calculated line profiles for different $\tau$ and viewing angles. As expected, the blue peak of the line is enhanced. The higher $\tau$ the stronger the blue peak. Figure 4b shows the corresponding two-dimensional spectra obtained with the same slit width and position as in the static case. For low optical depth, the shape of the contour lines is a triangle. Photons changing frequency in order to escape via the blue wing are also scattered in space. The consequence is the greater spatial extent of the blue wing. Apart from a growing gap between the two peaks with higher optical depth, the general triangular shape in conserved.

Next, we investigated the elliptical model configuration with Keplerian rotation (l=0.5), where the z-axis is the axis of rotation. The results are plotted in Fig. 5 for $\tau =1$ and $\tau =10$. The emergent line profiles are symmetric with respect to the line center and show the same behavior with growing optical depth as in the static case. However, the extension of the line wings towards higher and lower frequencies is strongly increasing with increasing viewing angle because of the growing effect of the velocity field. Rotation is clearly visible in the two-dimensional spectra (Fig. 5b). The shear of the contour lines is the characteristical pattern indicating rotational motion. For an edge-on view at $\tau =10$, Keplerian rotation produces two banana-shaped emission regions.

In spite of the relatively low optical depth of the simple model configurations, our results reflect the form of line profiles and the patterns in two-dimensional spectra of many high redshift galaxies. For example, the two-dimensional spectra of the Ly$\alpha$blobs discovered by Steidel et al. (2000, Fig. 8) are comparable to the results obtained for infalling (Fig. 4b) and rotating (Fig. 5b) halos. In the sample of van Ojik et al. (1997) are many high redshift radio galaxies with single-peaked and double-peaked Ly$\alpha$ profiles. The corresponding two-dimensional spectra show asymmetrical emission regions which are more or less extended in space. De Breuck et al. (2000) find in their statistical study of emission lines from high redshift radio galaxies that the triangular shape of the Ly$\alpha$ emission is a characteristical pattern in the two-dimensional spectra of high redshift radio galaxies. Since the emission of the blue peak of the line profile is predominately less pronounced, most of the associated halos should be in the state of expansion.

6.2 Multiple sources

High redshift radio galaxies are found in the center of proto clusters. In such an environment, it could be possible that the Ly$\alpha$ emission of several galaxies is scattered in a common halo. To investigate this scenario, we started with a spherically symmetric distribution for the extinction coefficient and with three source regions located at x₁=[0.5,0.25,0], x₂=[-0.5,0.25,0] and x₃=[0,-0.25,0] forming a triangle in the x-y-plane. In the following figures, the results for $\tau =10$ are presented.

We begin with the static model. Figure 6a shows Ly$\alpha$ images for four selected viewing directions. The specified angle is the angle between the viewing direction and the x-y-plane. Note that the orientation of the source positions within the plane is different for each image. Viewing the configuration almost perpendicular to the x-y-plane (70$^\circ$), renders all three source regions visible, because the source regions are situated in the more optically thin, outer parts of the halo. Remember that most of the scattering matter is in and around the center of the system. For other angles, some of the source regions are located behind the center and therefore less visible on the images. Figure 6b shows images at different frequencies. In the line wings at $(\nu-\nu_0)/\Delta\nu_{\rm D}=\pm2$, the number and position of the source regions can be determined for all viewing angles. However, at frequencies around the line center at $(\nu-\nu_0)/\Delta\nu_D=\pm1\;\mbox{and}\;0$, the number of visible sources depends on the viewing angle. Some of the images show only one, other two or three sources.

The corresponding line profiles and two-dimensional spectra are displayed in Fig. 7 for the static case as well as for a halo with global inflow and Keplerian rotation. Width and position of the slits are depicted in Fig. 6a. We get double-peaked line profiles for almost all cases, except for the rotating halo, where the line profiles are very broad for viewing angles lower than 70$^\circ$. Additional features, dips or shoulders, are visible in the red or blue wing. They arise because the three sources have significantly different velocities components with respect to the observer.

The slit for a viewing angle of 70$^\circ$ contains two sources. They show up as four emission regions in the two-dimensional spectra (Fig. 7b). The pattern is very symmetric, even for the moving halos. For a viewing angle of 10$^\circ$, the slit covers all source regions. Nevertheless, the two-dimensional spectra are dominated by two pairs of emission regions resulting from the sources located closer to the observer. The third source region only shows up as a faint emission in the blue part in the case of global inflow. In the case of rotation, emission regions from a third source are present. But the overall pattern is very irregular and prevents a clear identification of the emission regions.

This example demonstrates the complexity of three dimensional problems. In a clumpy medium, the determination of the number and position of Ly$\alpha$ sources would be difficult by means of frequency integrated images alone. Two-dimensional spectra may help, but could prove to be too complicated. More promising are images obtained from different parts of the line profile or information from other emission lines of OIII or H$\alpha$ in a manner proposed by Kurk et al. (2001).

7 Summary

We presented a finite element code for solving the resonance line transfer problem in moving media. Non-relativistic velocity fields and complete redistribution are considered. The code is applicable to any three-dimensional model configuration with optical depths up to 10^3-4.

We showed applications to the hydrogen Ly$\alpha$ line of slightly optically thick model configurations ( $\tau\le 10^2$) and discussed the resulting line profiles, Ly$\alpha$ images and two-dimensional spectra. The systematic approach from very simple to more complex models gave the following results:

• An optical depth of $\tau\ga 1$ leads to the characteristic double peaked line profile with a central absorption trough as expected from analytical studies (e.g. Neufeld 1990). This form of the profile is the result of scattering in space and frequency. Photons escape via the line wings where the optical depth is much lower.
• Global velocity fields destroy the symmetry of the line profile. Generally, the blue peak of the profile is enhanced for models with infall motion and the red peak for models with outflow motion. But there are certain velocity fields (e.g. with steep gradients) and spatial distributions of the extinction coefficient, where the formation of a prominent peak is suppressed.
• Double-peaked line-profiles show up as two emission regions in the two-dimensional spectra. Global infall or outflow leads to an overall triangular shape of the emission. Rotation produces a shear pattern resulting in banana-shaped emission regions for optical depths $\ga$10.
• For non-symmetrical model configurations, the optical depth varies with the line of sight. Thus, the total flux, the depth of the absorption trough and the pattern in the two-dimensional spectra strongly depend on the viewing direction.

The applications demonstrate the capacity of the finite element code and show that the three-dimensional structure and the kinematics of the model configurations are very important. Thus, beside exploring higher optical depths up to the limits of our code, we will consider clumpy density distributions as well as dust absorption and try to model the Ly$\alpha$ emission of individual high redshift galaxies. In addition, we intend to implement a second order method for the frequency discretization. Furthermore, we plan to overcome the difficulties in solving the line transfer problem for configurations with $\tau>10^4$. Therefore, we will extend our method and use the diffusion approximation in the most optically thick regions. In the other regions, the full frequency-dependent line transfer problem must be solved.

Acknowledgements

We would like to thank Rainer Wehrse and Guido Kanschat for useful discussions. This work is supported by the Deutsche Forschungsgemeinschaft (DFG) within the SFB 359 "Reactive Flows, Diffusion and Transport'' and SFB 439 "Galaxies in the young Universe''.

References

 

• Adams, T. F. 1972, ApJ, 174, 439 In the text NASA ADS  
• Ahn, S.-H., Lee, H.-W., & Lee, H. M. 2001, ApJ, 554, 604 In the text NASA ADS  
• Ahn, S.-H., Lee, H.-W., & Lee, H. M. 2002, ApJ, 567, 992 In the text  
• Baschek, B., & Wehrse, R. 1999, Phys. Rep., 311, 201 In the text NASA ADS  
• Bicknell, G. V., Sutherland, R. S., van Breugel, W. J. M., et al. 2000, ApJ, 540, 678 In the text NASA ADS  
• De Breuck, C., Röttgering, H., Miley, G., van Breugel, W., & Best, P. 2000, A&A, 362, 519 In the text NASA ADS  
• Fynbo, J. U., Thomsen, B., & Møller, P. 2000, A&A, 353, 457 In the text NASA ADS  
• Fynbo, J. U., Møller, P., & Thomsen, B. 2001, A&A, 374, 443 In the text NASA ADS  
• Hippelein, H., & Meisenheimer, K. 1993, Nature, 362, 224 In the text NASA ADS  
• Hu, E. M., Cowie, L. L., & McMahon, R. G. 1998, ApJ, 502, L99 In the text NASA ADS  
• Hummer, D. G., & Kunasz, P. B. 1980, ApJ, 236, 609 In the text NASA ADS  
• Kanschat, G. 1996, Ph.D. Thesis, Univ. of Heidelberg In the text  
• Kurk, J. D., Röttgering, H. J. A., Miley, G. K., & Pentericci, L. 2001, RMA&A, Ser. Conf., in press In the text  
• Kudritzki, R.-P., Mendez, R. H., Feldmeier, J. J., et al. 2000, ApJ, 536, 19 In the text NASA ADS  
• Mihalas, D., & Weibel-Mihalas, B. 1984, Foundation of Radiation Hydrodynamics (Oxford University Press, New York) In the text  
• Neufeld, D. A. 1990, ApJ, 350, 216 In the text NASA ADS  
• Richling, S., Meinköhn, E., Kryzhevoi, N., & Kanschat, G. 2001, A&A, 380, 776 (Paper I) In the text NASA ADS  
• Rhoads, J. E., Malhotra, S., & Dey, A. 2000, ApJ, 545, L85 In the text NASA ADS  
• Steidel, C. C., Adelberger, K. L., Sharply, A. E., et al. 2000, ApJ, 532, 170 In the text NASA ADS  
• van Ojik, R., Röttgering, H. J. A., Carilli, C. L., et al. 1996, A&A, 313, 25 In the text NASA ADS  
• van Ojik, R., Röttgering, H. J. A., Miley, G. K., & Hunstead, R. W. 1997, A&A, 317, 358 In the text NASA ADS  
• Villar-Martín, M., Binette, L., & Fosbury, R. A. E. 1999, A&A, 346, 7 In the text NASA ADS  
• Wehrse, R., Baschek, B., & von Waldenfels, W. 2000, A&A, 359, 780 In the text NASA ADS