text:
# A code for optically thick and hot photoionized media

## 1 Introduction.

A number of codes were built in the past to compute the structure and the emission of photoionized media. With time, these codes became more and more sophisticated, and able to treat a larger number of situations. In the sixties they were specially designed for planetary nebulae and HII regions, i.e. for dilute and optically thin media (except in the Lyman continuum) ionized by the thermal radiation of a hot star, and were rapidly used also for ionization by a non thermal continua extending in the X-ray range, i.e. for the Narrow Line Region of Active Galactic Nuclei (AGN). At the same time another generation of codes was developed for optically thin hot collisionally ionized media, and used for the solar corona, for supernova remnants, and for the hot intergalactic gas. At the end of the seventies, observations of compact X-ray sources implying the presence of a Thomson thick medium incited Ross (1979) to develop a radiative transfer method for the continuum using a modified version of the Kompaneets equation. In the same paper he introduced the so-called “escape probability” approximation to take into account diffusion and absorption into the lines. This approximation was thereafter amply used for the Broad Line Region: this was the beginning of the local escape probability era, with the fiducial paper of Kwan & Krolik (1981) aimed at studying the Broad Line Region of AGN (BLR), immediately followed by a similar computation for the atmospheres of X-ray binary stars (Kallman and McCray 1982) from which xstar derived. The most popular of these codes is cloudy, designed by Ferland (cf. for instance Ferland & Rees 1988) and continuously updated since this time (cf. Ferland 1996, Ferland et al. 1998). These codes are not only used for relatively “cold” clouds like the BLR, but also for warm photoionized media (Broad Absorption Line region in quasars and “Warm Absorbers” in AGN, shock heated media, etc…).

At the other extreme (i.e. for dense, non irradiated, and semi-infinite media) model atmosphere codes have been developed since the fifties with an emphasis put onto the computation of line transfer, which was treated completely, but with the LTE approximation. Non LTE effects were introduced in the sixties, and extensively studied in the seventies, with the state of the art being beautifully exposed in Mihalas book (1978). Since this time many improvments have been made in the numerical methods, in particular aiming at taking into account a large number of atoms and levels, cf. for introductory reviews Rutten 1995, Hubeny 1997.

In the nineties appeared the urgent need for “intermediate” codes, valid for thick or semi-infinite dense media, eventually hot, irradiated by a non thermal continuum extending in the hard X-rays, to cover the whole range of situations encompassed in AGN and in binary stars. Several codes were then built, mainly to compute the em-ission spectrum of accretion disks, irradiated or not by an X-ray continuum. For irradiated disks, they were either of “photoionization type”, using the Kompaneets equation (Ross & Fabian, 1993, based on the Ross 1979 code), or coupling an existing photoionization code to a Monte Carlo computation (Zycki et al. 1994), to take into account Compton diffusions. On the other hand, sophisticated model atmosphere codes were transformed to be applied to accretion disk structure and emission (Hubeny, 1990, Hubeny & Hubeny 1997 and 1998), but without the external irradiation by an X-ray continuum.

Owing to the high optical thickness of the medium in several frequency ranges, such codes require that the transfer of both the continuum and the lines be solved in an “exact” way, i.e. avoiding approximations such as local escape approximation (for the lines) or one stream approximation (for the continuum). Since the medium is generally dense and sometimes close to LTE, they require that all processes and inverse processes be carefully handled. Being irradiated by an X-ray continuum, the medium contains a large number of ionic species, from low to high ionization, which should all be introduced in the computation. Finally, the medium being hot and thick, not only Thomson, but also Compton scattering, should be taken into account.

We have undertaken to build a code in order to satisfy these requirements. Precisely we have built several interconnected codes, which allow more flexibility. The ensemble is far from being perfect and still contains several approximations which restrict its use, but we intend to improve it in the future.

We have presently four codes. One of them (titan) is designed to study the structure of a warm or hot thick photoionized gas, and to compute its emission - reflection - transmission spectrum from the infrared up to about 20 KeV. It solves the energy balance, the ionization and the statistical equilibria, the transfer equations, in a plane-parallel geometry, for the lines and continuum. Then, given the thermal and ionization stratification, the computation of the emitted spectrum from 1 KeV to a few hundre-ds KeV is performed with noar which uses a Monte-Carlo method taking into account direct and inverse Compton scattering (it allows also to study various geometries). pegas adresses the case of more diluted and thin media, and is similar to cloudy with however some differences, and iris is specially devoted to the computation of the line fluxes, including very weak lines, using the most recent available atomic data and a full treatment of the statistical equilibrium equations for a great number of levels. These codes can be used to model a wide variety of astrophysical media, particularly those in which the energy transport is purely radiative, but also plasmas heated by other mechanisms (viscosity, shocks, energetic particules…). They have already been used in several published papers (Collin-Souffrin et al. 1996, Czerny & Dumont 1998, Porquet et al. 1998, Abrassart 1999).

The three codes noar, pegas and iris, are described in other papers (Abrassart 2000, Dumont & Porquet 2000). In this paper we present titan and we discuss a few simple cases, schematized by a slab of gas irradiated on one side or on both sides. As far as possible we will perf-orm comparisons with cloudy in the range of parameters where both codes can be used, to assess the validity of titan. However detailed comparisons are not possible, as they use quite different transfer methods in particular. We will also show some examples of the capabilities of titan in the parameter range which is not accessible to cloudy.

titan is mainly designed to determine the structure (temperature and ionization state) and the continuum em-ission spectrum of a thick hot photoionized slab of gas. Owing to the representation of the ions made in the code it does not give an accurate determination of the weak lines. For a detailed determination of the line spectrum of some ions, it should be complemented by iris.

We briefly summarize below the physical processes (Sect. 2), the transfer method and the iteration procedure (Sect. 3), the thermal equilibrium, (Sect. 4), focussing only on the aspects which are not treated in a standard way. The coupling of titan and noar is briefly described in Sect. 5, and some applications are presented in Sect. 6.

## 2 Physical processes

In titan the physical state of the gas (temperature, ion abundances and level populations of all ionic species) is computed at each depth, assuming stationary state, i.e. local balance between ionizations and recombinations of ions, excitations and deexcitations, local energy balance (equality of heating due to absorption and cooling due to local emission), and finally total energy balance (equality between inward and outward fluxes).

Due to the large range of density and temperature inside the medium, many physical processes play a role at some place and should therefore be taken into account. Ionization equilibrium equations include radiative ionizations by continuum and line photons, collisional ionizations and recombinations, radiative and dielectronic recombinations, charge transfer by H and He atoms, the Auger effects, and ionizations by high energy electrons arising from ionizations by X-ray photons. Energy balance equations include free-free, free-bound and line cooling, and Compton heating/cooling.

The emission-absorption mechanisms for the continuum include free-free and free-bound processes, two-photon process, and Thomson scattering. Special care is given to recombinations to ground state, which are very important from an observational point of view in the X-ray range. They are treated differently according to the relative values of $`kT`$ to the photon energy bin, in order to get an accurate frequency dependence.

Hydrogen and hydrogen-like ions are treated as 6-level atoms. Levels 2s and 2p are treated separately, while full l-mixing is assumed for higher levels. All processes including collisional and radiative ionizations and recombinations are taken into account for each level (cf. Mihalas 1978). Recombinations onto levels $`n>5`$ are not taken into account, which amounts to assuming that the higher levels are in LTE with the continuum, which is generally true in the conditions for which this code is presently used. In the future we plane to add several other levels, and to sum the contributions of the higher levels as it is done for instance in cloudy. Level populations are then obtained as usual by matrix inversion.

In order to save computation time, non H-like ions are presently treated with a rough approximation: interlocking between excited levels is neglected and populations of the excited levels are computed separately using a two-level approximation. This approximation does not predict correctly the details of the line spectrum, since it neglects subordinate lines. Nevertheless recombinations onto excited states are taken into account in the ionization equilibrium and the transfer of these photons is treated in an approximative way as proposed by Canfield & Ricchiazzi (1980). We assume also that each recombination produces a resonant photon after cascades.

The gas composition include 10 elements (H, He, C, N, O, Ne, Mg, Si, S, Fe), and all their ionic species are taken into account. Photoionization cross sections are fitted from Reilman & Manson (1979) and Band L.M. et al. (1990), these values being correct as far as neutral and once ionized ions are not concerned, and it is the case in hot media. For total radiative and dielectronic recombination rates, we use Aldrovandi & Péquignot’s data (1973). When possible, collisional excitation rates are taken from the Daresbury Report (1985). Most of data for iron come from Arnaud & Raymond (1992), Kaastra & Mewe (1993), and from Fuhr et al. (1988). Ionizations by high energy electrons arising from ionizations by X-range photons are taken from Bergeron & Souffrin (1973). Inverse processes (except dielectronic recombinations) are computed through the equations of detailed balance. In the case of a gas close to LTE, we neglect dielectronic recombinations for consistency to insure the balance for each process. All induced processes are taken into account.

The equations are not recalled as they have been given in previously quoted papers.

Iron K lines

These lines require special attention as they are intense in Seyfert nuclei, and they will be observed in detail in the future with Chandra and XMM. Though the iron K lines constitute a complex system described in Band D.L. et al. (1990), we assume presently only one “mean” line per ion with an oscillator strength equal to 0.4, as suggested by Band. More detailed computations are not required as far as the Doppler broadening of the lines (or Compton broadening as well, see Abrassart 2000) is much larger than the distance between the lines, as it is the case in AGN.

In other computations of the line fluxes it is generally assumed that resonant trapping of K$`\alpha `$ photons of FeXVII to Fe XXIII (Ross & Fabian 1993) suppress completely these lines when the Thomson thickness of the emitting medium is larger than a given value, of the order of 0.02 (Zycki & Czerny 1994). Here the transfer of these lines is handled in a standard way, with an additional term included in the statistical and ionization equations to take into account the competition between radiative deexcitation and the Auger process. The population $`N_i^k`$ of the upper level of a transition K$`\alpha `$ of the ion $`i`$ is given by:

$$N_i^kA_{ki}=y_{i1}(N_{i1}K_{i1}+N_iB_{ik}J_{ik})$$
(1)
while the ionization equation for the same ion $`i`$ writes:

$`N_i[(P_i`$ $`+`$ $`K_i)+B_{ik}J_{ik}(1y_{i1})]`$
$`=`$ $`N_{i+1}\alpha _iN_{i1}K_{i1}(1y_{i1}).`$
$`A_{ki}`$ and $`B_{ik}`$ are the Einstein coefficients of the K$`\alpha `$ line, $`J_{ik}`$ is the mean intensity integrated over the line profile, $`y_i`$ is the fluorescent yield, $`\alpha _i`$ is the recombination coefficient, and $`P_i`$ (respt. $`K_i`$) is the photoionization rate of the ground state (respt. of the K-shell) of the ion $`i`$. The second term on the left side of Eq. 2 is generally only of the order of $`10\%`$ of the total photoionization rate.

## 3 Radiation transfer

An important quantity will be used all along this paper, the illumination (or “ionization”) parameter. Among several definitions used in the literature, we adopt the following:

$$\xi =\frac{4\pi F_{\mathrm{inc}}}{n_\mathrm{H}},$$
(3)
where $`F_{\mathrm{inc}}`$ is the frequency integrated flux incident on one side of the slab and $`n_\mathrm{H}`$ is the hydrogen number density at this surface. Note that it does not preclude the possibility of having also a flux incident on the back side of the slab. We call the attention on the fact that all along this paper we will mimic a semi-isotropic and not a mono-directional illumination.

### 3.1 Present method for the transfer of the continuum

We want to study, on one hand media with an optical thickness larger than unity in a range of frequency rich in emission or in absorption lines, on the other hand media with a very inhomogeneous structure. For instance if the illuminated side has a temperature of a few 106 K, while the temperature is only a few 104 K on the back side of the cloud, the spectrum emitted by the illuminated side is mainly formed in hot layers where the absorption coefficient is weak at all wavelengths and Thomson diffusion is dominant (we shall call it the “reflected spectrum” though it does not correspond to pure reflection). On the contrary the spectrum emitted by the backside of the slab is formed in cold layers where absorption must be carefully taken into account (we shall call it the “outward spectrum”). Note also that, even if the incident illumination is mono-directional, the “transmitted spectrum” is almost isotropic since the Thomson thickness, $`\tau _\mathrm{T}`$, is larger than unity. It is therefore not differentiated from the “outward spectrum” even when the incident source is not located on the line of sight. It is possible to estimate an approximate transmitted spectrum using the expression $`F_\nu ^{\mathrm{transm}}=F_\nu ^{\mathrm{inc}}exp(\tau _\nu )`$, where $`F_\nu ^{\mathrm{inc}}`$ is the incident flux, in an aim of comparison.

The transfer is treated with the Eddington two-stream approximation, i.e. the intensity is assumed to be constant in each hemisphere. The transfer equations can then be written :

$`{\displaystyle \frac{1}{\sqrt{3}}}{\displaystyle \frac{dI_\nu ^+}{dz}}`$ $`=`$ $`(\kappa _\nu +{\displaystyle \frac{\sigma }{2}})I_\nu ^++{\displaystyle \frac{\sigma }{2}}I_\nu ^{}+ϵ_\nu `$ (4)
$`{\displaystyle \frac{1}{\sqrt{3}}}{\displaystyle \frac{dI_\nu ^{}}{dz}}`$ $`=`$ $`(\kappa _\nu +{\displaystyle \frac{\sigma }{2}})I_\nu ^{}+{\displaystyle \frac{\sigma }{2}}I_\nu ^++ϵ_\nu `$
where $`z`$ is the distance to the illuminated edge, $`\kappa _\nu `$ is the absorption coefficient, $`\sigma `$ is the diffusion coefficient - here it is due to Thomson scattering - and $`ϵ_\nu `$ is the emissivity (all these quantities are local and depend on the frequency $`\nu `$ except $`\sigma `$). Note that this approximation is closer to the semi-isotropic case than to the normal case. We have already mentioned that it is appropriate for the reflected and for the transmitted flux when the Thomson thickness is larger than unity. It is also appropriate for a semi-isotropic illumination, if the source of radiation is extended (like for instance in the disk-corona model of AGN, or in the blob model of Collin-Souffrin et al. 1996).

The mean intensity $`J_\nu `$ is equal to $`(I_\nu ^++I_\nu ^{})/2`$, while the flux F defined as usual by $`I_\nu \mathrm{cos}\theta d\omega `$ is equal to $`(I_\nu ^+I_\nu ^{})2\pi /\sqrt{3}`$ ; so that the reflected flux is equal to $`I_\nu ^{}(0)2\pi /\sqrt{3}`$, and the outward flux to $`I_\nu ^+(H)2\pi /\sqrt{3}`$. The optical depth and the total optical depth are defined as :

$$\tau _\nu (z)=_0^z\sqrt{3}(\kappa _\nu +\sigma )𝑑z^{}\mathrm{and}T_\nu =\tau _\nu (H)$$
(5)
where $`H`$ is the total thickness of the slab, and the source function as:

$$S=\frac{(ϵ_\nu +\sigma J_\nu )}{(\kappa _\nu +\sigma )}$$
(6)
Both sides of the slab can be illuminated by an external radiation, so the boundary conditions are, at $`z=0`$:

$$I_\nu ^+(0)=\frac{\sqrt{3}}{2\pi }F_\nu ^{\mathrm{inc}},$$
(7)
and at $`z=H`$:

$$I_\nu ^{}(H)=\frac{\sqrt{3}}{2\pi }F_\nu ^{\mathrm{back}}$$
(8)
where $`F_\nu ^{\mathrm{back}}`$ is the back-side flux (it is equal to zero in many of the following computations).

The formal solution of the transfer equations between $`z\delta z`$ and $`z`$ gives:

$$I_\nu ^+(z)=I_\nu ^+(z\delta z)e^{\delta \tau _\nu }+e^{\tau _\nu }_{\tau _\nu \delta \tau _\nu }^{\tau _\nu }S_\nu (t)e^{+t}𝑑t.$$
(9)
Assuming that $`S_\nu (t)`$ is nearly constant in the interval $`z,z+\delta z`$ , one gets:

$$I_\nu ^+(z)=I_\nu ^+(z\delta z)e^{\delta \tau _\nu }+\frac{1e^{\delta \tau _\nu }}{2}[S_\nu (z\delta z)+S_\nu (z)]$$
(10)
with a similar equation for $`I_\nu ^{}(z)`$.

Guessing the initial value of $`I_\nu ^{}(0)`$ fails because calculations diverge, unless this initial value is provided with an extreme precision: the predicted value should differ from the real one at all frequencies by less than $`10^6`$. This problem has also been pointed out by Coleman (1993). We therefore start to compute $`I_\nu ^{}(z)`$ from the back side with:

$$I_\nu ^{}(z)=I_\nu ^{}(z+\delta z)e^{\delta \tau _\nu }+\frac{1e^{\delta \tau _\nu }}{2}[S_\nu (z+\delta z)+S_\nu (z)]$$
(11)
o Thus we calculate $`S_\nu (z)`$ and the mean intensity $`J_\nu (z)`$ from Eq. 6 and:

$`J_\nu (z)`$ $`=`$ $`{\displaystyle \frac{I_\nu ^+(z\delta z)+I_\nu ^{}(z+\delta z)}{2}}e^{\delta \tau _\nu }`$
$`+`$ $`{\displaystyle \frac{(1e^{\delta \tau _\nu })}{4}}[S_\nu (z\delta z)+2S_\nu (z)+S_\nu (z+\delta z)]`$
assuming a constant $`\delta z`$ (actually it is not constant so the formulae are more complicated but do not deserve to be given here).

The procedure is then the following. We divide the sl-ab into a set of plane-parallel layers with the density $`n_\mathrm{H}`$ given in each layer. Note that the code can easily be coupled with a prescription for the density or for the pressure, such as a constant pressure, or a pressure given by hydrostatic equilibrium. Presently the slab is divided in about 300 layers. The geometrical thickness of each layer is sm-aller close to the surfaces than in the middle of the slab because of the rapid variation of the physical paramet-ers (the optical thickness, the temperature, the ionization state).

\- For each layer, starting from the illuminated side, where $`I_\nu ^+(0)`$ is given by Eq. 7, the ionization and thermal balance equations are solved. The source function, the opacity and the emissivity, are computed. As they depend on $`J_\nu `$, we use Eq. 3.1 with the value of $`I_\nu ^+(z\delta z)`$ of the previous layer and the value of $`I_\nu ^{}(z+\delta z)`$ provided by the previous iteration and given by Eq. 11 (see below).

\- $`I_\nu ^+(z)`$ is transferred through the layer according to Eq. 10, while $`S_\nu (z)`$ and $`\tau _\nu (z)`$ are buffered for each layer and each frequency.

\- when the back side of the cloud is reached, new values of $`I_\nu ^{}(z)`$ are calculated, starting from the back side where $`I_\nu ^{}(H)`$ is given by Eq. 8, and using Eq. 11 with the previous values of $`S_\nu (z)`$ and $`\tau _\nu (z)`$ ; actually we use the values given by several previous iterations to accelerate the procedure.

\- the whole calculation is repeated until convergence. It is stopped when energy balance is achieved for the whole slab (i.e. when the flux entering on both sides of the slab is equal to the flux coming out from both sides).

The drawback of this method is that it requires a large number of iterations to get a complete convergence of the model. The iteration starts by assuming given values for $`I_\nu ^{}(z)`$. If this initialisation is not adequate, convergence problems are encountered. We have chosen for the first iteration $`I_\nu ^{}(z)`$ equal to $`I_\nu ^+(z)`$ multiplied by a factor of the order of unity, which works reasonably well. Actually the convergence of $`I_\nu ^{}(z)`$ (as well as of the temperature and the ionization state) is reached within less than 1$`\%`$ after only a relatively small number of iterations in the layers corresponding to about one Thomson thickness. This property is interesting as far as the outward spectrum is negligible, which is the case in the EUV and X-ray range. It can also be used to determine the ionization state and the temperature of the layers where Compton diffusion takes place (cf. Sect. 5 about coupling with the code noar).

Fig. 1 displays the reflected continuum obtained after different numbers of iterations (for clarity we prefer to show here the continuum and not the complete spectrum). The model is a slab with the following characteristics:

\- a constant density 1012 cm-3,

\- a column density 1026 cm-2;

\- it is illuminated on one side by an incident power law continuum proportional to $`\nu ^1`$ from 0.1 eV to 100 KeV,

\- the illumination parameter is $`\xi =10^3`$ erg cm s-1;

\- there is no illumination on the other side;

\- the abundances are (in number): H: 1, He: 0.085, C: 3.3 10-4, N: 9.1 10-5, O: 6.6 10-4, Ne: 8.3 10-5, Mg: 2.6 10-5, Si: 3.3 10-5, S: 1.6 10-5, Fe: 3.2 10-5.

In the following we will call this our “reference model”. All along the paper we will use the same density and the same spectral distribution of the illuminating radiation, so we will not specify these parameters anymore. Note that, for this peculiar spectral distribution, $`\xi `$ and the ratio of incident ionizing photon number to hydrogen densities, U, are related by : $`\xi =8.21U\mathrm{ln}(\nu _{\mathrm{max}}/\nu _{\mathrm{min}})`$ where $`\nu _{\mathrm{max}}`$ and $`\nu _{\mathrm{min}}`$ are the maximum and minimum photon frequencies.

We want to stress here the importance of the backward intensity. It cannot be neglected as it participates to the ionization, to the level population, and to the energy balance. For instance in a purely scattering medium the returning flux varies as $`F_{\mathrm{inc}}\frac{\tau _\mathrm{T}/2}{1+\tau _\mathrm{T}/2}`$, so when $`\tau _\mathrm{T}`$ is much larger than unity, all the incident flux is reflected. Thus the temperature of the first layer at the illuminated side (which depends on the sum of incident plus returning radiation) should increase with the column density for a given incident flux (except possibly in cases of very high values of the illumination parameter).

Fig. 2 shows the dependence of the surface temperature on the column density $`N`$. The model is similar to the reference one except that the illumination parameter is equal to 100 erg cm s-1 and the column density is varying. We see that the surface temperature is constant up to about $`N=10^{22}`$ cm-2, but after it increases from 2 105 K to 2.7 105 K for $`N`$ increasing from $`10^{22}`$ to $`10^{24}`$ cm-2. It saturates when a large fraction of the slab becomes close to LTE, and the returning flux becomes constant. With cloudy (version 9004) the surface temperature does not depend on the column density (for a correct comparison the illumination parameter in cloudy must be multiplied by a factor $`\sqrt{3}`$ to account for the semi-isotropic approximation used in titan).

Fig. 3 displays the dependence of the surface temperature on an assigned ratio $`I^{}(0)/I^+(0)`$ independent of the frequency. We see that a modest returning flux (20$`\%`$) is sufficient to account for the 60$`\%`$ increase of the surface temperature reached for $`N=10^{24}`$ cm-2, with respect to a very small column density. To understand this behaviour, one should remember that the spectral distribution of the reflected flux contains a much larger proportion of EUV radiation than the incident continuum, in particular in the lines. So the illumination parameter is increased by a larger amount than the ratio of the integrated reflected over incident continuum. We can also guess that the increase of the temperature at the illuminated surface with respect to cloudy will be compensated by a decrease at the back side, because of energy conservation (cf. Fig. 14).

### 3.2 Present method for line transfer

Radiative transfer of the lines is treated in the same way as the continuum, the line profile being represented by several frequencies distributed symmetrically around the line center.

We assume that even if large macroscopic velocities are present, they do not play any role in the line transfer, in other words that the thickness of the slab is much smaller than the scale length of the velocity gradient. We also assume that the lines do not overlap, which is valid as far as multiplets are treated globally, and if the microscopic turbulent velocity is not much larger than the thermal velocity (doing this we neglect line fluorescence due to the overlap of two lines, like the HeII Balmer$`\beta `$ line and the hydrogen Ly$`\alpha `$ line, whose wavelengths differ only by a fraction 4 10-4).

Finally we assume presently complete redistribution in all the lines. We intend to include partial redistribution for resonance lines like Ly$`\alpha `$ in the near future. For the moment partial redistribution in these lines can be mimicked by assuming complete redistribution in a pure Doppler profile, as the line photons diffuse more in space than in frequency, owing to the absence of enlargment of the lower level. A Doppler profile can also be used for strongly interlocked lines, and in any case, it does not differ strongly from a Voigt profile when $`\tau _0`$ is smaller than $`1/a`$ so the medium is optically thin in the damping wings. Note that partial redistribution is less important in dense media, like irradiated accretion disks in AGN, than in dilute media, owing to collisional redistribution.

The line profile $`\varphi _\nu `$ is a symmetrical Voigt function. The Doppler width $`\mathrm{\Delta }\nu _\mathrm{d}`$ and the damping constant depend on the temperature and on the density, and consequently vary with $`z`$. It may happen that the temperature decreases by two orders of magnitude from the illuminated to the back side, so $`\mathrm{\Delta }\nu _\mathrm{d}`$ decreases by one order of magnitude. This must be taken into account to choose the points in the profile. To solve the statistical equilibrium equations one needs to know the total line intensity weighted by the profile $`J_\nu \varphi (\nu )𝑑\nu `$, and to treat the thermal and ionization equilibria one simply integrates the line intensity over the frequencies $`J_\nu 𝑑\nu `$. For intense lines the first integral is dominated by the Doppler core, and the second by the wings. In the case of a saturated line with strong wings, the set of frequencies chosen to describe the line must cover a range from a few tenths $`\mathrm{\Delta }\nu _\mathrm{d}`$ to a few hundreds $`\mathrm{\Delta }\nu _\mathrm{d}`$ and at least 10 points are necessary for a correct representation of the profile.

Fig. 4 displays for the reference model the variation with $`z/H`$ of the OVIII Ly$`\alpha `$ profile in both directions, $`I_\nu ^{}`$ and $`I_\nu ^+`$, $`J_\nu `$ being the half sum of these profiles. We see that the profile varies considerably from the surface to the deeper layers.

Now we are able to show the difference between the transfer treatment and the escape probability approximation for the lines discussed in the Appendix. In this aim, we have compared the divergence flux $`\rho `$ ( Eq. 17) with the escape probability towards the illuminated side $`P_e(\tau _0)`$, and with an approximate expression proposed by Rees et al. (1989) to take into account continuum absorption in the statistical and in the energy balance equations :

$$\beta _{\mathrm{eff}}=\frac{\kappa _c}{\kappa _c+\kappa _l}+\frac{\kappa _l}{\kappa _c+\kappa _l}\frac{[P_e(\tau _0)+P_e(T_0\tau _0)]}{2}$$
(13)
where $`\beta _{\mathrm{eff}}`$ is the total effective escape probability, $`\kappa _l`$ and $`\kappa _c`$ are the absorption coefficients respectively at the line center and in the underlying continuum, $`T_0`$ is the total optical depth of the slab at the line center, and $`P_e(\tau _0)`$ (respt. $`P_e(T_0\tau _0)`$) are the escape probabilities towards the illuminated (respt. the back) side.

$`P_e(\tau _0)/2`$ and $`\beta _{\mathrm{eff}}`$ have been computed with titan and compared to $`\rho `$ in Fig. 5, as functions of the distance $`z`$ to the illuminated surface, for the reference model.

We see on this figure that $`\beta _{\mathrm{eff}}`$ and $`P_e(\tau _0)/2`$ are almost always identical. It is first due to the fact that the slab is very thick, so $`P_e(\tau _0)`$ is much larger than $`P_e(T_0\tau _0)`$ except very close to the back side, and second because $`\kappa _c`$ is always $`\kappa _l`$, as all the lines displayed on the figure are intense. But $`\beta _{\mathrm{eff}}`$ differs considerably from $`\rho `$ in the region where the lines are formed (i.e. the region where $`\rho `$ is not equal to zero and where the ion is present). In particular $`\rho `$ is negative in an important fraction of the line formation region for H$`\alpha `$, L$`\alpha `$, and OVIII L$`\alpha `$, while $`\beta _{\mathrm{eff}}`$ is a positive quantity.

As a consequence the line spectrum is not correctly computed when escape probabilities are used. This can be seen on Table 1 which gives the equivalent widths (negative if the line is in absorption, positive if it is in emission) of some intense lines from both sides of the slab. Several of these lines are in absorption, either from the illuminated or from the dark side, which is impossible with the escape probability approximation.

Presently the radiative transfer is solved for 300 frequencies in the continuum from 0.01 eV to 25 keV, including all the ionization edges, and for 440 spectral lines including the fluorescence Fe K lines, the integrals over the line profiles being achieved using 15-point Gauss-Legendre quadrature, together with a change of variable.

Remarks on the convergence of the method

It is well-known that the lambda iteration method used here (actually mainly for historical reasons, because we st-arted to build a nebular code with the line escape probability approximation) converges extremely slowly for intense lines in stellar atmospheres, and we intend to replace it by the Accelerated Lambda Iteration method as soon as possible. The convergence of $`J_\nu `$ is achieved after 600 iterations, for all lines except a few ones and for all continua. This can be checked on Fig. 6 which displays for the same lines as Fig. 5 the profiles of the lines emitted by the bright side, after different number of iterations, for the reference model. As a consequence, the structure of the slab (T, ionization) does not change at all after a few hundreds iterations. Indeed our case is more favorable than stellar atmospheres. The explanation is double: 1. the relatively small density of the medium compared to a stellar atmosphere, which decreases comparatively the influence of collisional excitations, and 2. the fact that in a photoionized plasma the temperature is low for a high ionization degree, so the excited levels lie at very high energies compared to the thermal energy. Consequently the populations of these highly ionized ions are smaller than at Boltzman equilibrium with the ground level, and the lack of convergence for excited level populations does not have a strong influence on the ionization state of the corresponding ion. We have also checked that the same structure is obtained for different initial boundary conditions. All this implies that the continuum emitted from both sides (including for instance the Lyman discontinuity) is completely converged, and only the few non converged lines are not well computed, in particular those emitted by the back side (but they are negligible in the whole energy balance). Since we do not claim to compute a correct line spectrum before the subordinate lines are implemented in the code, we think that this approximation is sufficient for our present purpose.

## 4 Energy balance

titan has been designed for stationary radiatively heated media, so the energy conservation implies the balance of radiative flux. However it can be used also for other cases out of radiative equilibrium, for instance by imposing the temperature of the medium, or by adding a non radiative energy in each layer, such as the viscous flux of an accretion disk. We describe here how the code works in a purely radiative case.

We first stress the fact that the energy balance should be realized both locally and globally.

The local energy balance equation writes:

$$\frac{dF_\nu }{dz}𝑑\nu =0=n_en_H[\mathrm{\Gamma }_{\mathrm{tot}}+\mathrm{\Lambda }_{\mathrm{tot}}],$$
(14)
where $`n_en_H\mathrm{\Gamma }_{\mathrm{tot}}`$ and $`n_en_H\mathrm{\Lambda }_{\mathrm{tot}}`$ are the total heating and cooling rates per unit volume. The temperature of a layer is computed through an iteration process until $`\mathrm{\Gamma }_{\mathrm{tot}}`$ and $`\mathrm{\Lambda }_{\mathrm{tot}}`$ are equal.

As an example, Fig. 7 gives the dependence of the cooling and heating rates, and of the equilibrium temperature, on the distance from the illuminated side $`z`$, for the reference model. We note that in this model Compton heating-cooling is not very important, and that bound-free contributes as a heating term, compensated by a cooling term due to lines and free-free processes. Note also that a quasi LTE is reached at a few Thomson depths in this model (it is reached at a smaller optical depth for a sm-aller value of $`\xi `$). Above this value of $`z`$ the temperature is almost constant, except close to the back side where it decreases more rapidly due to the leakage of photons.

The global energy balance is reached through an iteration procedure, when the total flux emitted by both sides of the slab, $`F_{\mathrm{tot}}^{\mathrm{em}}`$, is equal to the total flux incident on both sides, $`F_{\mathrm{tot}}^{\mathrm{inc}}`$. Actually this is the most important convergence criterium and it requires many iterations. To illustrate the problem, Fig. 8 displays the flux of a few intense emission lines and of the total flux in the lines and continuum, for the reference model, as functions of the number of iterations. We see that after about 200 iterations the global energy balance is reached within about one percent.

In the same way, Fig.9 shows that the energy balance with titan is reached within one percent after 200 iterations for two models with a column density equal to 1025 cm-2, and with an illumination parameter equal to $`\xi =100`$ and $`\xi =1000`$ erg cm s-1. cloudy or xstar have been used several times in this range of parameters (Martocchia & Matt, 1996, Zycki et al. 1994, for instance), but it is not sure that the global energy balance is achieved owing to the escape treatment of the lines, as discussed above.

Finally we must recall that the neglect of subordinate lines (except in hydrogen-like ions) could lead to a substantial increase of the energy losses and consequently to a change of the thermal balance, as these lines have a smaller optical thickness than resonant lines, and therefore can escape more easily from the deep layers. Since in a photoionized plasma the temperature is relatively low and the excited level populations of highly ionized ions are generally small, one could expect only a small influence on the overall spectrum.

To check this point we compare in Figs. 10 and 11 the temperature and the emitted spectrum for the reference model, with and without taking into account interlocking between excited levels in hydrogen like ions (note that interlocking is included for hydrogen and helium in both cases). These ions are the most abundant ones in the region emitting the reflected spectrum, in our reference model. Fig. 10 shows that the temperature is slightly shifted in layers corresponding to a Thomson thickness of a few. As a consequence the emitted spectrum is also different, but not by an important amount. We can infer a similar behaviour for smaller illumination parameters, even if the abundant ions are less ionized, as the temperature is also lower. For a higher illumination parameter the lines contribute less to the energy balance.

Since the detailed computed line spectrum of hydrogen-like ions is different with the two treatments, contrary to the overall spectrum, it is mandatory to treat correctly the excited levels for all the ions in order to get a correct line spectrum. Note that in this respect cloudy is presently doing much better than titan, as it includes subordinate lines and more levels for each ion.

Compton heating and cooling

The Compton heating-cooling term is equal to:

$$\mathrm{\Gamma }\mathrm{\Lambda }=\frac{\sigma _T}{m_ec^2}\frac{1}{n_H}4k(T_{\mathrm{Comp}}T)4\pi J_\nu 𝑑\nu $$
(15)
where $`T_{\mathrm{Comp}}`$ is the Compton temperature. We give here this well-known equation only to recall that it depends (as well as $`T_{\mathrm{Comp}}`$) on the spectral distribution of the local flux, so it requires a careful computation of this flux (including the lines).

Below 25 KeV, $`J_\nu `$ is provided directly by titan in each layer. Above 25 KeV, $`J_\nu `$ is not calculated by titan, but by the code noar, as the energy shift of the photons due to Compton scattering is no more negligible. In this case the results from noar are also used to compute the exact value of the Compton heating-cooling term. This is discussed in more detail in the next section.

## 5 Coupling with the Monte Carlo code

In order to compute the effect of Compton scattering on the high energy part of the reprocessed radiation, and particularly on the reflected spectrum above 1 keV, and to determine precisely the Compton heating-cooling rate entering in the energy balance equation, we need another numerical approach. The transfer scheme in titan is based on frequency arrays, and inelastic scattering would introduce a redistribution on energy bins making the code much heavier. In this aim Abrassart (2000) developped a Monte Carlo code, noar. The asset of such an approach is also that it enables to investigate an arbitrary geometry and to determine the angular dependance of the observed spectra. These aspects were difficult to tackle with a code such as titan. Moreover, it allows to easily extract time variability information.

noar is described in detail in Abrassart (2000). We explain here only how it is coupled with titan. It shares the same opacity data and includes the same set of 102 ions, for the 10 elements of highest cosmic abundances. Given all the ionisation degrees provided by titan, it computes absorption cross sections in each layer. Free-free absorption is also taken into account. noar includes the recombination continua of hydrogen and helium like ions, but does not include line emission, except fluorescence lines. The proper yields to account for the competition with the Auger effect are used for these lines. Fluorescence of Fe XVII-XXIII is suppressed by resonnant trapping. For the pseudo fluorescence of H and He-like species, a simplifying prescription is adopted which includes only the two most probable outcomes of a K shell photoionization, i.e. direct recombination on ground level or L-K transition. The spectral distribution of the recombination continuum is determined by the local temperature in the “current” slab.

noar takes into account direct and inverse Compton scattering. The method used for modelling Compton scattering is basically the one described in Pozdniakov, Sobol & Sunyaev (1983) and in Gorecki & Wilczewski (1984), with a different use of statistical weight, because of the competition with photoelectric absorption. When the high energy cut-off of the incident continuum is of the order of 100 KeV, the reflected spectrum above 10 KeV exhibits a “Compton hump” peaking at about 30 KeV (Lightmann & White 1988). This spectral region of optimal continuum albedo occurs where the sum of photoelectric absorption and Compton losses is minimum. This feature depends weakly on the ionization state, although its prominence in an observed X-ray spectra, and notably its level with respect to the iron line, is sensitive to $`\xi `$ (see Fig. 19) as the soft X-ray albedo. Note that the shape and extension of the hump depends on the high energy cut-off of the incident spectrum, which is 100 KeV all along this paper.

Fluorescent lines are significantly Compton broadened for rather high illumination parameter, above $`\xi 10^3`$ erg cm s-1. The broadening is asymmetyric, the profile is skewed toward the red. At still higher $`\xi `$ and temperature, a blue wing appears, due to the inverse Compton effe-ct. The extend and importance of both wings depends on the optical depth of the scattering medium. The red wing depends only on the frequency of the line (i.e. on the ionization state of iron) whereas the blue one only depends on the temperature of the scattering medium.

titan and noar lead to almost identical reflected spectra in the 1-20 KeV range (they differ only by Compton broadening of lines and of photoelectric edges), and noar allows to obtain the spectrum in the higher energy range. This spectrum is used for the global energy balance.

Another use of noar is to provide titan with the local Compton gains and losses in each layer. This is necessary, because Compton heating-cooling is dominated by energy losses of photons $`>`$ 25 Kev, in particular for high values of $`\xi `$ when Compton heating-cooling is dominating the energy balance and should be computed in an exact way.

As an example Fig. 12 gives the dependence of the ratio $`(\mathrm{\Gamma }_{\mathrm{Comp}}\mathrm{\Lambda }_{\mathrm{Comp}})/\xi `$ on the Thomson depth, for different values of $`\xi `$, the other parameters being the same as in the reference model. We see that the ratio does not depend on the illumination parameter if $`\xi 10^4`$ erg cm s-1. This is expected since the Compton heating-cooling is dominated by high energy photons, which do not care about the ionization state. It also means that inverse Compton is negligible as a cooling process (though it has an influence on the line spectrum), otherwise $`\frac{\mathrm{\Gamma }\mathrm{\Lambda }}{\xi }`$ would depend on the gas temperature and therefore on $`\xi `$. This is due to the relatively low temperature of the gas compared to the mean energy of the photons. This property breaks for very high illumination parameters, since low energy photons are then able to penetrate in the deepest layers and to play a role in the cooling rate, and the temperature of the gas is close to $`T_{\mathrm{Comp}}`$.

The deposition of energy shown on Fig. 12 corresp-onds to a semi-isotropic distribution of the incident radiation. It is different in the case of a mono-directional radiation. This is shown on Fig. 13 which displays $`(\mathrm{\Gamma }_{\mathrm{Comp}}\mathrm{\Lambda }_{\mathrm{Comp}})/\xi `$ for a neutral gas in the case of an isotropic and of a perpendicular illumination. When the radiation field is semi-istropic there is a rapid decrease of $`\frac{\mathrm{\Gamma }\mathrm{\Lambda }}{\xi }`$ for small values of $`\tau _\mathrm{T}`$, due to the small distance covered by the photons. This is compensated by a smaller (by about a factor two) deposition of energy in the deeper layers.

The Compton heating-cooling rate obtained with noar is then fitted analytically as a function of $`z`$, which is transferred to titan.

## 6 Some applications

A full grid of models will be published elsewhere (Dumont & Abrassart 2000). Here we discuss only a few interesting cases.

### 6.1 Optically thin hot medium: the Warm Absorber in AGN

The Warm Absorber (WA) is a hot medium located on the line of sight of the X-ray source in AGN. It is responsible of the absorption edges of OVII and OVIII observed in their soft X-ray spectrum in about 50$`\%`$ of them. It might also constitute the “mirror” invoked in the Unified Scheme to account for the polarized Broad Lines observed in Seyfert 2 galaxies (Antonnucci & Miller 1985). It has been modelled in many papers, using ion (Netzer 1993, 1996), xstar (Krolik & Kriss 1995), cloudy (Nicastro, Fiore & Matt 1999 for the most recent paper), pegas and iris (Porquet et al. 1999). This medium has a column density of 1021 to 1024 cm-2, and a typical illumination parameter of $`\xi =100`$ erg cm s-1. Porquet et al. have shown that its density should be at least 1010 cm-3 to avoid the emission of too strong coronal lines.

As the WA is relatively thin ($`\tau _\nu 1`$ for the continuum in a large range of frequencies), it is possible to use the line escape probability formalism. This is however limited to the case where $`\tau _\nu `$ is smaller than unity in the continuum underlying intense lines. This specific problem will be discussed in detail in the paper describing the code pegas (Dumont & Porquet 2000). A correct computation of the backward flux is also required. We have shown indeed that the temperature depends on the intensity of this backward flux.

We have run several models corresponding to the conditions of the WA, i.e. for $`\xi =100`$ erg cm s-1 and column densities varying from 1022 to 1024 cm-2 (the other parameters being the same as in the reference model). We have compared the results to those obtained with cloudy for the same models, to check that they are similar at low column densities, and to determine at which column densities they begin to differ strongly. Note that a detailed comparison between titan and cloudy is not possible, as these codes differ both in the transfer method and in the atomic data (for instance owing to the semi-isotropic illumination, a given column density corresponds to a larger optical thickness with titan than with cloudy), so only trends can be obtained.

Fig. 14 displays the dependence of the temperature on $`z`$, obtained with titan and with cloudy. As discussed previously, it is higher with titan at the illuminated side for a high column density, owing to the returning flux, and it is lower at the dark side.

The ionization state (of fundamental importance for the study of the WA) and the emission spectrum, are also different with titan and with cloudy at high column densities. Fig. 15 shows the ionization state of Carbon, Oxygen and Iron for a column density of 1024 cm-2, computed with titan and with cloudy. While the ionization state close to the illuminated side is quite similar (in particular if one remembers that there is a factor $`\sqrt{3}`$ in the optical depth), it is very different in the deepest layers, where cloudy leads to a smaller degree of ionization than titan (CIII and OIII instead of CIV and OIV). This effect cannot be explained by the difference in optical thickness, which acts in the opposite direction. It is certainly due to the different treatments of the diffuse continuum, and perhaps of the lines. With titan we have checked that the diffuse continuum at large depths is sufficient to ionize OIII into OIV. To illustrate this point Fig. 16 displays the variation of the flux at 54.9 eV (the OIII-OIV ionization edge) as a function of the depth. It shows that the diffuse flux is still very intense and sufficient to ionize OIII at large depths (contrary to the transmitted flux which is very weak). Note that this result depends strongly on the way the spectral distribution of the diffuse flux is computed. If it would be reduced to a unique frequency at the recombination frequency of He++, 54.4 eV, (as it is generally the case with the “On The Spot” approximation), there will be no photons able to ionize OIII.

Fig. 17 displays the reflected and outward spectra, for column densities 1023 and 1024 cm-2. We see that the overall shape of the continuum is very similar with titan and with cloudy. However, some differences in the detailed features appear for the largest column density. For instance the Lyman edge is in emission in the reflected spectrum with both codes, but it is larger with titan (presumably because the temperature is higher in the em-ission layers). In the outward spectrum it is present as a very weak absorption with cloudy, while it is in emission with titan.

The line spectra are obviously also different. There are much more lines with cloudy, but on the other hand the fewer lines of titan are more intense. This is expected since subordinate lines are not taken into account exce-pt for hydrogen-like ions in titan. As mentioned before, some lines appear in absorption with titan, both in the reflected and in the outward spectrum.

### 6.2 Hot medium optically thick to Compton scattering: the UV-soft X-ray emitting medium in AGN

titan is specially designed for Compton thick photoionized media like those commonly assumed to emit the UV/soft X-ray continuum in AGN and to produce through Compton reflection the FeK $``$7 KeV line and the 30 KeV hump observed in many Seyfert 1 galaxies.

Although the exact nature of this medium is not known (an irradiated accretion disc, as proposed by Ross & Fabian 1993, Zycki et al. 1994, to quote only the first papers on the subject, or a clumpy Compton thick medium, cf. Collin-Souffrin et al. 1996), its characteristics and physical state are comparable: a shell of gas with a column density of at least 1025 cm-2, a temperature of 105-106 K due to radiative heating, and a density spanning a range from 1012 to 1015 cm-3. The corresponding value of $`\xi `$ is typically 300 to 3000 erg cm s-1. Our reference model is therefore representative of this medium. In the case of an irradiated accretion disk, it has been mimicked as a sh-ell of constant density irradiated from above by a power law continuum, and from below by a black body radiation (actually it is not the best way to take into account the presence of the underlying viscously heated disk, cf. Rozanska et al. 2000).

Let us consider first a slab illuminated on one side by a semi-isotropic radiation. Fig. 18 displays $`T`$ versus $`z`$, for different values of $`\xi `$. The other parameters are the same as in the reference model. It is interesting to notice that the quasi constant temperature regime (corresponding to LTE) is reached for increasing values of $`z`$ when $`\xi `$ increases. For instance it is reached for $`z10^{13}`$ cm, corresponding to $`\tau _\mathrm{T}6`$ for $`\xi =1000`$ erg cm s-1. It means that any computation aiming at giving the reflected spectrum of such a model should solve correctly the thermal and ionization equilibrium until these deep layers.

Fig. 19 shows the reflected spectrum for the same models. Note that the resolution adopted here (100) does not correspond to a region where the lines are broadened by large velocity field, either turbulent or organized (rotation). We use this high resolution only to show more clearly the details of the spectrum.

We see that the overall shape of the spectrum computed with titan and noar is very similar in the range 1-20 KeV. However the detailed spectral features are different, particularly for the higher values of $`\xi `$. For instance the complex feature near 7 KeV (which is made of a mixture of several Iron edges and of Iron lines dominated by low ionization stages for $`\xi =100`$ erg cm s-1, and by FeXXV and FeXXVI for $`\xi =10^4)`$ erg cm s-1 is strongly modified by Compton scattering. For $`\xi =10^4`$ erg cm s-1 the line has a large red wing due to direct Compton scattering, and a weak blue wing due to inverse Compton scattering. The absorption edge is also completely erased. The effect is smaller for $`\xi =10^3`$ erg cm s-1, and almost absent for $`\xi =100`$ erg cm s-1, as discussed above. Similar results have been obtained in many papers dealing with the X-ray spectrum of Seyfert galaxies (cf. for instance Ross et al. 1999).

Let us now compare the results obtained by us with titan, and by Zycki et al. (1994), for a slab of constant density 1014 cm-3, illuminated by a power law $`\nu ^{0.9}`$ between 10 eV and 100 KeV on one side, and by a black body $`T_{\mathrm{BB}}=10^5`$K on the other side. We will compare the results for $`\xi `$(Zycki et al.)=300, which corresponds to $`\xi `$(titan)$``$ 200 erg cm s-1, owing to our semi-isotropic illumination, and to a slightly different definition of $`\xi `$.

Like us, Zycki et al. used a Monte-Carlo method to take into account Compton scattering. The thermal and ionization state of the gas is calculated apart with the code xstar, using like cloudy the escape probability approximation. Fig. 20 gives $`T`$ versus $`\tau _\mathrm{T}`$ obtained by us and by Zicky et al. The same result is observed as with cloudy: our temperature is larger near the surface of the slab. It reaches $`T_{\mathrm{BB}}`$ at $`\tau _\mathrm{T}`$=0.7. For a larger value of $`\xi `$, the connection with the underlying black body would be reached in even deeper layers, as discussed above. For $`\tau _\mathrm{T}0.7`$ the temperature becomes smaller than $`T_{\mathrm{BB}}`$, and joins $`T_{\mathrm{BB}}`$ only at $`\tau _\mathrm{T}20`$. This behaviour is however dependent on the thickness of the slab, it is why this representation of the disk is not meaningful (Rozanska et al. 2000). On the contrary Zycki et al. impose the gas temperature to be equal to $`min(T_{\mathrm{BB}},T_{\mathrm{thermal}\mathrm{balance}})`$.

Fig. 21 displays the reflected spectrum, computed by us and by Zycki et al. The resolution is degraded to 20, to fit the energy bins of Zycki et al. The overall shape of the spectra are roughly similar. The major discrepancy is the OVIII edge, in emission with our computation and in absorption with Zycki et al.’s one. We think that it is due to the higher surface temperature and to the lower temperature in the deeper layers obtained in our computation.

## 7 Conclusion

We have shown that the coupling of titan and noar allows to compute the structure and the emission of hot Compton thick irradiated media in an unprecendented way. First it solves consistently both the global and the local energy balance, which is impossible in a thick medium with codes handling the line transfer with the escape probability approximation, as all present photoionization codes do. We have also shown the importance of the returning flux (which is neglected in photoionization codes) even for relatively low column densities. Second it takes into account in an exact way inverse and direct Compton scattering, both in the energy balance and in the computation of the emitted spectrum. Finally, it allows to treat any geometry, open or closed.

Although the problem of the convergence process is not as drastic as in stellar atmospheres, since we are able to get complete converged structure and spectrum in a still reasonable computing time, a most urgent improvement of the code is to accelerate the convergence process through the use of the Accelerated Lambda Iteration method. This will not only allow to get the results in a much smaller time, but also to get convergence for the few lines which are still not converged after about 103 iterations.

Then the following improvements of titan will be to take into account subordinate lines in solving a multi-level atom for all ions and to bring up to date the atomic data. More detailed line spectra for some abundant ions will be obtained through coupling with the code iris. The L shells of Iron which are already taken into account for FeXVII to FeXXII will be implemented as well as a better representation of the Iron K lines. A few elements will be added to the already ten existing ones.

###### Acknowledgements.

We are grateful to the referee, G. Ferland, for his comments which have led to substantially clarify the paper, to M-C. Artru for a very careful reading of the manuscript, and to I. Hubeny for enlightning discussions.

## Appendix : Previous methods

### 7.1 Transfer of the continuum

It is not always clear which approximations are used in the different codes. Basically the computation of the diffuse radiation field seems to be similar in all codes, except in Ross & Fabian (1993) and in subsequent works using this code, where it is computed using the Kompaneets equation. cloudy (Ferland 1996), xstar (Kallman & Krolik, 1995 for the last version), ion (Netzer 1990, 1993), use a modified version of the “on the spot” approximation, which amounts to assuming a kind of escape probability for the diffuse continuum, and one stream approximation in the transfer equation. These approximations are correct only for a continuum optical thickness smaller than unity (effective or total opacity), and when the properties of the cloud do not vary considerably between the point where a photon is emitted and the point where it is reabsorbed. In particular one very important requirement for studying hot and thick media is to take into account the radiation returning from the backside of the slab, even when it is not illuminated. As the medium is optically thick in a large frequency range, this radiation is intense and modifies the physical state of the whole slab, including the illuminated side(see Fig. 2).

### 7.2 Line transfer: the escape probability approximation

All previous photoionization codes, except that of Collin-Souffrin and Dumont (S.) (1986), treat line transfer by the so-called local escape probability formalism, which uses the probability that a line photon emitted at a given point can escape in a single flight from the cloud, assuming that the rest of the cloud is homogeneous and has the same properties as the emitting layer. It amounts to identifying the divergence flux of the statistical equilibrium equations with the escape probability intervening in the line emerging flux. This approximation is valid only in the case of complete frequency redistribution, absence of line interlocking, and if the medium is homogeneous. Its use can have severe consequences on the emission line spectrum and on the energy balance when these conditions are not fulfilled as it was often discussed, see for example by Collin & S.Dumont (1986) or Elitzur (1984). We introduce here the divergence flux in the aim to show that our results are completely different from those obtained by the escape probability formalism.

To simplify the discussion, let us consider in the rest of this section a simple two-level atom. The level population balance writes:

$`n_u(A_{ul}+B_{ul}{\displaystyle J_\nu \psi _\nu 𝑑\nu }+n_\mathrm{e}C_{ul})`$ (16)
$`=n_l(B_{lu}{\displaystyle J_\nu \varphi _\nu 𝑑\nu }+n_\mathrm{e}C_{lu})`$
where $`\varphi _\nu `$ and $`\psi _\nu `$ are the absorption and emission line profiles, $`A_{ul}`$, $`B_{ul}`$ and $`B_{lu}`$, $`C_{ul}`$ and $`C_{lu}`$ are the usual radiative (Einstein) and collisional excitation and deexcitation coefficients, and $`J_\nu `$ is the angle averaged intensity.

Let us now define the divergence flux of the transition, $`\rho _{ul}`$:

$$\rho _{ul}=\frac{\left\{n_u(A_{ul}+B_{ul}J_\nu \psi _\nu 𝑑\nu )n_l(B_{lu}J_\nu \varphi _\nu 𝑑\nu )\right\}}{n_uA_{ul}}.$$
(17)
With this definition Eq. 7.2 becomes simply:

$$n_u(n_\mathrm{e}C_{ul}+\rho _{ul}A_{ul})=n_ln_\mathrm{e}C_{lu}.$$
(18)
Let us assume the line is the only contributor to radiation. The energy emitted locally in the line in a small volume of surface 1 cm-2 and length $`dz`$ writes:

$`dF`$ $`=`$ $`h\nu [n_u(A_{ul}+B_{ul}{\displaystyle J_\nu \psi _\nu 𝑑\nu })n_lB_{lu}{\displaystyle J_\nu \varphi _\nu 𝑑\nu }]dz`$ (19)
$`=`$ $`\rho _{ul}n_uAulh\nu dz=\rho _{ul}4\pi ϵdz`$
where $`ϵ`$ is the emissivity coefficient $`n_uA_{ul}h\nu /4\pi `$. Integrating this equation over depth, one gets the emerging line flux:

$$F=4\pi \rho _{ul}\epsilon 𝑑z$$
(20)
According to this equation the divergence flux seems to be equal to the probability that, once emitted at a distance $`z`$ from the surface, a photon can escape from the medium, which is called the “escape probability” $`P_e`$.

From this fact, if one makes the approximation that the emission and the absorption profiles are equal, one computes this escape probability by integrating the attenuation over the line profile:

$$P_e\frac{1}{2}\frac{_0^1\frac{h\nu }{4\pi }n_uA_{ul}\varphi _\nu exp(\tau \varphi _\nu /\mu )𝑑\nu 𝑑\mu }{\frac{h\nu }{4\pi }n_uA_{ul}\varphi _\nu 𝑑\nu }$$
(21)
where $`\mu `$ = cos$`\theta `$, and it gives the usual expressions for complete redistribution in a Voigt profile:

$$P_e(\tau _0)max(\frac{1}{1+2\tau _0\sqrt{\pi ln(\tau _0+1)}},\frac{2}{3}\sqrt{\frac{a}{\tau _0\sqrt{\pi }}}),$$
(22)
where $`\tau _0`$ is the optical thickness at the line center and $`a`$ is the usual damping constant. This expression can thus replace $`\rho _{ul}`$ in Eqs. 18 and 20.

With this method the frequency dependent transfer eq-uations are successfully replaced by frequency integrated ones, and solved consistently with the statistical equat-ions.

Several estimates of the escape probability are used in the literature (see for example Kwan & Krolik, 1981, Rees, Netzer & Ferland, 1989, Kallman & McCray, 1982, Ko & Kallman, 1994), to account for partial redistribution, for line or continuum interlocking, and for continuum absorption. For instance several ways to take into account the continuum opacity (or the overlapping of two lines) have been proposed (Elitzur & Netzer, 1984, Netzer, Elitzur & Ferland, 1985). In the case of the Broad Line Region of AGN, which is optically thin for the continuum underlying the main lines, Collin-Souffrin et al. (1981), Avrett & Loeser (1988), have discussed the influence of different approximations for partial or complete redistribution, and how they compare to an exact transfer treatment, and Collin-Souffrin & Dumont S. (1986) have shown that the escape probability approximation is roughly valid (it can however lead to a decrease by a factor as large as two of the ratio of a lower to a higher transition, such as H$`\alpha `$/H$`\beta `$ or L$`\alpha `$/H$`\alpha `$).

However these approximations are not valid if $`\tau _{\mathrm{cont}}`$ is of the order of unity at the line frequency and if the line photons are created in one place and absorbed in another place where the physical conditions are different. This happens for instance for EUV and soft X-ray lines created in the hot region and absorbed in ionizing He+ ions in a cooler region.

If the local escape probability approximation breaks down, it has a severe consequence: the emitted line spectrum and the thermal and the ionization balance are not correctly computed. For instance the divergence flux takes frequently negative values in our computations, even for intense lines, and even in the reflected spectrum. This is not allowed with escape probabilities. We illustrate this discussion with some examples in Sect. 3.2.


format: text

text:
# The decay 𝜔(782)→5𝜋 in chiral theory.

## Abstract

The arguments are put forward that the many pion decays $`\omega 2\pi ^+2\pi ^{}\pi ^0`$ and $`\pi ^+\pi ^{}3\pi ^0`$ provide an ideal test site for testing the predictions of chiral models of the vector meson decays into many pions. Using the approach based on the Weinberg Lagrangian or, in a new language, the Lagrangian of hidden local symmetry added with the term induced by the anomalous Lagrangian of Wess and Zumino, the partial widths of these decays are evaluated, and their excitation curves in $`e^+e^{}`$ annihilation are obtained. The discussed are the perspectives of the experimental study of the decays $`\omega 5\pi `$ in $`e^+e^{}`$ annihilation and photoproduction.

At present, the unconventional, from the point of view of the chiral pion dynamics, sources of soft pions are feasible. Indeed, the progress in increasing the luminosity of low energy $`e^+e^{}`$ colliders ($`\varphi `$ factories)  could offer the naturally controlled sources of soft pions. Other possible source of such pions could be the intense photon beams  , provided the sufficiently low invariant mass regions of the many pion systems are isolated. The yield of pions is considerably enhanced when they are produced through proper vector resonance states. Then, choosing the many pion decays of sufficiently low lying resonances, one can obtain the soft pions in quantities sufficient for testing the predictions of chiral models that include vector mesons.

The decay $`\omega (782)5\pi `$ whose final state pions possess the momenta $`|𝐪_\pi |74`$ MeV, is just of this kind. The latter value is sufficiently small to expect the manifestation of chiral dynamics in most clean form. By this we mean that the higher derivative and loop terms in the effective Lagrangian are severely suppressed. The present paper is devoted to the evaluation of the partial width of this decay and plotting its excitation curve in $`e^+e^{}`$ annihilation.

The $`\rho \pi `$ sector is considered here on the basis of the Weinberg Lagrangian  revived later as the Lagrangian of hidden local symmetry (HLS) . The former looks as

$``$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(_\mu 𝝆_\nu _\nu 𝝆_\mu +g[𝝆_\mu \times 𝝆_\nu ]\right)^2`$ (3)
$`+{\displaystyle \frac{m_\rho ^2}{2}}\left[𝝆_\mu +{\displaystyle \frac{[𝝅\times _\mu 𝝅]}{2gf_\pi ^2(1+𝝅^2/4f_\pi ^2)}}\right]^2`$
$`+{\displaystyle \frac{(_\mu 𝝅)^2}{2\left(1+𝝅^2/4f_\pi ^2\right)^2}}{\displaystyle \frac{m_\pi ^2𝝅^2}{2(1+𝝅^2/4f_\pi ^2)}},`$
where $`𝝅,m_\pi `$ and $`𝝆,m_\rho `$ stand for the isovector field and mass of pion and $`\rho `$ meson, respectively, cross denotes the vector product in the isovector space, and $`f_\pi =92.4`$ MeV is the pion decay constant. The $`\rho \rho \rho `$ coupling constant $`g`$ and the $`\rho \pi \pi `$ coupling constant $`g_{\rho \pi \pi }`$ are related to the $`\rho `$ mass and pion decay constant $`f_\pi `$ via the parameter of hidden local symmetry $`a`$ as

$`g`$ $`=`$ $`m_\rho /f_\pi \sqrt{a},`$ (4)
$`g_{\rho \pi \pi }`$ $`=`$ $`\sqrt{a}m_\rho /2f_\pi .`$ (5)
Note that $`a=2`$, if one demands the universality condition $`g=g_{\rho \pi \pi }`$ to be satisfied. Then the so called Kawarabayashi-Suzuki-Riazzuddin- Fayyazuddin (KSRF) relation  arises

$$2g_{\rho \pi \pi }^2f_\pi ^2/m_\rho ^2=1$$
(6)
which beautifully agrees with experiment. The $`\rho \pi \pi `$ coupling constant resulting from this relation is $`g_{\rho \pi \pi }=5.89`$. The inclusion of the interaction of the $`\omega (782)`$ with the $`\rho \pi `$ state is achieved upon adding the term induced by the anomalous Lagrangian of Wess and Zumino ,

$$_{\omega \rho \pi }=\frac{N_cg^2}{8\pi ^2f_\pi }\epsilon _{\mu \nu \lambda \sigma }_\mu \omega _\nu \left(𝝅_\lambda 𝝆_\sigma \right),$$
(7)
where $`N_c=3`$ is the number of colors, $`\omega _\nu `$ is the field of $`\omega `$ meson.

One may convince oneself that the $`\omega \rho \pi 5\pi `$ decay amplitude unambiguously results from the Weinberg Lagrangian Eq. (3) and the anomaly induced Lagrangian (7). This amplitude is represented by the diagrams shown in Fig. 1. As one can foresee, its general expression looks cumbersome. However, it can be considerably simplified upon noting that the small pion momenta permit one to use the nonrelativistic expressions,

$`M(\rho ^0\pi _{q_1}^+\pi _{q_2}^+\pi _{q_3}^{}\pi _{q_4}^{})`$ $``$ $`{\displaystyle \frac{g_{\rho \pi \pi }}{2f_\pi ^2}}(\epsilon ,q_1+q_2q_3q_4),`$ (8)
$`M(\rho ^0\pi _{q_1}^+\pi _{q_2}^{}\pi _{q_3}^0\pi _{q_4}^0)`$ $``$ $`{\displaystyle \frac{g_{\rho \pi \pi }}{4f_\pi ^2}}(\epsilon ,q_1q_2),`$ (9)
$`M(\rho ^+\pi _{q_1}^+\pi _{q_2}^+\pi _{q_3}^{}\pi _{q_4}^0)`$ $``$ $`{\displaystyle \frac{g_{\rho \pi \pi }}{4f_\pi ^2}}(\epsilon ,q_1+q_22q_4),`$ (10)
$`M(\rho ^+\pi _{q_1}^+\pi _{q_2}^0\pi _{q_3}^0\pi _{q_4}^0)`$ $``$ $`{\displaystyle \frac{g_{\rho \pi \pi }}{f_\pi ^2}}(\epsilon ,q_1),`$ (11)
for the $`\rho 4\pi `$ decay amplitudes in the diagrams Fig. 1(a), where $`\epsilon `$ stands for the $`\rho `$ meson polarization four-vector The notation for the Lorentz invariant product of two four-vectors $`a`$ and $`b`$ is $`(a,b)=a_0b_o\mathrm{𝐚𝐛}`$. . They are valid with the accuracy $`5\%`$ in the $`4\pi `$ mass range relevant for the present purpose. Likewise, the expression for the combination $`D_\pi ^1M(\pi 3\pi )`$ standing in the expression for the diagrams in Fig. 1(b) can be replaced, with the same accuracy, by $`(8m_\pi ^2)^1`$ times the nonrelativistic $`\pi 3\pi `$ amplitudes. The latter look as

$`M(\pi ^+\pi ^+\pi ^+\pi ^{})`$ $`=`$ $`2m_\pi ^2/f_\pi ^2,`$ (12)
$`M(\pi ^+\pi ^+\pi ^0\pi ^0)`$ $`=`$ $`m_\pi ^2/f_\pi ^2,`$ (13)
$`M(\pi ^0\pi ^+\pi ^{}\pi ^0)`$ $`=`$ $`m_\pi ^2/f_\pi ^2,`$ (14)
$`M(\pi ^0\pi ^0\pi ^0\pi ^0)`$ $`=`$ $`3m_\pi ^2/f_\pi ^2.`$ (15)
Note that, in the nonrelativistic limit, the $`\rho 4\pi `$ decay amplitudes depend on the HLS parameter $`a`$ only through an overall factor $`g_{\rho \pi \pi }/f_\pi ^2=\sqrt{a}m_\rho /2f_\pi ^3,`$ while the $`\pi 3\pi `$ amplitudes do not depend on it at all.

The LHS approach  permits one to include the axial mesons as well §§§The problem of inclusion of vector, axial mesons, and photons in the framework of chiral theories has demanded considerable efforts. See . It is solved in an elegant way in the HLS approach .. An ideal treatment would consist of that under the assumption of $`m_\rho Em_{a_1}`$, the difference between the models with and without $`a_1`$ meson would be reduced to the allowing for the higher derivatives and would be small Allowing for the higher derivatives demands also allowing for the chiral loops, the task which is not yet fulfilled for vector mesons.. In particular, the correction to the $`\rho ^04\pi `$ decay width due to $`a_1`$ meson is estimated at the level of 30-40 $`\%`$ . However, since the invariant mass of four pions in the diagrams Fig. 1(a) is less than 642 MeV, the $`a_1`$ contribution in the $`\omega 5\pi `$ decay is severely suppressed. Furthermore, the introducing $`a_1`$ meson results in the terms $`(_\mu 𝝅)^4`$ in the Lagrangian, but their contribution to the $`\pi 3\pi `$ transition amplitude is inessential in the relevant three pion invariant mass region less than 500 MeV.

Note also that the maximum pion momentum $`|𝐪_\pi |`$ is about 115 MeV, but the contribution of such regions of the phase space is negligible. The dominant contribution comes from the momentum $`|𝐪_\pi |70`$ MeV. In the kinematical situation of the $`\omega 5\pi `$ decay we also neglect the contributions of other higher derivatives and quantum corrections due to chiral loops. Thus, our approach is zeroth order approximation to the $`\omega 5\pi `$ decay amplitude in the framework of the Lagrangians (3) and (7), in a close analogy with the Weinberg amplitudes in classical $`\pi \pi `$ scattering. In this approximation all chiral models of the vector meson interactions with pions indistinguishable, and any difference could manifest upon including the higher derivatives and chiral loops.

Then one obtains, upon neglecting the corrections of the order of $`O(|𝐪_\pi |^4/m_\pi ^4)`$ or higher, the expression for decay amplitudes:

$`M(\omega 2\pi ^+2\pi ^{}\pi ^0)`$ $`=`$ $`{\displaystyle \frac{N_cg_{\rho \pi \pi }g^2}{8(2\pi )^2f_\pi ^3}}\epsilon _{\mu \nu \lambda \sigma }q_\mu ϵ_\nu \{(1+P_{12})q_{1\lambda }[{\displaystyle \frac{(q_2+3q_4)_\sigma }{D_\rho (qq_1)}}{\displaystyle \frac{2q_{4\sigma }}{D_\rho (q_1+q_4)}}]`$ (18)
$`(1+P_{35})q_{3\lambda }\left[{\displaystyle \frac{(q_5+3q_4)_\sigma }{D_\rho (qq_3)}}{\displaystyle \frac{2q_{4\sigma }}{D_\rho (q_3+q_4)}}\right]`$
$`(1+P_{12})(1+P_{35})q_{3\lambda }[{\displaystyle \frac{2q_{4_\sigma }}{D_\rho (qq_4)}}+{\displaystyle \frac{q_{1\sigma }}{D_\rho (q_1+q_3)}}]\},`$
with the final momentum assignment according to $`\pi ^+(q_1)\pi ^+(q_2)\pi ^{}(q_3)\pi ^{}(q_5)\pi ^0(q_4)`$, and

$`M(\omega \pi ^+\pi ^{}3\pi ^0)`$ $`=`$ $`{\displaystyle \frac{N_cg_{\rho \pi \pi }g^2}{8(2\pi )^2f_\pi ^3}}(1P_{12})(1+P_{34}+P_{35})\epsilon _{\mu \nu \lambda \sigma }q_\mu ϵ_\nu q_{1\lambda }`$ (21)
$`\times \{q_{3\sigma }[{\displaystyle \frac{1}{D_\rho (qq_3)}}{\displaystyle \frac{1}{D_\rho (q_1+q_3)}}]`$
$`q_{2\sigma }[{\displaystyle \frac{4}{3D_\rho (qq_1)}}{\displaystyle \frac{1}{2D_\rho (q_1+q_2)}}]\},`$
with the final momentum assignment according to $`\pi ^+(q_1)\pi ^{}(q_2)\pi ^0(q_3)\pi ^0(q_4)\pi ^0(q_5)`$. Note that the first term in each square bracket refers to the specific diagram shown in Fig. 1(a) while the second one does to the diagram shown in Fig. 1(b). In both above formulas, $`ϵ_\nu `$, $`q_\mu `$ stand for four-vectors of polarization and momentum of $`\omega `$ meson, $`P_{ij}`$ is the operator of the interchange the pion momenta $`q_i`$ and $`q_j`$, and

$$D_\rho (q)=m_\rho ^2q^2i\sqrt{q^2}\mathrm{\Gamma }_{\rho \pi \pi }(q^2)$$
(22)
is the inverse propagator of $`\rho `$ meson whose width is dominated by the $`\pi \pi `$ decay mode:

$$\mathrm{\Gamma }_{\rho \pi \pi }(q^2)=\mathrm{\Gamma }_{\rho \pi \pi }(m_\rho ^2)\frac{m_\rho ^2}{q^2}\left(\frac{q^24m_\pi ^2}{m_\rho ^24m_\pi ^2}\right)^{3/2}.$$
(23)
Yet even in this simplified form the expressions for the $`\omega 5\pi `$ amplitudes are not easy to use for evaluation of the branching ratios. To go further, it should be noted the following. One can check that the invariant mass of the $`4\pi `$ system on which the contribution of the diagrams shown in Fig. 1(a) depends, changes in very narrow range $`558\text{ MeV}<m_{4\pi }<642\text{ MeV}`$. Hence, one can set it in all the $`\rho `$ propagators standing as the first terms in all square brackets in Eqs. (21) and (18), with the accuracy $`20\%`$ in width, to the equilibrium value $`\overline{m_{4\pi }^2}^{1/2}=620`$ MeV evaluated for the pion energy $`E_\pi =m_\omega /5`$ which gives the dominant contribution. The same is true for the invariant mass of the pion pairs on which the $`\rho `$ propagators standing as the last terms in square brackets of the above expressions depend. This invariant mass varies in the narrow range $`280\text{ MeV}<m_{2\pi }<360\text{ MeV}`$. With the same accuracy, one can set it to $`\overline{m_{2\pi }^2}^{1/2}=295`$ MeV in all relevant propagators. On the other hand, the amplitude of the process $`\omega \rho ^0\pi ^0(2\pi ^+2\pi ^{})\pi ^0`$ is

$`M[\omega \rho ^0\pi ^0(2\pi ^+2\pi ^{})\pi ^0]`$ $`=`$ $`4{\displaystyle \frac{N_cg_{\rho \pi \pi }g^2}{8(2\pi )^2f_\pi ^3}}`$ (26)
$`\times \epsilon _{\mu \nu \lambda \sigma }q_\mu ϵ_\nu (q_1+q_2)_\lambda `$
$`\times {\displaystyle \frac{q_{4\sigma }}{D_\rho (qq_4)}},`$
where the momentum assignment is the same as in Eq. (18). The other relevant amplitude corresponding to the first diagram in Fig. 1(b) is

$`M[\omega \rho ^0\pi ^0(\pi ^+\pi ^{})(\pi ^+\pi ^{}\pi ^0)]`$ $`=`$ $`{\displaystyle \frac{N_cg_{\rho \pi \pi }g^2}{8(2\pi )^2f_\pi ^3}}\epsilon _{\mu \nu \lambda \sigma }q_\mu ϵ_\nu `$ (28)
$`\times (1+P_{12})(1+P_{35}){\displaystyle \frac{q_{1\lambda }q_{3\sigma }}{D_\rho (q_1+q_3)}}`$
Then, taking into account the above consideration concerning the invariant masses, and comparing Eqs. (18), (26), and (28), one can see that

$`M(\omega 2\pi ^+2\pi ^{}\pi ^0)`$ $``$ $`{\displaystyle \frac{5}{2}}M[\omega \rho ^0\pi ^0(2\pi ^+2\pi ^{})\pi ^0]`$ (30)
$`\times \left[1{\displaystyle \frac{D_\rho (\overline{m_{4\pi }^2})}{2D_\rho (\overline{m_{2\pi }^2})}}\right].`$
The same treatment shows that

$`M(\omega \pi ^+\pi ^{}3\pi ^0)`$ $``$ $`{\displaystyle \frac{5}{2}}M[\omega \rho ^+\pi ^{}(\pi ^+3\pi ^0)\pi ^{}]`$ (32)
$`\times \left[1{\displaystyle \frac{D_\rho (\overline{m_{4\pi }^2})}{2D_\rho (\overline{m_{2\pi }^2})}}\right],`$
where

$`M[\omega \rho ^+\pi ^{}(\pi ^+3\pi ^0)\pi ^{}]`$ $`=`$ $`4{\displaystyle \frac{N_cg_{\rho \pi \pi }g^2}{8(2\pi )^2f_\pi ^3}}`$ (34)
$`\times {\displaystyle \frac{\epsilon _{\mu \nu \lambda \sigma }q_\mu ϵ_\nu q_{1\lambda }q_{2\sigma }}{D_\rho (qq_2)}},`$
and the final momenta assignment is the same as in Eq. (21). The numerical values of $`\overline{m_{4\pi }^2}^{1/2}`$ and $`\overline{m_{2\pi }^2}^{1/2}`$ found above are such that the correction factor in parentheses of Eqs. (28) and (30) amounts to $`20\%`$ in magnitude. In what follows, the above correction will be taken into account as an overall factor of 0.64 in front of the branching ratios of the decays $`\omega 5\pi `$. When making this estimate, the imaginary part of the $`\rho `$ propagators in square brackets of Eq. (28) and (30) is neglected. This assumption is valid with the accuracy better than $`1\%`$ in width.

It would be useful to fulfill the model estimate of the $`\omega 5\pi `$ partial widths as follows. The corresponding equilibrium pion momenta are $`|𝐪_{\pi ^+}|=70`$ MeV and $`|𝐪_{\pi ^0}|=79`$ MeV. The integrations over angles of final pions can be fulfilled assuming them independent. Using the nonrelativistic expression for phase space of five pions ,

$$R_5=\frac{\pi ^6(\underset{i}{\overset{5}{}}m_{\pi i})^{1/2}}{60(_i^5m_{\pi i})^{3/2}}(m_\omega \underset{i}{\overset{5}{}}m_{\pi i})^5$$
(35)
whose numerical value coincides with the accuracy $`1\%`$ with the numerically evaluated exact expression, and introducing the branching ratio at the $`\omega `$ mass as

$$B_{\omega 5\pi }=\mathrm{\Gamma }_{\omega 5\pi }(m_\omega )/\mathrm{\Gamma }_\omega ,$$
(36)
one finds

$`B(\omega 5\pi )`$ $``$ $`\left[{\displaystyle \frac{5N_c}{2\pi ^2}}\left({\displaystyle \frac{g_{\rho \pi \pi }g^2}{8f_\pi ^3}}\right){\displaystyle \frac{m_\omega |𝐪_{\pi ^+}|}{|D_\rho (\overline{m_{4\pi }^2})|}}\right]^2`$ (39)
$`\times {\displaystyle \frac{R_5}{18(2\pi )^{11}m_\omega \mathrm{\Gamma }_\omega }}\left|1{\displaystyle \frac{D_\rho (\overline{m_{4\pi }^2})}{2D_\rho (\overline{m_{2\pi }^2})}}\right|^2`$
$`\times \{{\displaystyle \genfrac{}{}{0pt}{}{|𝐪_{\pi ^0}|^2(2\pi ^+2\pi ^{}\pi ^0)}{\frac{|𝐪_{\pi ^+}|^2}{3}(\pi ^+\pi ^{}3\pi ^0)}}`$
The calculation gives $`B(\omega 2\pi ^+2\pi ^{}\pi ^0)=2.5\times 10^9`$ and $`B(\omega \pi ^+\pi ^{}3\pi ^0)=1.0\times 10^9`$.

The evaluation of the partial widths valid with accuracy $`20\%`$ can be obtained upon taking the amplitude of each considered decays as 5/2 times the $`\rho \pi `$ production state amplitude with the subsequent decay $`\rho 4\pi `$, and calculate the partial width using the calculated widths of the latter:

$`B_{\omega 2\pi ^+2\pi ^{}\pi ^0}`$ $`=`$ $`\left|1{\displaystyle \frac{D_\rho (\overline{m_{4\pi }^2})}{2D_\rho (\overline{m_{2\pi }^2})}}\right|^2\left({\displaystyle \frac{5}{2}}\right)^2{\displaystyle \frac{2}{\pi \mathrm{\Gamma }_\omega }}`$ (43)
$`\times {\displaystyle _{4m_{\pi ^+}}^{m_\omega m_{\pi ^0}}}dm`$
$`\times {\displaystyle \frac{m^2\mathrm{\Gamma }_{\omega \rho ^0\pi ^0}(m)\mathrm{\Gamma }_{\rho 2\pi ^+2\pi ^{}}(m)}{|D_\rho (m^2)|^2}}`$
$`=1.1\times 10^9`$
where

$$\mathrm{\Gamma }_{\omega \rho ^0\pi ^0}(m)=g_{\omega \rho \pi }^2q^3(m_\omega ,m,m_{\pi ^0})/12\pi ,$$

$$g_{\omega \rho \pi }=\frac{N_cg^2}{8\pi ^2f_\pi }=14.3\text{ GeV}^1.$$

Note also the $`a^1`$ dependence of the $`\omega 5\pi `$ width on the HLS parameter $`a`$. The branching ratio $`B_{\omega \pi ^+\pi ^{}3\pi ^0}`$ is obtained from Eq. (43) upon inserting the lower integration limit to $`m_{\pi ^+}+3m_{\pi ^0}`$, $`m_{\pi ^0}m_{\pi ^+}`$ in the expression for the momentum $`q`$ and substitution of the $`\rho ^+\pi ^+3\pi ^0`$ decay width corrected for the mass difference of charged and neutral pions. Of course, the main correction of this sort comes from the phase space volume of the final 4$`\pi `$ state. One obtains

$`B_{\omega \pi ^+\pi ^{}3\pi ^0}`$ $`=`$ $`\left|1{\displaystyle \frac{D_\rho (\overline{m_{4\pi }^2})}{2D_\rho (\overline{m_{2\pi }^2})}}\right|^2\left({\displaystyle \frac{5}{2}}\right)^2{\displaystyle \frac{2}{\pi \mathrm{\Gamma }_\omega }}`$ (47)
$`\times {\displaystyle _{m_{\pi ^+}+3m_{\pi ^0}}^{m_\omega m_{\pi ^+}}}dm`$
$`\times {\displaystyle \frac{m^2\mathrm{\Gamma }_{\omega \rho ^+\pi ^{}}(m)\mathrm{\Gamma }_{\rho ^+\pi ^+3\pi ^0}(m)}{|D_\rho (m^2)|^2}}`$
$`=8.5\times 10^{10}`$
where

$$\mathrm{\Gamma }_{\omega \rho ^+\pi ^{}}(m)=g_{\omega \rho \pi }^2q^3(m_\omega ,m,m_{\pi ^+})/12\pi .$$

As is pointed out in Ref. , the inclusion of the direct $`\omega \pi ^+\pi ^{}\pi ^0`$ vertex reduces the 3$`\pi `$ decay width of the $`\omega `$ by $`33\%`$. This implies that one should make the following replacement to take into account the effect of the pointlike diagrams in Fig. 1(b) in the expression for the suppression factor:

$`\left|1{\displaystyle \frac{D_\rho (\overline{m_{4\pi }^2})}{2D_\rho (\overline{m_{2\pi }^2})}}\right|^2`$ $``$ $`|1{\displaystyle \frac{D_\rho (\overline{m_{4\pi }^2})}{2}}[{\displaystyle \frac{1}{D_\rho (\overline{m_{2\pi }^2})}}`$ (50)
$`{\displaystyle \frac{1}{3m_\rho ^2}}\left]\right|^2|1{\displaystyle \frac{D_\rho (\overline{m_{4\pi }^2})}{3D_\rho (\overline{m_{2\pi }^2})}}|^2`$
$`0.75,`$
instead of 0.64, which results in the increase of the above branching ratios by the factor of 1.17.

The numerical value of the $`\omega 5\pi `$ decay width changes by the factor of two when varying the energy within $`\pm \mathrm{\Gamma }_\omega /2`$ around the $`\omega `$ mass. In other words, the dependence of this partial width on energy is very strong. This is illustrated by Fig. 2 where the $`\omega 5\pi `$ excitation curves in $`e^+e^{}`$ annihilation,

$`\sigma _{e^+e^{}\omega 5\pi }(s)`$ $`=`$ $`12\pi \left({\displaystyle \frac{m_\omega }{E}}\right)^3\mathrm{\Gamma }_{\omega e^+e^{}}(m_\omega )`$ (52)
$`\times {\displaystyle \frac{\mathrm{\Gamma }_\omega B_{\omega 5\pi }(E)}{\left[(sm_\omega ^2)^2+(m_\omega \mathrm{\Gamma }_\omega )^2\right]}},`$
are plotted. Here $`B_{\omega 2\pi ^+2\pi ^{}\pi ^0}(E)`$ \[$`B_{\omega \pi ^+\pi ^{}3\pi ^0}(E)`$\] is given by Eq. (43) \[(47)\], respectively, with the substitution $`m_\omega E`$. The mentioned strong energy dependence of the partial width results in the asymmetric shape of the $`\omega `$ resonance and the shift of its peak by $`+0.7`$ MeV . As is seen from Fig. 2, the peak value of the $`5\pi `$ state production cross section is about 1.5-2.0 femtobarns. Yet the decays $`\omega 5\pi `$ can be observable on $`e^+e^{}`$ colliders. Indeed, with the luminosity $`L=10^{33}\text{cm}^2\text{s}^1`$ near the $`\omega `$ peak, which seems to be feasible, one may expect about 2 events per week for the considered decays to be detected at these colliders.

As for the angular distributions of the final pions are concerned, they, of course, should be deduced from the full amplitudes Eqs. (18) and (21). However, some qualitative conclusions about the angular distributions can be drawn from the simplified expressions Eqs. (26), (30). Since helicity is conserved, only the states of the $`\omega (782)`$ with the spin projections $`\lambda =\pm 1`$ on the $`e^+e^{}`$ beam axes given by the unit vector $`𝐧_0`$ are populated. In what follows the suitable notation for the vector product of the pion momenta are used:

$`[𝐪_i\times 𝐪_j]`$ $`=`$ $`|𝐪_i||𝐪_j|\mathrm{sin}\theta _{ij}`$ (54)
$`\times (\mathrm{sin}\mathrm{\Theta }_{ij}\mathrm{cos}\mathrm{\Phi }_{ij},\mathrm{sin}\mathrm{\Theta }_{ij}\mathrm{sin}\mathrm{\Phi }_{ij},\mathrm{cos}\mathrm{\Theta }_{ij}).`$
In other words, $`\theta _{ij}`$ is the angle between the pion momenta $`𝐪_i`$ and $`𝐪_j`$, $`\mathrm{\Theta }_{ij}`$, $`\mathrm{\Phi }_{ij}`$ being the polar and azimuthal angles of the normal to the plane spanned by the momenta $`𝐪_i`$ and $`𝐪_j`$. Choosing $`𝐧_0`$ to be the unit vector along z axes, the probability density of the emission of two $`\pi ^+`$’s with the momenta $`𝐪_1`$, $`𝐪_2`$, and $`\pi ^0`$ with the momentum $`𝐪_4`$ is represented as

$`w`$ $``$ $`\left[𝐪_4\times \left(𝐪_1+𝐪_2\right)\right]^2\left(𝐧_0\left[𝐪_4\times \left(𝐪_1+𝐪_2\right)\right]\right)^2`$ (58)
$`=𝐪_4^2[𝐪_1^2\mathrm{sin}^2\theta _{41}\mathrm{sin}^2\mathrm{\Theta }_{41}+𝐪_2^2\mathrm{sin}^2\theta _{42}\mathrm{sin}^2\mathrm{\Theta }_{42}`$
$`+2|𝐪_1||𝐪_2|\mathrm{sin}\mathrm{\Theta }_{41}\mathrm{sin}\mathrm{\Theta }_{42}\mathrm{sin}\theta _{41}\mathrm{sin}\theta _{42}`$
$`\times \mathrm{cos}(\mathrm{\Phi }_{41}\mathrm{\Phi }_{42})]`$
in the case of the final state $`2\pi ^+2\pi ^{}\pi ^0`$. Here the momentum assignment is the same as in Eq. (18). The angular distribution of two $`\pi ^{}`$’s with the momenta $`𝐪_3`$, $`𝐪_5`$, and $`\pi ^0`$ is obtained from Eq. (58) upon the replacement $`𝐪_{1,2}𝐪_{3,5}`$, because the identity $`\epsilon _{\mu \nu \lambda \sigma }q_\mu ϵ_\nu (q_1+q_2)_\lambda q_{4\sigma }=\epsilon _{\mu \nu \lambda \sigma }q_\mu ϵ_\nu (q_3+q_5)_\lambda q_{4\sigma }`$ is valid. Since another identity $`\epsilon _{\mu \nu \lambda \sigma }q_\mu ϵ_\nu (q_1+q_2)_\lambda q_{4\sigma }=\epsilon _{\mu \nu \lambda \sigma }q_\mu ϵ_\nu (q_1+q_2)_\lambda (q_3+q_5)_\sigma `$ is valid, one can write the angular distribution that includes four charged pions:

$`w`$ $``$ $`[(𝐪_1+𝐪_2)\times (𝐪_3+𝐪_5)]^2\left(𝐧_0[(𝐪_1+𝐪_2)\times (𝐪_3+𝐪_5)]\right)^2`$ (63)
$`=(1+P_{12})(1+P_{35})𝐪_1^2𝐪_3^2\mathrm{sin}^2\theta _{13}\mathrm{sin}^2\mathrm{\Theta }_{13}`$
$`+2|𝐪_1||𝐪_2|(1+P_{35})𝐪_3^2\mathrm{sin}\theta _{13}\mathrm{sin}\theta _{23}\mathrm{sin}\mathrm{\Theta }_{13}\mathrm{sin}\mathrm{\Theta }_{23}\mathrm{cos}(\mathrm{\Phi }_{13}\mathrm{\Phi }_{23})`$
$`+2|𝐪_3||𝐪_5|(1+P_{12})𝐪_1^2\mathrm{sin}\theta _{13}\mathrm{sin}\theta _{15}\mathrm{sin}\mathrm{\Theta }_{13}\mathrm{sin}\mathrm{\Theta }_{15}\mathrm{cos}(\mathrm{\Phi }_{13}\mathrm{\Phi }_{15})`$
$`+2|𝐪_1||𝐪_2||𝐪_3||𝐪_5|(1+P_{35})\mathrm{sin}\theta _{13}\mathrm{sin}\theta _{25}\mathrm{sin}\mathrm{\Theta }_{13}\mathrm{sin}\mathrm{\Theta }_{25}\mathrm{cos}(\mathrm{\Phi }_{13}\mathrm{\Phi }_{25}).`$
Here $`P_{ij}`$ interchanges the indices $`i`$ and $`j`$. In the case of the final state $`\pi ^+\pi ^{}3\pi ^0`$ the corresponding probability density can be obtained from Eqs. (32) and (34) and looks as

$`w`$ $``$ $`[𝐪_1\times 𝐪_2]^2(𝐧_0[𝐪_1\times 𝐪_2])^2`$ (65)
$`=𝐪_1^2𝐪_2^2\mathrm{sin}^2\theta _{21}\mathrm{sin}^2\mathrm{\Theta }_{21}.`$
Here the momentum assignment is the same as in Eq. (21). The corresponding angular distribution of one charged, say $`\pi ^+`$, and three neutral pions can be obtained from Eqs. (32) and (34) upon using the identity $`\epsilon _{\mu \nu \lambda \sigma }q_\mu ϵ_\nu q_{1\lambda }q_{2\sigma }=\epsilon _{\mu \nu \lambda \sigma }q_\mu ϵ_\nu q_{1\lambda }(q_3+q_4+q_5)_\sigma `$ and looks as

$`w`$ $``$ $`\left[𝐪_1\times {\displaystyle \underset{i}{}}𝐪_i\right]^2\left(𝐧_0\left[𝐪_1\times {\displaystyle \underset{i}{}}𝐪_i\right]\right)^2`$ (69)
$`=𝐪_1^2[{\displaystyle \underset{i}{}}𝐪_i^2\mathrm{sin}^2\theta _{i1}\mathrm{sin}^2\mathrm{\Theta }_{i1}`$
$`+2{\displaystyle \underset{ij}{}}|𝐪_i||𝐪_j|\mathrm{sin}\theta _{i1}\mathrm{sin}\theta _{j1}\mathrm{sin}\mathrm{\Theta }_{i1}\mathrm{sin}\mathrm{\Theta }_{j1}`$
$`\times \mathrm{cos}(\mathrm{\Phi }_{i1}\mathrm{\Phi }_{j1})].`$
Here indeces $`i,j`$ run over 3,4,5.

The strong energy dependence of the five pion partial width of the $`\omega `$ implies that the branching ratio at the $`\omega `$ mass, Eq. (36), evaluated above, is slightly different from that determined by the expression

$$B_{\omega 5\pi }^{\mathrm{aver}}(E_1,E_2)=\frac{2}{\pi }_{E_1}^{E_2}𝑑E\frac{E^2\mathrm{\Gamma }_\omega B_{\omega 5\pi }(E)}{(E^2m_\omega ^2)^2+(m_\omega \mathrm{\Gamma }_\omega )^2}.$$
(70)
Taking $`E_1=`$ 772 MeV and $`E_2=`$ 792 MeV, one finds $`B_{\omega 2\pi ^+2\pi ^{}\pi ^0}^{\mathrm{aver}}(E_1,E_2)=9.0\times 10^{10}`$ and $`B_{\omega \pi ^+\pi ^{}3\pi ^0}^{\mathrm{aver}}(E_1,E_2)=6.7\times 10^{10}`$ to be compared to Eq. (43) and (47), respectively. In particular, the quantity $`B_{\omega 2\pi ^+2\pi ^{}\pi ^0}^{\mathrm{aver}}(E_1,E_2)`$ is the relevant characteristics of this specific decay mode in photoproduction experiments. The Jefferson Lab ”photon factory”  could also be suitable for detecting the five pion decays of the $`\omega `$. However, in view of the suppression of the $`\omega `$ photoproduction cross section by the factor of 1/9 as compared with the $`\rho `$ one, the total number of $`\omega `$ mesons will amount to $`7\times 10^8`$ per nucleon. Hence, the increase of intensity of this machine by the factor of 50 is highly desirable, in order to observe the decay $`\omega 5\pi `$ and measure its branching ratio. Evidently, the $`\omega `$ photoproduction on heavy nuclei is preferable in view of the dependence of the cross section on atomic weight $`A`$ growing as $`A^{0.80.95}`$ .

The conclusions about the angular distributions in photoproduction are the following. Of course, their general expression should be deduced from the full decay amplitudes Eqs. (18) and (21), together with the detailed form of the photoproduction mechanism. The qualitative picture, however, can be obtained upon noting that $`s`$-channel helicity conservation is a good selection rule for the photoproduction reactions. Then in the helicity reference frame characterized as the frame where the $`\omega `$ is at rest, while its spin quantization axes is directed along the $`\omega `$ momentum in the center-of-mass system, the expressions for the angular distributions coincide with the expressions Eqs. (58),(63), (65) and (69). Since, at high energies, the direction of the final $`\omega `$ momentum lies at the scattering angle less than $`0.5^{}`$ in the case of the photoproduction on heavy nuclei, the vector $`𝐧_0`$ can be treated as pointed along the photon beam direction.

Together with the $`e^+e^{}`$ annihilation experiments, the study of the photoproduction of the five pion states on heavy nuclei would also allow to measure the corresponding partial width of the $`\omega (782)`$. The comparison with theoretical expectations presented here would give the possibility of testing the predictions of chiral models that include the vector mesons in the situation where the decay amplitudes are determined by very low pion momenta.

We are grateful to G. N. Shestakov for discussion. The present work is supported in part by the grant RFBR-INTAS IR-97-232.