text:
# Measurement of the 𝐷⁺→𝐾̄^{∗0}𝑙⁺𝜈_𝑙 Branching Fraction

## Abstract

Using $`13.53`$ $`\mathrm{fb}^1`$ of CLEO data, we have measured the ratios of the branching fractions $`R_e^+=\frac{(D^+\overline{K}^0e^+\nu _e)}{(D^+K^{}\pi ^+\pi ^+)}`$, $`R_\mu ^+=\frac{(D^+\overline{K}^0\mu ^+\nu _\mu )}{(D^+K^{}\pi ^+\pi ^+)}`$ and the combined branching fraction ratio $`R_l^+=\frac{(D^+\overline{K}^0l^+\nu _l)}{(D^+K^{}\pi ^+\pi ^+)}`$. We find $`R_e^+`$ $`=`$ $`0.74\pm 0.04\pm 0.05`$, $`R_\mu ^+`$ $`=`$ $`0.72\pm 0.10\pm 0.06`$ and $`R_l^+`$ $`=`$ $`0.74\pm 0.04\pm 0.05`$, where the first and second errors are statistical and systematic, respectively. The known branching fraction $`(D^+K^{}\pi ^+\pi ^+)`$ leads to: $`(D^+\overline{K}^0e^+\nu _e)`$ $`=`$ $`(6.7\pm 0.4\pm 0.5\pm 0.4)\%`$, $`(D^+\overline{K}^0\mu ^+\nu _\mu )`$ $`=`$ $`(6.5\pm 0.9\pm 0.5\pm 0.4)\%`$ and $`(D^+\overline{K}^0l^+\nu _l)`$ $`=`$ $`(6.7\pm 0.4\pm 0.5\pm 0.4)\%`$, where the third error is due to the uncertainty in $`(D^+K^{}\pi ^+\pi ^+)`$.

preprint:   CLNS 02-1776  CLEO 02-02

G. Brandenburg,1 A. Ershov,1 D. Y.-J. Kim,1 R. Wilson,1 K. Benslama,2 B. I. Eisenstein,2 J. Ernst,2 G. D. Gollin,2 R. M. Hans,2 I. Karliner,2 N. Lowrey,2 M. A. Marsh,2 C. Plager,2 C. Sedlack,2 M. Selen,2 J. J. Thaler,2 J. Williams,2 K. W. Edwards,3 R. Ammar,4 D. Besson,4 X. Zhao,4 S. Anderson,5 V. V. Frolov,5 Y. Kubota,5 S. J. Lee,5 S. Z. Li,5 R. Poling,5 A. Smith,5 C. J. Stepaniak,5 J. Urheim,5 S. Ahmed,6 M. S. Alam,6 L. Jian,6 M. Saleem,6 F. Wappler,6 E. Eckhart,7 K. K. Gan,7 C. Gwon,7 T. Hart,7 K. Honscheid,7 D. Hufnagel,7 H. Kagan,7 R. Kass,7 T. K. Pedlar,7 J. B. Thayer,7 E. von Toerne,7 T. Wilksen,7 M. M. Zoeller,7 S. J. Richichi,8 H. Severini,8 P. Skubic,8 S.A. Dytman,9 S. Nam,9 V. Savinov,9 S. Chen,10 J. W. Hinson,10 J. Lee,10 D. H. Miller,10 V. Pavlunin,10 E. I. Shibata,10 I. P. J. Shipsey,10 D. Cronin-Hennessy,11 A.L. Lyon,11 C. S. Park,11 W. Park,11 E. H. Thorndike,11 T. E. Coan,12 Y. S. Gao,12 F. Liu,12 Y. Maravin,12 I. Narsky,12 R. Stroynowski,12 J. Ye,12 M. Artuso,13 C. Boulahouache,13 K. Bukin,13 E. Dambasuren,13 R. Mountain,13 T. Skwarnicki,13 S. Stone,13 J.C. Wang,13 A. H. Mahmood,14 S. E. Csorna,15 I. Danko,15 Z. Xu,15 G. Bonvicini,16 D. Cinabro,16 M. Dubrovin,16 S. McGee,16 A. Bornheim,17 E. Lipeles,17 S. P. Pappas,17 A. Shapiro,17 W. M. Sun,17 A. J. Weinstein,17 G. Masek,18 H. P. Paar,18 R. Mahapatra,19 R. A. Briere,20 G. P. Chen,20 T. Ferguson,20 G. Tatishvili,20 H. Vogel,20 N. E. Adam,21 J. P. Alexander,21 C. Bebek,21 K. Berkelman,21 F. Blanc,21 V. Boisvert,21 D. G. Cassel,21 P. S. Drell,21 J. E. Duboscq,21 K. M. Ecklund,21 R. Ehrlich,21 L. Gibbons,21 B. Gittelman,21 S. W. Gray,21 D. L. Hartill,21 B. K. Heltsley,21 L. Hsu,21 C. D. Jones,21 J. Kandaswamy,21 D. L. Kreinick,21 A. Magerkurth,21 H. Mahlke-Krüger,21 T. O. Meyer,21 N. B. Mistry,21 E. Nordberg,21 M. Palmer,21 J. R. Patterson,21 D. Peterson,21 J. Pivarski,21 D. Riley,21 A. J. Sadoff,21 H. Schwarthoff,21 M. R. Shepherd,21 J. G. Thayer,21 D. Urner,21 B. Valant-Spaight,21 G. Viehhauser,21 A. Warburton,21 M. Weinberger,21 S. B. Athar,22 P. Avery,22 H. Stoeck,22 and J. Yelton22

1Harvard University, Cambridge, Massachusetts 02138

2University of Illinois, Urbana-Champaign, Illinois 61801

3Carleton University, Ottawa, Ontario, Canada K1S 5B6

and the Institute of Particle Physics, Canada M5S 1A7

4University of Kansas, Lawrence, Kansas 66045

5University of Minnesota, Minneapolis, Minnesota 55455

6State University of New York at Albany, Albany, New York 12222

7Ohio State University, Columbus, Ohio 43210

8University of Oklahoma, Norman, Oklahoma 73019

9University of Pittsburgh, Pittsburgh, Pennsylvania 15260

10Purdue University, West Lafayette, Indiana 47907

11University of Rochester, Rochester, New York 14627

12Southern Methodist University, Dallas, Texas 75275

13Syracuse University, Syracuse, New York 13244

14University of Texas - Pan American, Edinburg, Texas 78539

15Vanderbilt University, Nashville, Tennessee 37235

16Wayne State University, Detroit, Michigan 48202

17California Institute of Technology, Pasadena, California 91125

18University of California, San Diego, La Jolla, California 92093

19University of California, Santa Barbara, California 93106

20Carnegie Mellon University, Pittsburgh, Pennsylvania 15213

21Cornell University, Ithaca, New York 14853

22University of Florida, Gainesville, Florida 32611

The transition amplitude of the semileptonic decay $`D^+\overline{K}^0l^+\nu _l`$ is proportional to the product of the leptonic and hadronic currents. The hadronic current is represented by three analytic functions called form factors, $`A_1(q^2)`$, $`A_2(q^2)`$, and $`V(q^2)`$ , where $`q^2`$ is the invariant mass squared of the virtual $`W`$ boson. The form factors cannot be easily computed in quantum chromodynamics (QCD) since they are affected by significant non-perturbative contributions. Further, unlike the case of $`BD^{}l\nu `$, Heavy Quark Effective Theory (HQET)  cannot be applied in a straightforward way since the transition is to the light strange quark. Precise experimental measurements are needed to guide theoretical progress in this area. Measurements of these form factors are also valuable since, using HQET, they can be related to the form factors for the important decay $`bul\nu `$, a heavy to light transition like the $`csl\nu `$ process we study here. This approach could reduce theoretical uncertainties on the value of $`|V_{ub}|`$ measured in $`bul\nu `$.

In this work, the ratio of the branching fraction of $`D^+\overline{K}^0l^+\nu _l`$ to that of $`D^+K^{}\pi ^+\pi ^+`$ ($`R_l^+`$) was measured, where $`l`$ stands for either an electron or a muon. In other experiments , the form factor ratios $`r_2=\frac{A_2}{A_1}`$ and $`r_V=\frac{V}{A_1}`$ of $`D^+\overline{K}^0l^+\nu _l`$ were obtained from angular correlations between the daughter particles. Combining these, we can calculate the form factors $`A_1`$, $`A_2`$ and $`V`$.

The data used in this study was collected with the two configurations of the CLEO detector  at the Cornell Electron Storage Ring (CESR). It consists of $`9.13`$ $`\mathrm{fb}^1`$ of integrated luminosity on the $`\mathrm{{\rm Y}}(4S)`$ resonance and $`4.40`$ $`\mathrm{fb}^1`$ below $`B\overline{B}`$ threshold. The investigated decay chain for this analysis was $`D^+D^+\pi ^0`$, $`D^+\overline{K}^0l^+\nu `$, and $`\overline{K}^0K^{}\pi ^+`$. The use of $`D^+`$ produced in $`D^+D^+\pi ^0`$ decays reduced the background. This analysis assumes that the thrust axis  of the event approximates the direction of the $`D`$ meson, as will be explained below. This approach does not work well for isotropic events. Thus, events whose ratio of the second and zeroth Fox-Wolfram moments  was less than or equal to $`0.2`$ were rejected. Specific ionization in the drift chamber and shower information in the electromagnetic calorimeter were used to select good electron candidates. Muon candidates were required to penetrate at least 5 interaction lengths in our muon detector. To select good $`\pi ^0`$’s, all $`\gamma \gamma `$ pairs for which $`|M_{\gamma \gamma }M_{\pi ^0}|`$ $`<`$ $`2.5\sigma `$ were accepted, where $`\sigma `$ is the standard deviation of the $`\pi ^0`$ mass measurement, obtained as a function of $`\pi ^0`$ momentum using clean $`\pi ^0`$ data samples. Other kinematic criteria chosen to optimize our sensitivity are given in Table I.

To reconstruct the momentum of the $`\nu `$, two methods were used to obtain up to three values of $`\stackrel{}{p}_\nu `$. In the first method, one or two values for the $`\nu `$ momentum were obtained assuming that the thrust direction of the event represents the $`D`$ direction. Given $`\stackrel{}{p}_{K\pi l}`$, the constraint that $`m_{K\pi l\nu }=m_{D^+}`$ provided an ellipsoid of allowed $`D`$ momenta. We then took the two intersections between the ellipsoid and the direction of the $`D`$. When there was no intersection, the point on the ellipsoid that lay closest to the $`D`$ direction was used. The values for the $`\nu `$ momentum were then the difference between these $`D`$ momenta and $`\stackrel{}{p}_{K\pi l}`$. The second method used the missing momentum of each event as an estimate of the $`\nu `$ momentum. If the missing momentum gave a value of $`m_{K\pi l\nu }`$ much greater than the mass of the $`D^+`$ (we used $`m_{D^+}`$ for the limit for convenience), then the $`\nu `$ momentum estimate from the second method was discarded. Among these three $`\nu `$ momentum estimates, the one that gave the value of $`\delta m=m_{K\pi l\nu \pi ^0}m_{K\pi l\nu }`$ closest to the known value of $`m_{D^+}m_{D^+}`$  was chosen. According to our simulation, this selection method gave the $`\nu `$ momentum closest to the generated $`\nu `$ momentum 75% of the time.

The previous CLEO analysis of this decay  superseded by this study did not use $`\nu `$ reconstruction. Instead it used the pseudo mass difference, $`\delta m_p=m_{K\pi l\pi ^0}m_{K\pi l}`$, and the resolution of $`\delta m_p`$ was much inferior to $`\delta m`$. The $`\nu `$ reconstruction reduced the statistical uncertainty in the $`(D^+\overline{K}^0l^+\nu _l)`$ measurement by 20% based on studies using signal and generic continuum Monte Carlo (MC) events. However, since our method for choosing among the three solutions for $`|\stackrel{}{p}_\nu |`$ tends to bias the $`\delta m`$ distribution of the background toward the signal region, it is important to understand the contribution of the background that peaks in the signal region. The uncertainty coming from this effect was included in the systematic errors.

Two quantities were fitted to obtain the number of signal events: $`m_{K\pi }`$ and $`\delta m`$. First, the data were divided into 50 bins (25 bins for the smaller muon data) over the $`\delta m`$ range from $`0.135`$ GeV/$`c^2`$ to $`0.235`$ GeV/$`c^2`$. For each $`\delta m`$ bin, the $`K\pi `$ mass distribution was plotted and the number of $`K^{}`$’s was extracted by a fit. Second, these yields were plotted as a function of $`\delta m`$, and the number of signal events was obtained by another fit. In the $`K\pi `$ mass fit, the signal shape was described by a $`p`$-wave $`K^{}`$ Breit-Wigner distribution. The mass and width of $`K^{}`$ obtained by the fit were consistent with the Particle Data Group (PDG) values . An analytic threshold function

$$(m_{K\pi }0.64)^\alpha \times e^{\beta (m_{K\pi }0.64)+\gamma (m_{K\pi }0.64)^2}$$
(1)
parameterized the shape of the background. The appropriateness of this function was tested using generic continuum MC events after discarding the signal among those events. The parameters $`\alpha `$, $`\beta `$ and $`\gamma `$ in this background function were unconstrained in the data fits. A typical $`K\pi `$ mass fit in the $`\delta m`$ signal region is shown in Fig. 1.

The signal shapes for the $`\delta m`$ fits were obtained from signal MC samples, and the background shape was described by a continuous function which was designed to accommodate the excess at $`m_{D^+}m_{D^+}`$ due to our neutrino momentum selection method. The shapes of these functions in the $`\delta m`$ fits for data were determined using a generic continuum MC sample. The $`\delta m`$ fit of the electron data is shown in Fig. 2.

The decay chain, $`D^+D^+\pi ^0`$, $`D^+K^{}\pi ^+\pi ^+`$, was used for the normalization mode. One advantage of this mode is that many systematic errors in its “measurement” are common with the $`D^+\overline{K}^0l^+\nu `$ mode and cancel in their ratio. Furthermore, the branching fraction of $`D^+K^{}\pi ^+\pi ^+`$ is well measured: $`(9.0\pm 0.6)\%`$ . Most of the kinematic requirements were the same as those in the semileptonic analysis, but the requirement on $`|\stackrel{}{p}_D^{}|`$ was tuned to minimize our sensitivity to the uncertainty in the fragmentation of the charm quarks.

For the normalization mode analysis, two successive fits of $`m_{D^+}`$ and $`\delta m`$ were used. The data were divided into 50 $`\delta m`$ bins from $`0.135`$ GeV/$`c^2`$ to $`0.185`$ GeV/$`c^2`$. For each $`\delta m`$ bin, a $`D^+`$ mass fit was used to extract the number of $`D^+`$’s. These yields of $`D^+`$’s were plotted as a function of $`\delta m`$. Then, this $`\delta m`$ plot was fitted to obtain the number of the normalization mode events.

Table II lists the important contributions to the systematic error. The dominant sources are the $`\delta m`$ signal shape estimation due to finite MC statistics, potential flaws in the estimate of the background shape obtained from the simulation of $`e^+e^{}q\overline{q}`$ and background arising from $`D^+\overline{K}^0\pi ^0l^+\nu _l`$.

The systematic error due to a potential difference in the background between the MC and data was estimated using a $`\pi ^0`$ mass sideband, $`3.5\sigma `$ $`<`$ $`|M_{\gamma \gamma }M_{\pi ^0}|`$ $`<`$ $`15.0\sigma `$. The $`\pi ^0`$ mass sideband data were analyzed as $`\pi ^0`$ signal band data and the difference in yields of data and generic continuum MC sample was used in the estimation of the systematic error due to the deficiency of the background MC in reproducing the data.

The upper limit of the branching fraction for $`D^+\overline{K}^0\pi ^0\mu ^+\nu _\mu `$ with respect to $`D^+K^{}\pi ^+\mu ^+\nu _\mu `$ is 0.042  at the $`90\%`$ confidence level. We used the product of this upper limit and the ratio of “efficiencies” that $`K^{}\pi ^0l\nu `$ and $`K^{}l\nu `$ are detected as signal in our analysis as the systematic error due to this source.

The systematic errors from each source were estimated individually for the electron and muon modes and then combined to obtain the overall lepton mode systematic errors.

The ratio $`R_l^+`$ was calculated using

$$R_l^+=\frac{N_{sl}}{ϵ_{sl}(\overline{K}^0K^{}\pi ^+)}\frac{ϵ_{had}}{N_{had}},$$
(2)
where $`N_{sl}`$ and $`ϵ_{sl}`$ are the observed number of events and the efficiency for $`D^+\overline{K}^0l^+\nu _l`$, respectively. Similarly, $`N_{had}`$ and $`ϵ_{had}`$ are the corresponding quantities for $`D^+K^{}\pi ^+\pi ^+`$. The measured branching ratios $`R_e^+`$, $`R_\mu ^+`$ and $`R_l^+`$, including their errors, are summarized in Table III. The values of $`R_e^+`$ and $`R_\mu ^+`$ agree well. The PDG value of $`(D^+K^{}\pi ^+\pi ^+)=(9.0\pm 0.6)\%`$  was used to calculate the branching fractions of $`D^+\overline{K}^0e^+\nu _e`$ ($`_e`$), $`D^+\overline{K}^0\mu ^+\nu _\mu `$ ($`_\mu `$), and $`D^+\overline{K}^0l^+\nu _l`$ ($`_l`$).

Table IV compares the results of this work with previous measurements  and the PDG averages . Note that the E691 measurement of $`R_e^+`$ and ours are the two most significant measurements, and they differ by about three standard deviations. We studied whether this difference (or part thereof) may arise from differences between the $`D^+\overline{K}^0e^+\nu _e`$ model used here (ISGW2)  and the one used in the E691 analysis (WSB) . Since the lepton momentum distribution expected from this decay has a significant effect on its detection efficiency in both experiments, we compared the electron energy and $`q^2`$ distributions from the two models, but found no significant differences. We did not identify any other explanations for the discrepancy.

Table V shows the decay rate for $`D^+\overline{K}^0e^+\nu _e`$ measured in this work and the theoretical predictions . Most predictions are consistent with both our result and that of E691.

Using our branching ratio measurements and the ratios between the form factors measured by E791 , we calculate $`A_1(0)=0.69\pm 0.07`$, $`A_2(0)=0.48\pm 0.08`$, and $`V(0)=1.25\pm 0.15`$, where the errors are the quadratic sums of the statistical and systematic errors in this analysis and in the E791 measurement.

The ratio of the vector and pseudo-scalar decay widths, $`r_{vp}=\frac{\mathrm{\Gamma }(DK^{}e^+\nu _e)}{\mathrm{\Gamma }(DKe^+\nu _e)}`$, tests quark models. This ratio was predicted to be in the range $`0.9`$ to $`1.2`$ by early quark models  and lattice gauge calculations . The E691 result implied that $`r_{vp}=0.48\pm 0.10`$. The newer ISGW2 model  accommodated this measurement well, predicting $`r_{vp}=0.54`$. Our measurement of $`R_l^+`$, on the other hand, implies that $`r_{vp}=0.99\pm 0.06\pm 0.07\pm 0.06`$, where the first, second and third errors are due to our statistical and systematic errors and the uncertainty in $`(D^+\overline{K}^0e^+\nu _e)`$, respectively, using the PDG value of $`(D^+\overline{K}^0e^+\nu _e)`$ $`=`$ $`(6.8\pm 0.8)\%`$ . Our ratio is consistent with older models and disagrees with the ISGW2 prediction by four standard deviations.

We gratefully acknowledge the effort of the CESR staff in providing us with excellent luminosity and running conditions. This work was supported by the National Science Foundation, the U.S. Department of Energy, the Research Corporation, the Natural Sciences and Engineering Research Council of Canada, and the Texas Advanced Research Program.


format: text

text:
# The X-ray spectrum of Fe16+ revisited with a multi-ion model

## 1 INTRODUCTION

The soft X-ray emission-line spectrum of Ne-like iron (Fe16+, Fe XVII) is one of the most extensively investigated spectroscopic systems in astrophysics. Owing to the high abundance of Fe and to the closed electronic shell structure of the Ne-like ground configuration 2s22p6, Fe16+ forms over a wide range of temperatures ($`kT`$ = 100 - 1000 eV). Hence, many spectra of hot astrophysical objects exhibit prominent Fe16+ lines. A better understanding of the spectrum of Fe16+ is even more important in light of the multitude of high-resolution X-ray spectra obtained with the grating spectrometers on board the Chandra and XMM-Newton observatories. The brightest lines of Fe16+ in the X-rays arise from the 2p-3d and 2p-3s transitions, which appear roughly at 15 and 17 Å, respectively. The analysis of X-ray spectra emitted during solar flares has shown the potential of these lines to serve as diagnostic tools, mainly for the diagnostics of electron temperature (Raymond & Smith, 1986; Smith et al., 1985), electron density (Schmelz, Saba, & Strong 1992; Waljeski et al., 1994; Phillips et al., 1996), and elemental abundances (Waljeski et al., 1994).

However, the usefulness of the Fe16+ lines has been hampered by two persistent problems, which emerge from extensive comparisons between observations, laboratory measurements, and theoretical models. The first is that the strongest line at 15.01 Å arising from the 2p6 1S0 \- 2p53d 1P1 transition appears to be significantly weaker than its predicted relative intensity. In particular, the intensity ratio of that line to the 2p6 1S0 \- 2p53d 3D1 line at 15.26 Å, commonly labeled $`I_{3C}/I_{3D}`$ (Parkinson, 1973), is consistently lower than the ratio calculated for a wide range of plasma parameters. This has led to the suggestion that the high oscillator-strength line at 15.01 Å is quenched by resonant scattering (Rugge & McKenzie, 1985; Schmelz et al., 1992; Waljeski et al., 1994; Schmelz et al., 1997; Phillips et al., 1996, 1997; Bhatia & Kastner, 1999; Saba et al., 1999). The continuing discrepancies between observations and theory have prompted several laboratory studies in which the plasma elemental composition, and to a large extent also the plasma conditions, can be controlled. In particular, the EBIT (electron beam ion trap) device offers the opportunity to isolate nearly entirely the Fe16+ ions from neighboring charge states and to measure their atomic properties with high accuracy. The EBIT measurements of Brown et al. (1998) first supported the claim that the line at 15.01 Å is affected by resonance scattering although to a lesser extent than originally believed. However, as suggested by Behar, Cottam, & Kahn (2001a) and later demonstrated in an EBIT measurement (Brown et al., 2001), a satellite line of Fe15+ that coincides in wavelength and blends with the 15.26 Å line can explain most of the observed $`I_{3C}/I_{3D}`$ ratios without needing to invoke resonant scattering. This conclusion is also supported by measurements in the PLT tokamak (Beiersdorfer et al., 2001).

The second problem with the Fe16+ lines is the overall intensities of the 2p-3s lines, which often appear enhanced relative to theoretical predictions, when compared with the intensities of the 2p-3d lines. This has led to suggestions that 2p inner-shell ionization of Na-like Fe15+ (ground configuration 2p63s) might contribute to the population of the 2p53s levels (Feldman, 1995). Other explanations suggested possible contributions to these level populations from dielectronic recombination (DR) of F-like Fe17+ (2p5) followed by radiative cascades (Saba et al., 1999; Liedahl, 2000; Laming et al., 2001). An EBIT measurement that focused on the overall intensity ratio of the 2p-3s lines relative to the 2p-3d lines (Laming et al., 2001) yielded a ratio that generally agrees with theoretical predictions, but is significantly lower than most astrophysical measurements. This further suggested that ionization and recombination processes, which are by and large absent from EBIT experiments and are also neglected in the theoretical models, play an important role in producing the Fe16+ lines in cosmic sources.

The recent high-resolution observations of stellar coronae with Chandra (Brinkman et al., 2000; Canizares et al., 2000; Behar et al., 2001a; Ayres et al., 2001; Drake et al., 2001; Huenemoerder, Canizares, & Schulz, 2001; Ness et al., 2001; Mewe et al., 2001) and XMM-Newton (Brinkman et al., 2001; Güdel et al., 2001a, b; Audard, Güdel, & Mewe, 2001a; Audard et al., 2001b) produced a variety of 2p-3s / 2p-3d intensity ratios. However, none of the observed ratios are fully consistent with existing theoretical models, further demonstrating that the above problems occur for a wide range of coronal sources. Additionally, thanks to the high dispersion and efficiency of the Reflection Grating Spectrometer (RGS) on board XMM-Newton, Fe16+ line spectra can now be obtained not only for stellar coronae (point sources), but also for extended sources. In observations of supernova remnants (Behar et al., 2001b) and elliptical galaxies (Xu et al., 2002), the relative Fe16+ line intensities show the same discrepancies with available models. In fact, for the analysis of the X-ray spectrum of NGC 4636, Xu et al. (2002) were required to invoke the empirical Capella spectrum into their model, because even the most updated theoretical models could not reproduce the observed emission. The pervasive discrepancies between the calculated spectra and those measured for a wide variety of hot astrophysical sources strongly suggest that the problems lie in the shortcomings of the atomic models for line formation rather than in the uncertainties associated with the astrophysical environments.

In an effort to improve the existing atomic models, which neglect atomic processes that involve neighboring ions, we construct a 3-ion collisional-radiative model that includes levels of F-like, Ne-like, and Na-like Fe. Such a model allows one to consider collisional excitations as well as the effect of recombination and ionization simultaneously. The main drawback of a naive, all-inclusive 3-ion model, in which rate equations for the three ionization states are solved simultaneously, is that the fractional ion abundances, which self-consistently result from the model, would be clearly erroneous. In order to obtain correct results for the ionization balance in a collisional-radiative model for a wide range of temperatures, one would have to include more than three ionization states, as well as numerous levels for each charge state (to adequately account for DR), which makes this approach impractical. Since we are interested in the line spectrum and not the ionization balance, we can circumvent this difficulty by adopting the approach formerly used in Doron et al. (1998, 2000), by which an independent ionization balance is imposed on the level-by-level, multi-ion model. The details of the current 3-ion method are given in the following section (§2). In §3, we give the results of the 3-ion model for the line-emission of Fe16+ and compare them with results obtained with a basic, single-ion model typical of those that have been used until now.

## 2 THEORETICAL METHOD

Collisionally ionized, astrophysical plasmas are characterized by relatively high temperatures and low densities. The most common atomic model used to describe the line emission from a particular ionic species (say Fe16+) under these conditions includes only electron impact excitations from the ground state and subsequent radiative decays, either directly to the ground state or via cascades. In the present work, we wish to extend this standard model and take into account also atomic processes involving neighboring ionization states. Specifically, we include dielectronic and radiative recombination from F-like Fe17+, inner-shell ionization from Na-like Fe15+, and resonant excitation through doubly-excited levels of Fe15+ (i.e., dielectronic capture followed by autoionization to excited levels). We assume the plasma is optically thin and collisionally ionized, i.e., photoexcitation and photoionization are unimportant. Assuming that the free electron energy distribution is Maxwellian, corresponding to an electron temperature $`T_e`$, the general set of rate equations for the density (population) $`n_j^{q+}`$ of an ion with charge $`q+`$ in a level $`j`$ in a steady state can be written as:

$`{\displaystyle \frac{d}{dt}}n_j^{q+}=n_e{\displaystyle \underset{kj}{}}n_k^{q+}\left[Q_{kj}(T_e)+Q_{kj}^{RE}(T_e)\right]+{\displaystyle \underset{k>j}{}}n_k^{q+}A_{kj}`$
$`n_j^{q+}\left\{n_e{\displaystyle \underset{kj}{}}\left[Q_{jk}(T_e)+Q_{jk}^{RE}(T_e)\right]+{\displaystyle \underset{k<j}{}}A_{jk}+{\displaystyle \underset{j^{\prime \prime }}{}}A_{jj^{\prime \prime }}^a\right\}+`$
$`n_e\left\{{\displaystyle \underset{j^{}}{}}n_j^{}^{(q1)+}S_{j^{}j}(T_e)+{\displaystyle \underset{j^{\prime \prime }}{}}n_{j^{\prime \prime }}^{(q+1)+}\left[\alpha _{j^{\prime \prime }j}^{DC}(T_e)+\alpha _{j^{\prime \prime }j}^{RR}(T_e)\right]\right\}`$
$`n_j^{q+}n_e\left\{{\displaystyle \underset{j^{\prime \prime }}{}}S_{jj^{\prime \prime }}(T_e)+{\displaystyle \underset{j^{}}{}}\left[\alpha _{jj^{}}^{DC}(T_e)+\alpha _{jj^{}}^{RR}(T_e)\right]\right\}=0`$ (1)
where $`j^{}`$ and $`j^{\prime \prime }`$ are levels pertaining to ions with charges $`(q1)+`$ and $`(q+1)+`$, respectively. The terms $`Q(T_e)`$ and $`Q^{RE}(T_e)`$, respectively, represent the rate coefficients for direct and resonant (RE) electron-impact excitation or de-excitation. The terms $`A`$ and $`A^a`$ denote the rate coefficients for spontaneous radiative decay and for autoionization, respectively. The terms $`S(T_e)`$ stand for the rate coefficients for collisional ionization, and $`\alpha ^{DC}(T_e)`$ and $`\alpha ^{RR}(T_e)`$ are the rate coefficients for dielectronic capture and for radiative recombination, respectively. The electron density is denoted by $`n_e`$. The entire set of equations that needs to be solved consists also of the appropriate rate equations for $`n_j^{}^{(q1)+}`$ and $`n_{j^{\prime \prime }}^{(q+1)+}`$, consistent with the atomic processes that appear in eq. (2). Throughout this work, we use $`q`$ = 16. Densities for ions other than the three considered are taken as zero. In order to obtain a unique solution, one of the rate equations in eq. (2) is substituted by a normalization condition, such as $`_j^{}n_j^{}^{(q1)+}+_jn_j^{q+}+_{j^{\prime \prime }}n_{j^{\prime \prime }}^{(q+1)+}=1`$, which leaves us with $`N1`$ independent equations for the $`N`$ levels included in the model.

A straightforward solution for the set of equations (2) necessarily yields the ionization balance through the results for the level populations $`n_j^{}^{(q1)+}`$, $`n_j^{q+}`$, and $`n_{j^{\prime \prime }}^{(q+1)+}`$. As argued above, this resulting ionization balance would clearly be wrong for most temperatures, as the set of equations consists of only three ionization states, and practically, an insufficient number of levels to determine the correct ionization balance. In order to avoid this error, the set of rate equations is solved by imposing the following constraints:

$`{\displaystyle \underset{j^{}}{}}n_j^{}^{(q1)+}/{\displaystyle \underset{j}{}}n_j^{q+}=f_{q1}/f_q`$ (2)
$`{\displaystyle \underset{j}{}}n_j^{q+}/{\displaystyle \underset{j^{\prime \prime }}{}}n_{j^{\prime \prime }}^{(q+1)+}=f_q/f_{q+1}`$ (3)
where $`f_{q1}`$, $`f_q`$, and $`f_{q+1}`$ represent the relevant ion fractions, and are calculated independently. In other words, we impose the ionization balance and only seek for the consequential population distribution among the individual levels. We are primarily interested in levels $`j`$ of the ion with charge $`q+`$. In the coronal approximation, where the ions are essentially in their ground state, and in the absence of metastable levels, it is sufficient to apply the constraints only through the ground levels of the three charge states. Moreover, collisional transitions among excited levels can be neglected. For Fe16+, this approximation holds for $`n_e10^{12}`$ cm-3. For the present study, we have adopted the fractional abundances calculated by Mazzotta et al. (1998). Table 1 presents the ion fractions of Fe15+ through Fe17+ for several chosen electron temperatures in the relevant range of 200 \- 1000 eV. The last column in the table gives the combined fraction of these three ions as a function of temperature. These total values are provided in order to give an idea of the collective importance of these ions in astrophysical plasmas at each temperature.

For Fe16+ we include in the model the configurations 2s22p6, 2s22p$`{}_{}{}^{5}nl`$ ($`n`$ = 3 to 7), 2s2p$`{}_{}{}^{6}nl`$ ($`n`$ = 3 to 7), 2s22p43$`l`$3$`l^{}`$, 2s2p53$`l`$3$`l^{}`$, 2s22p43$`l`$4d, and 2s22p43d4$`l`$. For Fe17+, the configurations 2s22p5 and 2s22p43$`l`$ are taken into account. Finally, the ground configuration 2s22p63s of Fe15+ is also included. All of these configurations add up to 2111 levels. In fact, the inclusion of the RE processes ($`Q^{RE}`$) in the set of rate equations (2) is equivalent to including the doubly-excited configurations of Fe15+. Since we are not interested specifically in the population of the Fe15+ levels, we have incorporated the RE rate coefficients from the literature (instead of including all of the appropriate doubly-excited levels of Fe15+), and by doing so, we significantly reduce the number of rate equations to be solved. The singly-excited configurations of Fe16+ up to $`n`$ = 7 have been shown to account for most of the radiative-cascades effects (Liedahl, 2000) and are taken into account in all of our models. The Fe16+ doubly-excited configurations that are included here explicitly produce about 50% of the total DR effect in the entire 200 - 800 eV temperature range. The contribution of higher DR resonances is obtained by means of extrapolation based on the total DR rates published by Chen (1988). In order to account for the radiative recombination contribution onto high electronic shells with $`n>`$ 7, we have employed a hydrogenic extrapolation method.

The solutions for the rate equations (2), which are the normalized level populations $`n_j^{q+}/_kn_k^{q+}`$, are used to obtain the line power $`I_{ji}^P`$, which is defined here by:

$$I_{ji}^P=\frac{n_j^{q+}}{\underset{k}{}n_k^{q+}}f_q\frac{A_{ji}}{n_e}$$
(4)
Subsequently, the line emissivity $`I_{ji}^E`$, which gives the total number of photons emitted in the plasma per unit time and per unit volume as a result of the radiative decay from level $`j`$ to level $`i`$, can then be found by:

$$I_{ji}^E=I_{ji}^Pn_eA_{Fe}n_H$$
(5)
where $`A_{Fe}`$ is the iron elemental abundance and $`n_H`$ is the hydrogen density.

The basic atomic quantities used in this work were generated by means of the multiconfiguration, relativistic HULLAC (Hebrew University Lawrence Livermore Atomic Code) computer package developed by Bar-Shalom, Klapisch, & Oreg (2001). Resonance excitation rates are not calculated explicitly, but incorporated from Chen & Reed (1989) for the 2p6 to 2p53$`l`$ ($`l`$ = s, p) excitations of Fe16+. Following those authors, resonant excitation to 2p53d levels is assumed to be unimportant due to their high energies and has been neglected.

## 3 RESULTS AND DISCUSSION

We have employed the model described above to calculate the Fe16+ level populations and line intensities for various temperatures in the range 200 - 1000 eV. In the following, we focus our presentation on the results for the seven most intense Fe16+ lines, which fall in the range 13.5 - 17.5 Å. The corresponding transitions, their labels, and measured wavelengths ($`\lambda `$) from Brown et al. (1998) are given in Table 2. Since in optically thin plasma, the relative intensities of these lines show negligible sensitivity to the electron density, the present results are valid for electron densities below 1012 cm-3.

### 3.1 Line Powers

The calculated line powers \[defined in §2 eq. (4)\] as a function of the electron temperature for the seven transitions listed in Table 2 are given in Figs. 1 and 2 by the solid curves. For comparison, we also give the line powers calculated with a single-ion model, i.e. a model that takes into account only direct electron impact excitations and subsequent radiative decays. The data points in the plots throughout this paper mark the explicitly calculated values, while the curves represent a spline interpolation between these points. It can be seen from this comparison that the 2p-3s lines (Fig. 1), as well as the weak 2p-3d line at 15.45 Å (Fig. 2), are significantly affected by processes involving neighboring ions over a wide range of temperatures, while the two strong 2p-3d lines at 15.01 and 15.26 Å, as well as the 2s-3p line at 13.83 Å, are much less affected by these processes. In particular, the intensity of the strongest 2p-3d line at 15.01 Å is almost identical in the two types of models, while the intensity of the 2p-3s line at 17.10 Å roughly doubles when contributions from neighboring ionization states are considered.

### 3.2 Contribution of the Various Mechanisms to Upper Level Population

In order to obtain an insight into the dramatic differences between the results of the 3-ion model and the single-ion model, we investigate the relative importance of each atomic process separately. This is done by running a sequence of 3-ion models, where in each run we turn off one of the atomic processes: resonant excitation (RE), dielectronic recombination (DR), radiative recombination (RR), or collisional ionization (CI). The resulting line powers in each of these models are presented in Table 3 for several electron temperatures in the range 200 - 800 eV. The results obtained with a traditional single-ion model, which includes direct collisional excitation (CE) alone, are also given for comparison. The first column in the table indicates the electron temperature and the second column describes the model. The following columns give the line power for each of the seven lines. In general, it can be seen from Table 3 that among the additional atomic processes, which are included in the 3-ion model, RE and DR are the most important.

At an electron temperature of 400 eV, the power of the 17.10 Å line, which is the most affected, receives contributions of 16%, 11%, 0.5%, and 3% from RE, DR, RR, and CI, respectively, which together enhance the line by 47% compared with the value obtained when calculating direct CE alone. At 600 eV, where the absolute Fe16+ line powers are still very high, the additional processes enhance the CE value for the 17.10 Å line by 57%, where RE, DR, RR, and CI, respectively, produce 10%, 21.5%, 0.5%, and 2.5% of the total line power. The other 2p-3s lines at 16.78 and 17.05 Å are also enhanced by the non-CE processes, albeit to a somewhat lesser degree (20% - 30%). The main reason for this is DR, which preferentially populates the upper level of the 17.10 Å line. Additionally, CI processes to that level are about twice as effective as CI to any other upper level of the 2$`l`$-3$`l^{}`$ transitions. The weak 2p-3d line at 15.45 Å is also enhanced, mostly by DR (RE processes towards the 2p53d levels are neglected). At 600 eV, 13.5% of the 15.45 Å line power is due to DR.

The trends of the various contributions with temperature are also interesting to note. The RE effect is predominant at low temperatures. It represents about a third of the 2p-3s line powers at 200 eV, but declines rather sharply with temperature. At 800 eV it contributes less than 7% to the 2p-3s line powers. The effect of DR, on the other hand, is small at low temperatures, but increases dramatically as the temperature increases. The enhancement of the 2p-3s lines due to DR reaches almost 30% at 800 eV. These two opposing effects are a direct consequence of the strong variations in the fractional ionic abundances with temperature and the fact that the RE effect is correlated with the Fe16+ abundance, while that of DR depends strictly on the Fe17+ abundance. As can be inferred from Table 1 and from Mazzotta et al. (1998), the abundance of Fe16+ dominates at temperatures up to $``$400 eV, while that of Fe17+ peaks at $``$600 eV. Interestingly, through most of the temperature range in which the Fe16+ line powers are relatively high, say 300 - 700 eV (see Figs. 1 and 2), both RE and DR are important, which means neither can be neglected when accurate line powers are needed. Recent unpublished calculations by Chen & Pradhan (2002) claim that Chen & Reed (1989) have underestimated the RE rate coefficients. Indeed, higher RE rates would imply even a more pronounced enhancement to the line powers than demonstrated here, especially in the low temperature regime, where RE is important. Unfortunately, the absence of full level-by-level RE rate coefficients in that paper precludes us from estimating this effect quantitatively. Notwithstanding, the DR and CI enhancements calculated in this work are by and large independent of the RE effect.

It is important to distinguish between the DR enhancement of Fe16+ lines calculated in this work and the effect of Fe15+ DR satellite lines, which was traditionally incorporated in plasma models (see e.g., Raymond & Smith, 1986). The DR effect, which is discussed here, produces the exact same parent lines of Fe16+ (i.e., 2p53$`l`$ \- 2p6), by means of DR processes that populate the upper levels of these transitions (e.g., 2p43$`l`$3$`l^{}`$ \- 2p53$`l`$). This should not be confused with DR satellites emitted by lower charge states at slightly longer wavelengths (e.g., 2p53$`l`$3$`l^{}`$ \- 2p63$`l`$ in Fe15+). The latter, in many cases, can be resolved from the parent lines by high-resolution spectrometers. Note, however, that the radiative decays (2p43$`l`$3$`l^{}`$ \- 2p53$`l^{}`$) that populate the upper levels of the Fe16+ lines, are in fact associated with DR satellites of Fe17+. This line emission is explicitly included in our calculations and is illustrated in the spectra shown below in §3.4.

### 3.3 Line Intensity Ratios

The large impact of the additional atomic processes considered in the 3-ion model on some of the line intensities directly affects the line ratios and, thus, the diagnostic applications of the Fe16+ system. In Figs. 3 \- 5, we present the most interesting line intensities normalized in most cases to the 2p-3d line at 15.01 Å, which is the strongest Fe16+ line at most temperatures and also the least affected by the additional processes (see Table 3). As in Figs. 1 and 2, the solid curves represent the results of the (complete) 3-ion model and the dotted curves represent the single-ion model. The most pronounced effect is found for the ratios involving the 2p-3s lines in Fig. 3. Since the spectral resolution of contemporary spectrometers, such as the Reflection Grating Spectrometer (RGS) aboard XMM-Newton and the Low Energy Transmission Grating Spectrometer (LETGS) aboard Chandra, is hardly sufficient to resolve the two 2p-3s lines at 17.10 and 17.05 Å, we also show the corresponding curve for the summed intensities of these two lines. The ratio involving the third 2p-3s line at 16.78 Å is given separately in Fig. 4 to avoid overcrowding the plot. The calculated ratios shown in Fig. 3 shed new light on the discrepancy as yet between the observed and calculated values for the ratio of the 17.10 Å and 17.05 Å lines to the line at 15.01 Å. Using the single-ion model, one obtains (in the relevant temperature range) a maximum value of 1.2 for the ratio (Iλ17.05 \+ Iλ17.10) / Iλ15.01, while the observed ratios from various coronal sources, including many observations of the Sun, are in the range of $``$ 1.5 - 2.6. A list of various measured ratios for these lines is presented in Table 4. The inclusion of RE, DR, and CI processes in the calculations by means of the 3-ion model yields values of $``$ 1.4 - 1.9 for this ratio at temperatures of 200 eV and higher, where Fe16+ forms predominantly. Even higher values, which are observed mostly in active regions on the Sun (e.g., Phillips et al., 1997; Saba et al., 1999) can be attained for this ratio at lower temperatures, which are realistic for these environments. In the special case of spatially resolved observations of the Sun, resonant scattering of the 15.01 Å line is also plausible. Alternatively, non-equilibrium ionization (NEI) conditions can also play a role in enhancing the 2p-3s / 2p-3d ratios. In particular, NEI conditions make the effects of ionization and recombination on line formation more pronounced than in equilibrium and would tend to enhance the 2p-3s lines even further. In conclusion, the major part of the observed ratios can be explained in terms of the 3-ion coronal-equilibrium model. For the abnormally high ratios observed mostly in the Sun, resonant scattering and NEI are also plausible.

The present results can shed more light on the intepretation of laboratory measurements. Laming et al. (2001) found for the 2p-3s / 2p-3d ratios good agreement between their EBIT measurements and simple, single-ion models. In the EBIT experiment, the selectivity of the ionic charge Fe16+ is supposed to be very high and the presence of Fe15+ and Fe17+ can therefore be minimized. In addition, Laming et al. (2001) used a relatively high-energy ($``$1 keV) electron beam. Under such conditions, the contributions of RE, DR, and CI are expected to be very small, and a single-ion CE model could therefore be adequate for reproducing the experimental results. However, in a very recent EBIT experiment, also with a high-energy beam ($`>`$0.8 keV), Beiersdorfer et al. (2002) measure 2p-3s / 2p-3d ratios that are almost a factor of 2 higher than those measured by Laming et al. (2001). These recent results could imply that even the basic CE models may be inadequate, or that high-lying resonances do contribute to RE even at these high energies, or a combination thereof. In any event, the effects of the additional processes, which are put forward in this work, are relevant to all plasmas in collisional equilibrium and as such, they ubiquitously enhance the CE (+ additional RE) rates measured with EBIT. Another indication for additional processes in plasma sources that are not present in EBIT experiments comes from the relative intensities of the 17.05 and 17.10 Å lines. In the EBIT measurements of Brown et al. (1998) the line at 17.05 Å is significantly stronger than the line at 17.10 Å, while in many astrophysical observations the two lines have comparable intensities (e.g., Rugge & McKenzie 1985; Huenemoerder et al. 2001; Behar et al., 2001a). This difference is consistent with our results as presented in Table 3, where DR enhances the 17.10 Å line much more than it does the 17.05 Å line.

As can be seen in Fig. 3, for $`kT_e`$ 600 eV, the newly calculated 2p-3s / 2p-3d ratios are more sensitive to the electron temperature than those calculated with the single-ion model, while at higher temperatures, the little sensitivity that existed in the single-ion model is, in fact, lost in the new model. In particular, the major role played by the RE processes at low temperatures significantly increases the temperature sensitivity in this regime. The (Iλ17.05 \+ Iλ17.10) / Iλ15.01 ratio decreases from 1.9 to 1.4 from 200 to 600 eV. The temperature sensitivity of these ratios diminishes rapidly thereafter, but in any case the absolute power of these lines drops sharply beyond 600 eV (See Figs. 1 and 2). Therefore, in multi-temperature sources often unresolved in astrophysical observations, the high-temperature contribution to the integrated spectrum of Fe16+ would be very small. In Fig. 4, one can see that the values obtained for the ratio Iλ15.01 / Iλ15.26 in the 3-ion model are not very different ($`<`$ 3%) than the single-ion model values in the entire temperature range. We note that our single-ion calculations for the ratio Iλ15.01 / Iλ15.26 yield a value of 3.4. The present model includes configurations up to $`n`$=7. The inclusion of even higher configurations would have only a negligible effect on the ratio. A value of 3.4 is about 12% higher than the value measured by Brown et al. (1998) in EBIT with a monoenergetic beam at 1.15 keV, and about 15% higher than the EBIT ratio measured by Laming et al. (2001) with a beam energy of 0.9 keV. The EBIT values are associated with errors of about 5%. Both measurements are just within the 15% accuracy range expected from the CE rate coefficients produced by the HULLAC code for this kind of atomic system.

The ratio of the weak 2p-3d line at 15.45 Å to the 2p-3d resonance line at 15.01 Å is also enhanced, albeit to a lesser extent than the 2p-3s lines, as can be seen in Fig. 5. The 2s-3p line at 13.83 Å is only very slightly affected by the additional processes (Fig. 5). The 15.45 Å line is enhanced for the most part by DR, while the effects of CI and RR are very small. The DR effect increases with temperature due to the increase in the abundance of Fe17+. The CI effect, on the other hand, depends both on the CI rate coefficients that increase with temperature, but also on the decreasing Fe15+ abundance. Since at low temperatures, the rate coefficients for inner-shell ionization (CI) are very small while the Fe15+ abundance is high and vice versa at high temperatures, the CI effect is supressed in ionization equilibrium conditions throughout the relevant temperature range. In NEI, however, the CI effect could be more pronounced. Note that following Chen & Reed (1989), RE processes to 2p53d levels are neglected in our model. However, according to Chen & Pradhan (2002), RE processes do, in fact, populate the upper level of the 15.45 Å line and enhance this ratio to values of 0.08 and higher. On the other hand, the measurements of Beiersdorfer et al. (2002) show a value of 0.04 for this ratio, which is totally consistent with the present results. Either way, the DR effect on this ratio demonstrated in Fig. 5 remains. Admittedly, the ratios in Fig. 5 could be rather hard to use for diagnostics due to their weakness in most sources. However, in bright sources, such as Capella for instance, the 13.83 Å line could be intense enough to provide useful temperature diagnostics (Behar et al., 2001a).

We point out that the present line ratios have been compared with those calculated by Gu (private communications) using a similar model based on the FAC atomic code (Gu, 2002). The overall 2p-3s / 2p-3d intensity ratios obtained with the two codes agree to within 1%. The agreement between ratios for individual lines is $``$5%.

### 3.4 Synthetic Spectra

As seen in §3.2 and §3.3, additional processes to direct collisional excitation (CE) have a major role in producing the soft X-ray lines of Fe16+. In order to provide a more visual illustration of the overall effect on the spectrum, total and partial synthetic spectra calculated for 400 eV are presented in Fig. 6. The upper trace in the figure gives the total spectrum, while the lower trace shows the explicit contribution of the CE processes alone. The middle trace represents the complementary contribution of non-CE processes. The different plots in Fig. 6 are all on the same scale. The spectra clearly show the importance of the non-CE processes to the 2p-3s lines at $``$ 17 Å. The 2p-3d lines ($``$ 15 Å) are almost unaffected by these mechanisms. The spectra also include lines of Fe17+ around 14.2 Å (2p-3d) and 16 Å (2p-3s), which are produced self consistently by the model and are rather weak at 400 eV. The line intensities for Fe17+, however, should be viewed with caution, since the model for Fe17+ is limited. Particularly, DR processes from Fe18+ are not taken into account. Additionally, the model does not include the DR satellite lines of Fe15+ and therefore these are absent from the plots in Fig. 6. At 400 eV, these DR satellites are expected to produce weak line emission mostly around 15.26 Å, as mentioned above in §1.

## 4 CONCLUSIONS

We have constructed a 3-ion model to account for ionization and recombination processes that produce the Fe16+ emission lines, in addition to the collisional and radiative processes among Fe16+ levels. It is shown that recombination (DR and RR), inner-shell ionization (CI), and resonant excitation (RE) processes, contribute significantly to the line powers and can resolve the existing discrepancies between observations and previous calculations within the framework of a coronal equilibrium model. Among the additional processes considered, RE and DR (of Fe17+) are found to play a major role in producing the Fe16+ lines. The present calculations reproduce the systematically high values observed for the 2p-3s / 2p-3d intensity ratios in many astrophysical sources. Only abnormally high ratios observed mostly on the Sun may still require one to invoke resonant scattering of the 15.01 Å line, or less likely transient conditions that depart from ionization equilibrium. Since non-equilibrium conditions are not required to explain most of the observed ratios, we defer the study of those cases to future work. The fact that recombination and ionization are important makes the present predictions dependent upon the ionization balance of Fe15+ \- Fe17+. We expect similar effects of inner-shell ionization and of recombination on the 2$`l`$ \- 3$`l^{}`$ lines to be present in more highly ionized systems (Fe17+ \- Fe21+) as well. However, a reliable assessment of this effect will require detailed models similar to the one employed here. The results of the present, elaborate, 3-ion model should also be important for analyzing the intensities of other lines that arise from the decay of the 2p53$`l`$ levels, i.e., 3$`l`$ \- 3$`l^{}`$ transitions. In the case of Fe16+, these transitions fall in the EUV range. In lighter Ne-like ions, these 3$`l`$ \- 3$`l^{}`$ transitions are clearly observed in the UV solar spectra obtained by the Solar Ultraviolet Measurements of Emitted Radiation (SUMER) spectrometer on board the SOHO satellite.

EB acknowledges ongoing work and many useful discussions with Steven Kahn on the complexity of the line intensity patterns of Fe16+ and their manifestation in X-ray observations. We are grateful to Ming Feng Gu for useful comparisons of the present results with some of his own calculations. We thank Daniel Savin for reading the manuscript carefully prior to its submission. RD was supported by a NASA Solar Physics Guest Investigator Grant S137816.


text:
# Demixing in mesoscopic boson-fermion clouds inside cylindrical harmonic traps: quantum phase diagram and role of temperature

## I Introduction

After the realization of Bose-Einstein condensation (BEC) in atomic gases , one of the most challenging endeavours in experiments on cold atoms is the cooling of fermionic isotopes of alkali atoms down to expected the superfluid transition. Theoretical estimates based on the BCS model  indicate that the temperature for the superfluid transition should be much lower than the Fermi temperature $`T_f`$, which is of the order of the BEC transition temperature. In fact, cooling fermions is harder than cooling bosons: the main difficulty arises from Fermi statistics, as $`s`$-wave collisions between spin-polarized fermions in a magnetic trap are forbidden by the Pauli principle. A common strategy uses sympathetic cooling, which is based on $`s`$-wave collisions between the fermions and a second gaseous component made either of fermions in a different internal state or of bosons. The latter choice seems to minimize the effects of Pauli blocking, which limit the process of cooling two fermionic components . Several experiments are currently in progress on trapping and cooling various boson-fermion mixtures, i.e. 6Li-7Li , 6Li-23Na  and 40K-87Rb . The lowest temperature attained in these experiments so far is about $`0.2T_f`$.

From the theoretical point of view a mixture of condensed bosons and fermions in the normal state is already an interesting system to study, because it can show spatial demixing of the two components . In the homogeneous gas at $`T=0`$ this is an example of a quantum phase transition, that is a phase transition induced by the interactions : on increasing the boson-fermion repulsion the system minimizes its total energy by placing the bosons and the fermions (or boson-fermion mixtures of different compositions) in different regions of space , even though this implies a high cost of kinetic energy at the interface. Boson-fermion demixing in a quasi-spherical trap has been studied by Nygaard and Mølmer , while the conditions for demixing inside a spherical trap have been set out by two of us  for $`N_b=N_f`$ and by Miyakawa et al.  for $`N_fN_b`$, with $`N_b`$ and $`N_f`$ the numbers of bosons and fermions in the trap. The conditions for observing demixing in the Paris experiment on a 6Li-7Li mixture in a cigar-shaped trap  has been analyzed by Akdeniz et al. . There are, of course, important consequences of finite system size on demixing in trapped mixtures, relative to the case of a homogeneous one. In particular, the transition is spread out in a finite system and the anisotropy of the trap acts differently on the two types of atoms (expressions such as phase transition and phase diagram will nevertheless continue to be used for brevity in the following). Locating the onset of partial demixing is also relevant to the practicalities of fermion cooling, since at that point the diminishing overlap between the two clouds will start reducing the effectiveness of the collisional transfer.

In this paper we give the general conditions under which phase separation of bosons and fermions occurs inside a harmonic trap at varying anisotropy and relative numbers of the two components. Within a mean-field model we make an extensive study of the equilibrium density profiles of the bosonic and fermionic component, ranging from the fully mixed state at small values of the scattering lengths to the partially demixed state at intermediate values of the coupling constants and finally to the regime of full phase separation. We show that the results can be understood on the basis of scaling laws and simply summarized into a phase diagram, which is expressed in terms of scaling parameters of the system at zero temperature. In the phase-separated regime we find several configurations having different symmetry and topology, some of which are metastable, and investigate the role of the anisotropy of the confinement in determining the minimum-energy configurations. We then extend our analysis to finite temperature, finding that on increasing temperature some phase-separated configurations turn or decay into others of higher symmetry before phase separation disappears. Finally, we show that phase separation between the fermions and a bosonic thermal cloud is also possible in principle on further increase of the boson-fermion scattering length.

The paper is organized as follows. In Sec. II we describe the model that we have used and give the limits of its applicability, together with a schematic description of the numerical method employed to find the phase-separated configurations. Sec. III summarizes the conditions under which phase separation occurs at $`T=0`$, illustrates various configurations which can be found in the phase-separated regime and gives the phase diagram for the lowest energy configurations. Sec. IV illustrates the effect of temperature on the configurations obtained for both small and large boson-fermion scattering length. Finally, Sec. V gives a summary and some concluding remarks.

## II The method

We describe the boson-fermion mixture by means of the particle density profiles, which are $`n_c(𝐫)`$ for the condensate, $`n_{nc}(𝐫)`$ for the bosonic thermal cloud and $`n_f(𝐫)`$ for the fermions. The components are subject to axially symmetric confining potentials given by

$$V_{b,f}^{ext}(𝐫)=m_{b,f}\omega _{b,f}^2(r_{}^2+\lambda _{b,f}^2z^2)/2$$
(1)
where $`m_{b,f}`$ are the atomic masses, $`\omega _{b,f}`$ the trap frequencies and $`\lambda _{b,f}`$ the trap anisotropies. We evaluate the density profiles in a mean-field model using the Thomas-Fermi approximation for the condensate and the Hartree-Fock approximation for the other clouds .

The Thomas-Fermi approximation assumes that the number of condensed bosons is large enough that the kinetic energy term in the Gross-Pitaevskii equation may be neglected . It yields

$$n_c(𝐫)=[\mu _bV_b^{ext}(𝐫)2gn_{nc}(𝐫)fn_f(𝐫)]/g$$
(2)
for positive values of the function in the square brackets and zero otherwise. Here, $`\mu _b`$ is the chemical potential of the bosons and the coupling constants are $`g=4\pi \mathrm{}^2a_{bb}/m_b`$ and $`f=2\pi \mathrm{}^2a_{bf}/m_r`$, with $`a_{bb}`$ and $`a_{bf}`$ the boson-boson and boson-fermion scattering lengths and $`m_r=m_bm_f/(m_b+m_f)`$ the reduced mass. The Hartree-Fock approximation, on the other hand, treats the thermal boson cloud and the fermion cloud as ideal gases subject to effective potentials $`V_{b,f}^{eff}`$, that is

$$n_{nc,f}(𝐫)=\frac{d^3p}{(2\pi \mathrm{})^3}\left\{\mathrm{exp}\left[\left(\frac{p^2}{2m_{b,f}}+V_{b,f}^{eff}(𝐫)\mu _{b,f}\right)/k_BT\right]1\right\}^1,$$
(3)
where

$$V_b^{eff}(𝐫)=V_b^{ext}(𝐫)+2gn_c(𝐫)+2gn_{nc}(𝐫)+fn_f(𝐫)$$
(4)
and

$$V_f^{eff}(𝐫)=V_f^{ext}(𝐫)+fn_c(𝐫)+fn_{nc}(𝐫).$$
(5)
The chemical potentials $`\mu _{b,f}`$ are determined by requiring that the volume integrals of the densities $`n_c(𝐫)+n_{nc}(𝐫)`$ and $`n_f(𝐫)`$ should be equal to the numbers $`N_b`$ and $`N_f`$ of particles.

The model is valid for the bosons when the diluteness condition $`n_ca_{bb}^31`$ holds and the temperature of the mixture is sufficiently below the condensation temperature. The fermionic component has been taken as dilute spin-polarized gas, for which the fermion-fermion $`s`$-wave scattering processes are inhibited by the Pauli principle and $`p`$-wave scattering is negligible . The condition $`k_fa_{bf}1`$ with $`k_f`$ the Fermi wave number should hold in the mixed regime, but this is not a constraint in the regime of phase separation where the boson-fermion overlap is rapidly dropping.

The density profiles at fixed numbers $`N_b`$ and $`N_f`$ are determined numerically by a self-consistent solution of the above equations. The profiles are obtained in two self-consistency loops: first, Eqs. (2) and (3) are solved at fixed chemical potentials starting from some initial guesses for the densities; then the chemical potentials are found iteratively by a standard algorithm for multi-variable minimization. Whereas in the mixed state the density profiles at convergence do not depend on the initial choice, different configurations are found in the phase-separation regime for the same values of the parameters by varying the initial conditions. We have used this fact to search for several possible metastable configurations, which will be illustrated in Sec. III B below.

## III Phase diagram at zero temperature

As is evident from the preceding section, a boson-fermion mixture under confinement is characterized by a large number of parameters. The aim of this section is to analyze what is the effect of varying them independently and to provide a unified understanding of the results using scaling laws.

The interaction energy grows on increasing the values of the boson-boson and boson-fermion coupling constants. There exists a threshold beyond which the total energy is minimized by a configuration in which the two components are spatially separated. The transition to the phase-separated regime is smooth owing to the finite size of the confined gas, and three different states are indeed recognizable: mixed, partially demixed, and fully demixed. Depending on the value of the relative strength of the two coupling constants, it may happen that either the fermions or the bosons are pushed away from the center of the trap.

### A Partial demixing

The boson-fermion interaction energy is proportional to the overlap between bosonic and fermionic clouds and at $`T=0`$ we have

$$E_{int}=fd^3rn_c(𝐫)n_f(𝐫).$$
(6)
We locate the onset of partial demixing by looking at the point where the boson-fermion overlap starts to decrease with increasing $`a_{bf}/a_{bb}`$.

In the following we focus on the simplest case where only one scattering length is varied, which also seems realistic from an experimental point of view. Then the most efficient way to approach phase separation is by changing $`a_{bf}`$ (see Sec. III B below). The behavior of the interaction energy as a function of the ratio $`a_{bf}/a_{bb}`$ at fixed $`a_{bb}`$ is shown in Fig. 1.

The position of the maximum of the interaction energy at fixed $`a_{bb}`$ can be estimated analytically from the condition $`E_{int}/f=0`$, by using the approximate expressions $`n_{b,f}N_{b,f}/(4\pi R_{b,f}^3/3)`$ and the values for the cloud radii $`R_f`$ and $`R_b`$ as obtained in the absence of boson-fermion interactions, $`R_f=(48N_f/\lambda _f)^{1/6}a_f`$ and $`R_b=(15\lambda _bN_ba_{bb}/a_b)^{1/5}a_b/\lambda _b^{1/3}`$ with $`a_{b,f}=(\mathrm{}/m_{b,f}\omega _{b,f})^{1/2}`$. The expression for the position of the maximum is

$$\frac{a_{bf}}{a_{bb}}|_{max}=\left(c_1\frac{N_f^{1/2}}{N_b^{2/5}}+c_2\frac{N_b^{2/5}}{N_f^{1/3}}\right)^1$$
(7)
where

$$c_1=\frac{15^{3/5}}{48^{1/2}}\frac{\lambda _f^{1/2}}{\lambda _b^{2/5}}\frac{m_f^{3/2}\omega _f}{2m_rm_b^{1/2}\omega _b}\left(\frac{a_{bb}}{a_b}\right)^{3/5}$$
(8)
and

$$c_2=\frac{48^{1/3}}{15^{3/5}}\left(\frac{6}{\pi }\right)^{2/3}\frac{\lambda _b^{2/5}}{\lambda _f^{1/3}}\frac{m_b\omega _b}{2m_r\omega _f}\left(\frac{a_{bb}}{a_b}\right)^{2/5}.$$
(9)
We recognise in Eq. (7) a geometric combination of two scaling parameters: $`c_1N_f^{1/2}/N_b^{2/5}`$ is dominant in the case $`N_bN_f`$ while $`c_2N_b^{2/5}/N_f^{1/3}`$ was previously identified  in the regime $`N_fN_b`$. The predictions obtained from Eq. (7) are indicated in Fig. 1 by vertical arrows. There clearly is very good agreement between the analytical estimate and the numerical results.

### B Full demixing

The regime of full demixing is taken to be reached when the boson-fermion overlap becomes negligible. As is shown in Fig. 1, within the three-fluid model the transition point does not depend on the number of bosons in the trap. This can also be predicted from the results obtained by Viverit et al.  for the homogeneous mixture, by using the values of the densities taken at the center of the trap. Within the same approximation for the fermion density as in Sec. III A we obtain the condition for phase separation at $`T=0`$ as

$$\alpha k_fa_{bb}>\left(\frac{a_{bb}}{a_{bf}}\right)^2,$$
(10)
where $`k_f=(48\lambda _fN_f)^{1/6}/a_f`$ and $`\alpha =[3^{1/3}/(2\pi )^{2/3}]m_bm_f/(4m_r^2)`$. Equation (10) is valid for different numbers of bosons and fermions as well as for different atomic masses and trap frequencies, and agrees very well with the phase-separation criterion given earlier by two of us  for the case where $`m_f=m_b`$ and $`\omega _f=\omega _b`$. It confirms that the two relevant parameters for describing the transition to the fully demixed regime are $`a_{bf}/a_{bb}`$ and $`\alpha k_fa_{bb}`$. We also immediately find that the value of $`a_{bf}`$ at the transition point scales as $`N_f^{1/12}`$ at fixed $`a_{bb}`$: this result is illustrated in the inset in Fig. 1.

Equation (10) has been verified by the numerical solution of the three-fluid model for a variety of different sets of parameters. First of all we have studied the location of phase separation in parameter space by varying the scattering lengths over much wider ranges than allowed by the usual off-resonant values (e.g. in the 6Li-7Li mixture), as might be attained experimentally by exploiting optically or magnetically induced Feschbach resonances . For these calculations we have chosen the geometrical parameters of the traps and the numbers of atoms as in the Paris experiment on 6Li-7Li : namely $`N_bN_f=10^4`$, $`\omega _b/2\pi =4000`$ Hz, $`\omega _f/2\pi =3520`$ Hz and $`\lambda _b\lambda _f1/60`$. The results are summarized in a quantum phase diagram at $`T=0`$ in the plane $`\{a_{bf}/a_{bb},k_fa_{bb}\}`$, which is reported in Fig. 2. This contains all the lowest-energy configurations that we have found, with the inset showing an enlargement of the region near the origin of the above-mentioned plane. The curve drawn inside the inset in Fig. 2 is obtained from Eq. (10) and is in full agreement with the numerical results in separating the region of full demixing from the mixing or partial-demixing regions.

We turn to a detailed account of the configurations corresponding to the various symbols and letters in the phase diagram in Fig. 2. Below the phase separation line in the inset two different types of configurations are found, one with a core of bosons in partial overlap with an envelope of mainly fermions (triangles down) and a specular one with a core of fermions enveloped by mainly bosons (triangle up). An interesting configuration that we have observed in the latter case shows coexistence of a mixed state with a purely bosonic phase (see Fig. 3): this combination is not stable in the homogeneous mixture .

In the regime of phase separation a variety of different configurations are observed, including several that require a break of the symmetry imposed by the confining harmonic potentials. The letters a-d in the phase diagram in Fig. 2 give the locations of four such energetically stable structures, which are shown in Fig. 4 together with some other metastable structures. The symmetric configuration with boson inside and fermion outside, indicated by (a), is obtained in Fig. 4 with the choice $`a_{bb}=600`$ Bohr radii and $`a_{bf}=10a_{bb}`$, while the complementary configuration (b) formed from fermions inside and bosons outside is obtained with $`a_{bb}=64600`$ Bohr radii and $`a_{bf}=a_{bb}/2`$. Other configurations obtained with the same scattering lengths as the symmetric structure of type (a) are shown in Fig. 4: these are (c) a boson torus inside a fermion cloud, (d) a threefold symmetric structure formed from fermions surrounded by a shell of bosons inside a fermion envelope, (e) an asymmetric structure in which the bosons are shifted away from the center of the trap along the $`z`$ axis, and (f) a “sandwich” formed by bosons inside a fermion cloud. Finally the structure indicated by (g) in Fig. 4, which consists of a torus of fermions around an elongated boson core, is found with the choice $`a_{bb}=6250`$ Bohr radii and $`a_{bf}=2a_{bb}`$.

With the exception of (a) and (b) the configurations shown in Fig. 4 are very different from the density distributions that are met in the mixed phase, and thus could be identified in an experiment by looking at the column-density pictures of the atomic clouds. As an illustration we give in Fig. 5 the column-density images corresponding to the case (c) of a bosonic torus, and in Fig. 6 those for the case (f) of a boson sandwich.

Finally, we have examined the role of the anisotropy of the confinement on the configurations of lowest energy in the phase-separated regime. To this end we have evaluated the stability of the lowest-energy configuration by varying the anisotropy from $`\lambda _f=1/60`$ to $`\lambda _f=1`$ while keeping the ratio $`\lambda _b/\lambda _f1`$, with the choice $`a_{bb}=600`$ Bohr radii and $`a_{bf}=10a_{bb}`$ for the scattering lengths. The most stable configuration remains the symmetric structure (a) except that for $`\lambda _f1`$ the structure (d) has the lowest energy. For this particular choice of parameters the relative difference in energy between the metastable configurations does not vary much with $`\lambda _f`$.

## IV The role of temperature

At finite temperature an increasing fraction of the bosons populates the thermal cloud, thus depleting the condensate. In harmonic confinement the thermal cloud is wider and more dilute than the condensate, and thus is much less interacting with the fermions. As a result the thermal cloud is not yet in phase separation when the condensate already is.

At low temperatures the non-condensate fraction is small and its presence does not affect the description of the phase separation between the condensate and the fermions that we have given in the preceding section. In particular we have verified that at $`T`$ of the order of $`0.2T_f`$, as are reached in current experiments, all the exotic configurations (c) to (g) remain possible. They are not found at higher temperatures: we show as an example in Fig. 7 the behavior of the boson torus configuration (c) as a function of temperature. For $`0.2T_f<T0.3T_f`$ a part of the condensed bosons moves from the torus to the center of the trap and for $`T>0.3T_f`$ the phase separation is lost.

At still higher temperatures the thermal cloud starts to be consistently populated and we have found that phase separation between the fermions and the entire bosonic cloud becomes possible for very large values of $`a_{bf}/a_{bb}`$ . In this case, as is illustrated in Fig. 8, the fermion cloud exhibits a central hole having the size of the bosonic thermal cloud.

## V Summary and concluding remarks

We have performed an extensive analysis of the phenomenon of phase separation for a boson-fermion mixture under confinement on the basis of the atomic density profiles of the two components. First of all we have pointed out that three regimes can be identified owing to the finite size of the system: a mixed phase where the bosonic and fermionic clouds overlap, a partially demixed state where the overlap is decreasing, and a phase-separated regime where there is no overlap between the two components at zero temperature. We have given analytical expressions for the positions of the two cross-overs between these regimes in terms of the physical parameters of the system, and verified these conditions by an extensive numerical study based on varying the coupling constants, the numbers of atoms, and the geometry of the trap. General expressions for the transition lines are particularly useful in view of the several experiments in progress , which are run with different experimental parameters such as the relative numbers of bosons and fermions, the atomic masses, and the trap anisotropies.

The transition to the partially phase-separated regime, that we have located where the boson-fermion interaction energy starts to decrease on increasing the boson-fermion scattering, is important for the experiments on fermion cooling: the boson-fermion collision rate is decreased by a geometrical factor which is proportional to the overlap between the two clouds. At fixed boson-boson scattering length we have given a simple model for the transition in terms of a combination of two scaling parameters, which characterize the cases $`N_bN_f`$ and $`N_fN_b`$.

The transition to full demixing has been illustrated by a universal phase diagram expressed in terms of scaling parameters of the system. It is remarkable that in our model the location of this transition does not depend on the number of bosons. In the region of phase separation we have found several configurations with various topology and identified the minimum-energy configurations. We have also investigated the effect of increasing temperature on phase separation. Some exotic configurations that are found in the phase diagram at $`T=0`$ turn into others of higher symmetry before the condensate-fermion demixing is lost.

The main limitations of the present model come from the two basic approximations that we have employed: the Thomas-Fermi approximation for the condensate and the mean-field Hartree-Fock approximation for the other gaseous clouds. The Thomas-Fermi approximation breaks down for a very small condensate or for negative values of the boson-boson scattering length, and one would need to solve the full Gross-Pitaevskii equation. A description beyond mean field is required when the dilution parameters increase and for a trapped gas this is still an open problem. The analogue of the Beliaev expansion at $`T=0`$ has recently been derived for the homogeneous boson-fermion mixture .

Of course, the evolution of the mixture towards phase separation modifies not only its equilibrium properties, but also its dynamical properties and collective excitation spectrum. The approach to the transition will be signalled by a softening of the frequencies of surface or bulk modes having the appropriate symmetry for a phase-separated configuration to be attained. Experiments measuring dynamical properties could therefore provide an alternative method for revealing the onset of demixing. Calculations on these spectra are in progress.

## ACKNOWLEDGMENTS

This work was supported by INFM through the PRA2001-Photonmatter.


text:
# Absorbers and Globular Cluster Formation in Powerful High Redshift Radio Galaxies

## 1 Introduction

Collimated radio-luminous outflows are by now detected from the cosmological neighborhood up to redshifts in excess of five. According to the unified scheme these jet-sources come in two classes: Radio-loud quasars, where the current redshift record is $`z=5.8`$ (Fan et al. (2001)), and radio galaxies. Among the latter ones the highest known redshift is $`z=5.2`$ (van Breugel et al. (1999)). Infrared observations indicate that radio galaxies contain a luminous quasar hidden by a dusty torus, at least at redshift above $`z0.8`$ (Meisenheimer et al. (2001)). The maximum power of radio galaxies increases with redshift (De Breuck et al. (2000)), the suspicion being that the most powerful active galaxies at $`z>2`$ develop into BCGs (Carilli et al. (2000), and references therein). This idea is supported e.g. by the fact that the co-moving space density of clusters at low redshift is comparable to that of powerful AGN at $`z=2.5`$ (West (1994)).

There is good reason to believe that the jets of higher redshift sources bore into a medium of higher density. If one just considers the cosmological expansion, the number density of the intergalactic material (IGM) scales like $`n_{\mathrm{IGM}}(1+z)^3`$, on large scales. This means e.g. that if the well known radio galaxy Cygnus A, which has a surrounding IGM number density of a few times $`10^2\mathrm{cm}^3`$ (Carilli & Barthel (1996)), was located at redshift $`z=3`$, its jet had to bore into a medium with $`n1\mathrm{cm}^3`$. However, the environments of the most powerful radio galaxies, which are believed to highlight the highest density peaks in the early universe, are not a priori expected to simply follow the Hubble flow. For example, it has to be taken into account that the associated galaxies appear as bright ellipticals only at redshift $`z<1.5`$, whereas above this redshift one finds clumpy regions, distributed over upto 100 kpc (Pentericci et al. (2001)), and a considerable fraction of the gas is transformed into stars around that redshift. Direct estimates from emission line observations indicate IGM densities in the range of $`0.1\mathrm{cm}^3`$ (van Ojik et al. (1997)) upto $`10\mathrm{cm}^3`$, if one ascribes the dominant energy source for ionization to the bow shock of the jet (Bicknell et al. (2000)). If the jets themselves consist in electron proton plasma and are subrelativistic on large scales – e.g. in Cygnus A a kpc scale jet velocity ($`v_\mathrm{J}`$) of 0.4 $`c`$ is favoured (Carilli & Barthel (1996), $`c`$ being the velocity of light) – the number density in the jet ($`n_\mathrm{J}`$) is given by

$$n_\mathrm{J}=\frac{\dot{M}_\mathrm{J}}{\pi R_\mathrm{J}^2v_\mathrm{J}m_\mathrm{H}}=1.3\times 10^4\mathrm{cm}^3\frac{\dot{M}_\mathrm{J}c}{3v_\mathrm{J}}\left(\frac{\mathrm{kpc}}{R_\mathrm{J}}\right)^2.$$
(1)
Here $`\dot{M}_\mathrm{J}`$ is the jet’s mass flux rate (in solar masses per year), $`R_\mathrm{J}`$ is the jet radius (about 0.6 kpc in Cygnus A, Carilli & Barthel (1996)), and $`m_\mathrm{H}`$ the mass of hydrogen. Adopting an IGM number density of $`1\mathrm{cm}^3`$ for high redshift radio sources gives a density contrast $`\eta =n_\mathrm{J}/n_{\mathrm{IGM}}10^4`$. Due to the lower environment density (Daly (1995)), $`\eta `$ is thought to be higher in low redshift sources. Because the expansion velocity of the radio structure scales proportional to $`\sqrt{\eta }`$, the average extension of radio galaxies with $`z<1.5`$ of 100 kpc versus 10 kpc at $`z>2`$ (Carilli et al. (2000)) could be nicely explained by such a cosmological IGM density increase, if the average active time stays constant. Another indication for high density environments at high redshift is the observed bending of the kpc scale radio jets (e.g. Pentericci et al. (1997)). This underlines the clumpy nature of the IGM around high redshift radio galaxies (HZRG).

The optical emission of radio galaxies is often aligned with the jet axis. While below $`z=1.5`$ this emission is highly polarized, indicating mainly scattered light from a hidden AGN, above $`z=2`$ less then 2% polarization is found and stellar absorption lines are detected (Carilli et al. (2000)). This could indicate that – at high redshift – jets induce star formation. The aligned light could also reflect the large scale matter distribution (West (1994)).

Radio galaxies are surrounded by a Ly-$`\alpha `$ halo. The size of this halo grows with redshift (Carilli et al. (2000)). Whilst sources at $`z<1.5`$ have halos of $`10^8\mathrm{M}_{}`$ and an extension of 10 kpc, $`z>2`$ sources have 100 kpc halos of $`10^9\mathrm{M}_{}`$. The width of the emission lines is typically $`(1000\pm 500)`$ km/s (van Ojik et al. (1997); De Breuck et al. (2000)). HZRG with diameters less than 50 kpc show associated narrow absorption lines with a width of $`(40\pm 30)`$ km/s (van Ojik et al. (1997)). The absorption is often saturated with column densities of $`(10^{18}10^{20})\mathrm{cm}^2`$, and is preferentially blue-shifted by upto 250 km/s. These absorption systems were thought of to consist of $`10^{12}`$ solar system sized clouds by van Ojik et al. (1997). Advanced modeling with photoionization codes lead to the proposition that the absorber was a low density region, surrounding every HZRG (Binette et al. (2000)). None of these models take into account the hydrodynamic facet of the problem. This letter is therefore dedicated to an understanding of the above mentioned features based on hydrodynamic modeling of HZRGs.

## 2 A hydrodynamic model for high redshift radio galaxies

One of the basic differences between jets at low and high redshift is the environmental density. We consider a – typical – jet with a density contrast of $`\eta =10^4`$, a velocity of $`v_\mathrm{J}=10^{10}`$ cm/s (internal Mach number 85), a jet radius of 1 kpc, and a density of the homogeneous external medium of 1 $`\mathrm{cm}^3`$. This results in a kinetic jet luminosity of $`10^{45.4}`$ erg/s, which is transfered by some efficiency factor to radio luminosity. The velocity of advancement of the jet’s head – and the bow shock in front of it – is expected to be $`v_\mathrm{B}\sqrt{\eta }v_\mathrm{J}=1000`$ km/s. This velocity causes the external gas to get heated to a temperature of $`T=1.38\times 10^7\mathrm{K}(v_\mathrm{B}/1000\mathrm{km}/\mathrm{s})^2`$, according to standard hydrodynamic shock jump conditions, neglecting the internal energy of the pre-shock gas. The cooling time of this post-shock gas due to thermal bremsstrahlung with an emissivity of $`ϵ_{\mathrm{brems}}=2.1\times 10^{27}(n/\mathrm{cm}^3)^2\sqrt{T/\mathrm{K}}`$ $`\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1`$ will be:

$$t_{\mathrm{cool}}\frac{2nkT}{ϵ_{\mathrm{brems}}}3\times 10^6\frac{\sqrt{T/10^7\mathrm{K}}}{n/4\mathrm{cm}^3}\mathrm{yrs},$$
(2)
which is less than the typical lifetime of a jet. Therefore, cooling has a major effect on the evolution of the post-shock gas. This gas will cool down to high densities and a temperature of about $`10^4`$ K, where the cooling curve drops significantly (Sutherland & Dopita (1993)).

We have solved the hydrodynamic equations for the above mentioned parameters numerically with the code NIRVANA\_C (Ziegler & Yorke (1997)), upgraded due to optically thin cooling (Thiele (2000)). NIRVANA\_C is a second order accurate finite difference scheme. Cooling comes into play via the equation for the internal energy density ($``$):

$$\frac{}{t}+(v)=p(v)\mathrm{\Lambda },$$
(3)
where $`\mathrm{\Lambda }`$ is the global cooling function. NIRVANA\_C also solves the ionization and recombination equations, by a time implicit method, in order to allow for nonequilibrium processes. For the present simulation, three species were used, namely ionized and neutral hydrogen, and electrons. Below $`T=10^6`$ K, an average cooling curve for solar metalicity was applied (Stone & Norman (1993)). The simulation was carried out in axisymmetry (2.5D), and the jet is resolved with 10 points. With that resolution, NIRVANA\_C marginally resolves shocks in the jet beam and gets basic flow parameters right. Turbulence on small scales as well as details of the beam, especially of the head region, are unresolved (Krause & Camenzind (2001)). The simulation was stopped when the thickness of the shell behind the bow shock shrinked below the resolution, which was clearly the case after a simulation time of $`6.0\times 10^6`$ yrs.

The results are quite unusual compared to simulations of higher $`\eta `$: Most of the mass accumulates into an almost spherical shell behind the bow shock. The maximum density there is $`47\mathrm{cm}^3`$, and the gas has cooled to an average temperature of $`1.7\times 10^4`$ K in regions with density above $`10\mathrm{cm}^3`$, which still cover almost the whole shell surface. The mass accumulated in this shell is $`6.19\times 10^{10}\mathrm{M}_{}`$, 75 % of which are neutral. This gives an average neutral column for the shell of $`N_{\mathrm{HI}}=3.8\times 10^{21}\mathrm{cm}^2`$, which has to be compared to the amount of matter displaced by the shell: $`6.32\times 10^{10}\mathrm{M}_{}`$. It follows, that 98 % of the mass are accumulated in the thin shell. Therefore 2 % of the mass or $`1.4\times 10^9\mathrm{M}_{}`$ were entrained into the cocoon.

The cocoon itself has a very peculiar structure. It is well known, that jet cocoons broaden when lowering $`\eta `$ (Norman et al. (1983)). But in this simulation, the cocoon is not elongated along the jet beam, forming individual vortices. Instead, all the vortices join into one big vortex with a diameter of roughly 8 kpc. The vortex shedding is sometimes violent: while, at the time shown, only a small eddy is ejected at the beam head, at other times we observed upto half the jet beam splitting and forming an eddy. The big vortex turns around at a maximum velocity of $`0.20.25v_\mathrm{J}`$.

The velocity of the jet head differs considerably from higher $`\eta `$ simulations. At early times, it is about 10 times higher than the expectation above. Only after about two million years, the jet head moves, relatively constant, at maximum expected pace.

Within the cocoon, the sound speed is, on average, 20000 km/s. Therefore, pressure differences within the cocoon gas are rapidly communicated, with the result that the pressure is almost constant. The simulation shows that the bow shock is, at early times, perfectly spherical. It looks more like a supernova bubble than a jet cocoon. Indeed, if one compares the velocity of a spherical blast wave (e.g. Wünsch & Palous (2001), with constant energy ejection): $`v_{\mathrm{BW}}=1412\left(\frac{L/10^{46}\mathrm{erg}/\mathrm{s}}{n/\mathrm{cm}^3}\right)^{1/5}\left(\frac{t}{10^6\mathrm{yrs}}\right)^{2/5}\mathrm{km}/\mathrm{s}`$ to the velocity of the jet head, one finds, that the jet head outruns the blast wave only at times above $`1.9\times 10^6`$ yrs, for the parameters of our simulation. This is what we observe. We rewrite the above equation for the case, when the energy injection comes from a jet: $`v_{\mathrm{BW}}/v_\mathrm{J}=0.6(R_\mathrm{J}/R_{\mathrm{BW}})^{2/3}\eta ^{1/3}`$. The blast wave velocity depends on $`\eta ^{1/3}`$, whereas the jet head velocity depends on $`\eta ^{1/2}`$ (light jet limit). Because of that, the spherical state of the bow shock can be observed only at the earliest times, for higher $`\eta `$.

## 3 Comparison to observations

Because of the large neutral column, this shell is a good candidate for the absorber in HZRGs. The width of an absorption line caused by this absorber would be approximately the sound speed within it (Dyson et al. (1980)), which is in our simulation about 55 km/s, on average. This is also in remarkable agreement with observations (van Ojik et al. (1997)). Typically, the material that is just heated by the shell emits Lyman $`\alpha `$ at a few percent of the kinetic jet luminosity. This is less then the observed Ly $`\alpha `$ luminosity, indicating that the bow shock is not the main emitter. At late times compared to $`t_{\mathrm{cool}}`$, it is justified to approximate the bow shock with an isothermal shock description. Following Steffen et al. (1997), we derive for the thickness of the spherical thin shell:

$$d=27.7\frac{(\overline{m}/m_\mathrm{p})(r/10\mathrm{kpc})(T/10^4\mathrm{K})}{(v_{\mathrm{bw}}/100\mathrm{km}/\mathrm{s})^2}\mathrm{pc}.$$
(4)
Here, $`\overline{m}/m_\mathrm{p}`$ denotes the mean molecular weight per proton mass. The sound crossing time through this shell is about $`5\times 10^5`$ yrs, which is much longer than the cooling time (about $`10^4`$ yrs, Sutherland & Dopita (1993)). Therefore, it is possible, that the shell breaks up into individual subshells (Dyson et al. (1980)). This could well correspond to the multiple associated absorption systems, often observed in HZRGs. These absorbers are preferentially blue shifted, a natural feature of our model. But almost never, a velocity above 250 km/s is observed. Approximating the sideways expansion of the shell by the laws for the spherical shell also for some time after the jet head begins to shape the shell, one can calculate $`\eta `$, given the bubble size ($`>12`$ kpc) at observation time. This yields $`\eta =2\times 10^5`$, which means that the full range of currently known observations can be explained if the environmental density is $`5\mathrm{cm}^3`$. $`\eta `$ in the above formula could be replaced by $`\eta _0(R/R_\mathrm{J})^p`$ to account for a declining density distribution. The $`\eta _0`$ would be lower, in that case. It is interesting to note here, that very similar numbers have been found from quite different arguments for the environment of the radio galaxy 4C 41.17 (Bicknell et al. (2000)).

The bubble interior looses a considerable amount of internal energy via synchrotron radiation, which is not included in the simulation. This lowers the requirement for the external density. Although the observed hydrogen mass in these systems barly exceeds $`10^9\mathrm{M}_{}`$, it is likely that they contain significantly more mass. Given that HZRG host the most powerful quasars of the universe, we should expect black hole masses of the order $`10^9\mathrm{M}_{}`$. The gas mass should be more than that. If the gas mass would be of the order $`10^{1112}\mathrm{M}_{}`$, as is suggested here, the gas to black hole mass ratio would be similar to the bulge to black hole mass ratio found at low redshift. In that case, most of the shell’s mass would be – at any evolutionary state, the systems are observed at so far – in cold fragments just about to form stars.

The shell is not only thermally, but also gravitationally unstable. The time, when the instability on the shell surface occurs for the first time is (Wünsch & Palous (2001), again for constant energy injection rate and for the parameters of our model): $`t_\mathrm{i}=2.4\times 10^7\left(\frac{(c_\mathrm{s}/55\mathrm{km}/\mathrm{s})^5}{(r/\mathrm{kpc})^2(\eta /10^4)(3v_\mathrm{J}/c)^3(n/\mathrm{cm}^3)^5}\right)^{1/8}\mathrm{yrs}`$. Gravitational and thermal instability will support each other. The Jeans mass in the shell is $`10^610^7M_{}`$. Hence, the shell will form stars in globular clusters of that mass. Assuming an efficiency of $`10\%`$, the shell would form, roughly, $`10^4`$ globular cluster systems. This is in remarkable agreement to what is found in nearby BCGs, like e.g. M87 (Harris et al. (1998)): BCGs show an excess of globular cluster systems of about $`10^410^5`$ compared to non BCG ellipticals with comparable luminosity. Sometimes, two distinguished globular cluster populations are observed. Harris et al. (1998) consider various formation scenarios in detail and conclude that a kind of galactic wind must have driven out a large fraction of the galaxies gas, before it was able to form stars.

When the shell fragments, its covering factor decreases further and the absorption vanishes. This scenario could nicely explain, why no absorbers are observed for galaxies larger than 50 kpc. Afterwards the shell would become visible in the optical, due to the newly born stars, which will also increase the ionization of the remaining parts of the absorbing shell. A bubble shaped star forming region like that can be found around the radio galaxy 1243+036 (van Ojik et al. (1996)).

These associated absorption systems also seem to be observed in high redshift quasars (Baker (1997); Baker (2001)). There, they have velocities of a few thousand km/s, consistent with the model, if one assumes that the line of sight in quasars is not far from the jet axis. In addition, this part of the shell could be accelerated by photoabsorption of photons from the quasar (Falle et al. (1981)).

In the simulation, we find that $`10^9M_{}`$ of shocked IGM are entrained into the cocoon. This number should not be taken too literally, because small scale turbulence could change that and high resolution studies are required to get this number accurate. Nevertheless, it is slightly more than the typical observed Lyman $`\alpha `$ emitting gas mass. Higher resolution would also be essential in order to determine the mixing properties and density of this gas. The big vortex, we observe in the cocoon is an ideal accelerator for the entrained gas. In the simulation, the gas is swept along, and accumulates mainly on the left-hand side and along and next to the jet beam. There it is spun up to a velocity of $`1000`$ km/s. We expect, that at higher resolution, small scale Kelvin-Helmholtz instabilities will entrain mass at the boundary of the vortex. This is also observed in other simulations (Krause & Camenzind (2001)). Given the high velocity of this propeller, it seems quite likely that emission line gas could be accelerated to the observed velocity of $``$ 1000 km/s. All that should be studied in more detail with simulations of higher resolution.

###### Acknowledgements.

This work was supported by the Deutsche Forschungsgemeinschaft (Sonderforschungsbereich 437).