next previous
Up: Metal abundances and kinematics


3 Computational procedure

The complete description of our computational procedure is given in Paper I. Here we summarize briefly basic model assumptions and emphasize new details recently included in the MCI.

We assume that all lines observed in a metal system arise in a continuous absorbing gas slab of a thickness L (presumably the outer region of a foreground distant galaxy). The absorber exhibits a fluctuating gas density and a mixture of bulk motions such as infall and outflows, tidal flows etc., resulting in a stochastic velocity field. Metal abundances are assumed to be constant within the absorber and gas is supposed to be optically thin for the ionizing UV radiation.

Within the absorbing region the radial velocity v(s) and total hydrogen density $n_{\rm H}(s)$ distributions along the line of sight are the same for all ions. In the computational procedure these two random fields are represented by their sampled values at equally spaced intervals $\Delta x$ (x is dimensionless radial coordinate s/L), i.e. by the vectors $\{ v_1, \ldots, v_k \}$ and $\{ n_1, \ldots, n_k \}$ with k large enough to describe the narrowest components of complex spectral lines.

Further we suppose the thermodynamic and ionization equilibrium at each computational point along the sightline which means that fractional ionizations of different ions are determined exclusively by the gas density and vary from point to point. These fractional ionization variations are just the cause of the observed diversity of profile shapes. To calculate the kinetic temperature and fractional ionization of ions the photoionization code CLOUDY (Ferland 1997) was used.

The inputs to CLOUDY are the dimensionless ionization parameter $U = n_{\rm ph}/n_{\rm H}$ ( $n_{\rm ph}$ - the number density of photons with energies above 1 Ry), metallicity and the background ionizing spectrum for which the Haardt-Madau spectrum (HM) was adopted (Haardt & Madau 1996). The number density of the ionizing photons for this spectrum

 \begin{displaymath}%
n_{\rm ph} = \frac{4\pi}{c\,h}J_{912}\,\int^\infty_{\nu_{\r...
...}\,
\left(\frac{J_\nu}{J_{912}}\right) \frac{{\rm d}\nu}{\nu},
\end{displaymath} (1)

equals at z = 2 to $2.26\times10^{-5}$ cm-3 [here c, h, $\nu_{\rm c}$ and $J_\nu$ are the speed of light, the Planck constant, the frequency of the hydrogen Lyman edge, and the specific intensity (in ergs cm-2 s-1 sr-1 Hz-1) at frequency $\nu$, respectively].

Fractional ionization curves $\Upsilon (U)$ were computed with CLOUDY for solar abundance pattern and different metallicities and then included in the MCI code to calculate the optical depths for ions involved in the fitting. If the obtained solution revealed the abundance pattern different from the solar one, $\Upsilon (U)$ were recalculated for this new pattern and all computations repeated. It should be noted, however, that differences of $0.2\div0.3$ dex from solar values influence the fractional ionizations only weakly.

The values of velocity and density at subsequent computational points are considered to be correlated and are described by means of Markovian processes. In particular, the velocity is computed as follows:

 \begin{displaymath}%
v(x + \Delta x) = f_{\rm v}\,v(x) + \epsilon(x + \Delta x)\; ,
\end{displaymath} (2)

where $f_{\rm v} = R_{\rm v}(\Delta x)/\sigma^2_{\rm v}$, $R_{\rm v}$ being the correlation between the velocity values at points separated by a distance $\Delta x$, i.e. $R_{\rm v} = \langle v(x + \Delta x)\,v(x) \rangle$, $\sigma_{\rm v}$ the velocity dispersion of bulk motions, and $\epsilon$ a random normal variable with zero mean and dispersion

 \begin{displaymath}%
\sigma_{\epsilon,{\rm v}} = \sigma_{\rm v}\sqrt{1 - f^2_{\rm v}}.
\end{displaymath} (3)


  \begin{figure}
\par\includegraphics[height=12cm,width=8.6cm,clip]{H3203F2.PS}\end{figure} Figure 2: Computed velocity (panel a) and density of gas (panel  b) and ions (panels c-j) distributions along the line of sight for the system at z = 1.92595toward J2233-606. Shown are patterns rearranged according to the principle of minimum entropy production rate (see text).


  \begin{figure}
\par\includegraphics[height=12cm,width=8.4cm,clip]{H3203F3.PS}\end{figure} Figure 3: Density-weighted velocity distribution functions, p(v), for H I, C II, Si II, N III, Si III, C IV, and Si IV as restored by the MCI procedure in the z = 1.92595 system toward J2233-606.

Density is supposed to have a log-normal distribution (characterized by the mean n0 and the second central moment $\sigma_{\rm y}$ of dimensionless variable y = $n_{\rm H}$/n0) and is described with help of an auxiliary Markovian process $\nu (x)$:

 \begin{displaymath}%
\nu(x + \Delta x) = f_{\nu}\,\nu(x) + \epsilon(x + \Delta x),
\end{displaymath} (4)

where $f_{\nu} = R_\nu(\Delta x)/\sigma^2_\nu$, $R_\nu$ being the correlation between the values of $\nu$ at points separated by a distance $\Delta x$, $\sigma_\nu$ the logarithmic density dispersion, $\sigma_\nu = \sqrt{\ln(1 + \sigma^2_{\rm y})}$ and $\epsilon$ a random normal variable with zero mean and dispersion $\sigma_{\epsilon,\nu} = \sigma_\nu\sqrt{1 - f^2_{\nu}}$ .

Having defined $\nu (x)$, the total hydrogen density can be obtained as

 \begin{displaymath}%
n_{\rm H}(x) = n_0\,{\rm e}^{\nu(x) - \frac{1}{2}\sigma^2_\nu}.
\end{displaymath} (5)

The optical depth $\tau_{{\rm a},i}(\lambda)$ at wavelength $\lambda$ for element `a' in ith ionizing stage is calculated from the equation:
 
$\displaystyle %
{
\tau_{{\rm a},i}(\lambda) =
k_0 Z_{\rm a} N_0\,{\rm e}^{-\fra...
...ilon_{{\rm a},i}[U(x)]\,
\Phi_\lambda[\Delta\lambda_{\rm D}(x),v(x)]\,{\rm d}x,$ (6)

where k0 is a constant, $N_0 = n_0\,L$ - the expectation value of the total hydrogen column density, $Z_{\rm a} = n_{\rm a}/n_{\rm H}$ - the relative abundance of element "a'', $\Upsilon_{{\rm a},i}[U(x)]$ - the fractional ionization for stage i, $\Phi_\lambda$ - the profile function, $\Delta\lambda_{\rm D}$ - the Doppler width of the line.

Before being compared with the observed spectrum, the synthetic intensities $\exp[-\tau (\lambda )]$ are convolved with the spectrograph point-spread function.

Thus the proposed model is fully defined by specifying the following values: the velocity vector $\{v_i\}$, the total hydrogen density vector $\{n_i\}$, the total hydrogen column density N0, the mean ionization parameter U0, the radial velocity dispersion $\sigma_{\rm v}$, the density second central moment $\sigma_{\rm y}$, the element abundances $Z_{\rm a}$ and the correlation coefficients $f_{\rm v}$ and $f_\nu$. The common least-squares minimization (LSM) of the objective function is used to estimate the model parameters[*].

The computational procedure itself consists of two steps: firstly a point in the parameter space ( $N_0, U_0, \sigma_{\rm v}, \sigma_{\rm y},
Z_{\rm a})$ is chosen and then an optimal configuration of $\{v_i\}$ and $\{n_i\}$ for this parameter set is searched for. Correlation coefficients are considered as external parameters and remain fixed during the calculations.


  \begin{figure}
\par\includegraphics[height=13cm,width=17.6cm,clip]{H3203F4.PS}\end{figure} Figure 4: Same as Fig. 1 but for the $z_{\rm abs} \simeq 1.94$ system toward J2233-606. The zero radial velocity is fixed at z = 1.942616. The corresponding physical parameters are listed in Table 1. The normalized $\chi ^2_{\rm min} = 1.08$ (the number of degrees of freedom $\nu = 888$).

The optimization of $\{v_i\}$ and $\{n_i\}$ is the most time-consuming part of the procedure and needs an effective algorithm to achieve a quick and stable convergence of the computations. In the MCI we use the simulated annealing with Tsallis acceptance rule (Xiang et al. 1997) and an adaptive annealing temperature choice. Namely, the annealing temperature $T_{\rm a}$ at iteration k+1 is decreased according to following equation:

 \begin{displaymath}%
T^{k+1}_{\rm a} = \frac{ T^k_{\rm a} }{1 + \beta
\frac{\Delta E_k}{E_k} r},
\end{displaymath} (7)

where $\Delta E_k/E_k$ is the relative energy $(\equiv$ the sum of squares of the deviations, ${\cal S}$) variation on step k, r is the acceptance rate, i.e. the ratio of the accepted trials to the total trial number, and $\beta$ is a constant of order 1.

The calculation of the uncertainty ranges for the fitting parameters is in our method not so straightforward as a simple inversion of the Hesse matrix since the velocity and density distributions represent additional degrees of freedom and widen, in general, the confidence intervals. However, these $\{v_i\}$ and $\{n_i\}$distributions themselves are nuisance parameters and should be "integrated out'' when one computes the errors for the other parameters. To estimate the confidence levels, the following procedure can be applied: the values of the physical parameters in the vicinity of the global minimum of the objective function are chosen at random and then the optimal density and velocity distributions are computed. Assuming that the probability of each parameter set can be linked to the derived ${\cal S}$ value [e.g. as $\exp(-{\cal S})$] we can estimate from the obtained sample the joint probability density function for parameters and hence calculate all necessary statistical moments.

It should be noted, however, that the reliable estimation of this multidimensional function requires a very large sample which is quite time consuming. In our case the exact estimation of confidence levels is not very crucial taking into account the intrinsic uncertainties in atomic data (e.g. Savin 2000) or unknown shape and intensity of the local background ionizing radiation. Because of this we restricted our samples to a few dozens of points and estimated the accuracy of the fitting parameters only approximately.

The recovered density and velocity patterns are not unique - many configurations are possible with comparable probability. But all these configurations have the same density-weighted velocity distributions which actually determine the observed line shapes (see Paper I). As already mentioned above, we represent these random fields by their values sampled at equally spaced intervals $\Delta x$. In order to compare the calculated patterns we rearrange these values in such a way that the final configuration exhibits a lowest rate of entropy production: according to the Prigogine theorem (Prigogine 1967), this configuration has the minimal dissipation and, hence, is more stable and more probable. All necessary equations to calculate the entropy production are presented in Appendix. We stress, however, that configurations produced on the base of such rearrangement should in no case be considered as something final - they represent only the (most) probable case of the density and velocity distributions along the line of sight and are used here exclusively for illustrations.


next previous
Up: Metal abundances and kinematics

Copyright ESO 2002