A&A 392, 827-839 (2002)
E. Meinköhn1,2 - S. Richling1
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
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 103-4. Results of Ly 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 line profile for all models with an optical depth 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 spectra, whereas rotation creates a shear pattern. Frequency-resolved Ly images may help to find the number and position of multiple Ly sources located in a single halo. A qualitative comparison with observations of extended Ly 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
Hydrogen Ly 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 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 emission also bears information on the distribution and kinematics of the interstellar gas as well as the nature of the energy source. Ly 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 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 .
The interpretation of Ly 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 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.
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
operator form simply reads
The extinction coefficient
is the sum of the
frequency-dependence is given by a
normalized profile function .
Usually, we adopt a Doppler
For the source term
The "Doppler'' operator
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
The scattering operator
depends on the general redistribution
giving the probability that
is absorbed and a photon
emitted. In the following, we assume isotropic scattering
For the modeling of prominent resonance lines, in particular
we restrict the frequency discretization to a finite
and are located far out of the line center in the continuum. The function
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
In order to solve the radiative transfer Eq. (1), we
first carry out the frequency discretization including coherent
scattering. With the abbreviation
Equations (20) and (29) are of the same form as the monochromatic radiative transfer equation, , 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:
In Paper I (Sect. 2.3), we introduced various possibilities to determine . If we use the total flux leaving the computational domain in direction as the quantity whose error functional is used to calculate (dual method), our convergence criterion described above is perfectly reasonable. We found that this criterion is also useful for the much simpler method (L2 method), where is determined from the global error functional. In this case, the rate of convergence of the line profile in the directions and for the line profile in direction are very similar. We are aware that this result may not be valid for more complex model configuration.
Since we are predominantly interested in the transfer of radiation in
resonance lines like Ly,
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
each centered at a position :
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 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 (v0<0) or outflow (v0>0) and is of the form
A particular line and the corresponding frequency scale must be
defined when transforming to the observers frame. From the solution
in the comoving frame, we obtain the solution
in the observers frame by performing
In comparison to column densities deduced from Voigt profile fitting
procedures of Ly
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
but only with a
large amount of computational time. Still higher column densities are
|1.0||0.2||103||0.2||10-3 c||-10-3 c||0.2||1.0|
|Figure 1: Ly 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 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 (thin lines) and (thick lines) and in d) only for . 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 . b) Two-dimensional spectra for an edge-on view (90) and a nearly face-on view (20) 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 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 at the right hand side.|
|Open with DEXTER|
We start with a spherically symmetric model configuration (a=b=c=1) and a single source region located at . 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 and 80 directions. Starting with a grid of 43 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 . As expected, the Doppler profile is preserved and the flux is independent on . The deviation of the numerical results from the analytical solution indicated with crosses is very small. The line profiles for and are identical even in the little window which shows the peak of the line in more detail. For , 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 defined in Eq. (36). The line profiles are redshifted for (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 (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 . With increasing optical depth the photons more and more escape through the line wings. For , we get a double-peaked line profile with an absorption trough in the line center. The greater 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 .
Figure 1d shows the results for an infalling halo with complete redistribution for 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.
All calculations discussed in this section were performed with complete redistribution. We used 49 equidistant frequencies covering the interval , 80 directions and started from a grid with 43 cells and pre-refined source regions, which results in several 103 cells for the initial mesh. 3-7 refinement steps were performed leading to approximately 105 cells for the finest grid.
|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 . b) Two-dimensional spectra for an edge-on view (90) and a nearly face-on view (20) 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 . b) Two-dimensional spectra for an edge-on view (90) and a nearly face-on view (20) 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 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 at the right hand side.|
|Open with DEXTER|
|Figure 7: Results obtained for a spherically symmetric model configuration with 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) and a nearly edge-on view (10). 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 for several and for different velocity fields. The most optical thick direction 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 means, the disk is "observed'' edge-on. Again, we obtain the characteristical double-peaked line profile for . The optical thickness decreases for smaller viewing angles. For a viewing angle of the effective optical depth is only . 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 escaping along the z-axis is increasing with increasing . 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 and different viewing angles. At 90 the Ly 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), the halo is optically thinner and the intensity distribution appears spherically symmetric.
Figure 3b displays the spatial distribution of the Ly 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 the contours are elliptical and reproduce the disk-like form of the source region. The peaks of the line profile are at . 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 ( ), the direct view to the source region is blocked by optically thick material. There appear two knots of Ly emission directly above and below the disk, which are surrounded by an extended elliptical halo. The Ly knot below the disk vanishes with decreasing viewing angle. At 20 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 a single emission region is visible. But already at 90, 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 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 and viewing angles. As expected, the blue peak of the line is enhanced. The higher 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 and . 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 , 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 Lyblobs 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 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 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.
High redshift radio galaxies are found in the center of proto clusters. In such an environment, it could be possible that the Ly 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 x1=[0.5,0.25,0], x2=[-0.5,0.25,0] and x3=[0,-0.25,0] forming a triangle in the x-y-plane. In the following figures, the results for are presented.
We begin with the static model. Figure 6a shows Ly 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), 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 , the number and position of the source regions can be determined for all viewing angles. However, at frequencies around the line center at , 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. 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 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, 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 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 in a manner proposed by Kurk et al. (2001).
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 103-4.
We showed applications to the hydrogen Ly line of slightly optically thick model configurations ( ) and discussed the resulting line profiles, Ly images and two-dimensional spectra. The systematic approach from very simple to more complex models gave the following results:
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''.