text:
# Charmonium production in relativistic proton-nucleus collisions : What will we learn from the negative 𝑥_𝐹 region ?

\[

## Abstract

We study the nuclear medium effects on the $`c\overline{c}`$ time evolution and charmonium production, in a relativistic proton-nucleus collision. In particular, we focus on the fragmentation region of the nucleus where the formation length of the charmonium mesons is shorter than the size of the nucleus. Little is known on the nuclear effects in this region. We use a quantum-mechanical model which includes a realistic potential for the $`c\overline{c}`$ system and an imaginary potential to describe the collisions of the $`c\overline{c}`$ with the nucleons. The imaginary potential introduces a transition amplitude among the charmonium states and produces an interference pattern on the charmonium survival probability, which is particulartly important for $`\psi ^{}`$. Our results on the suppression factors are compared with data from the NA50 and E866/NuSea Collaborations. Predictions are given for the suppression of $`J/\psi ,\psi ^{},\chi _c`$ as a function of the nuclear mass and in the negative $`x_F`$ region, where data will be available soon.

\]

A very exciting era in the study of strong interactions has begun with the advent of high energy heavy ion collision experiments. In particular, this has opened the way to the exploration of new states of matter, such as the colour glass condensate, corresponding to the nonlinear regime of quantum chromodynamics, or the high temperature phase of nuclear matter, the so-called quark-gluon plasma, which is also important for our understanding of the early universe. In relativistic heavy ion collisions, evidence for the formation of the quark-gluon plasma may come from the combined information of different signals, such as strangeness, dilepton, photon or charmonium production.

Matsui and Satz  first suggested that the suppression of charmonia could signal that a phase transition has occurred. This pioneering idea has triggered a series of experiments at the CERN SPS , whose interpretation has been very controversial. It has become clear that, to use charmonium as a signal, all possible mechanisms which can affect charmonium production need to be well understood .

In this context, proton-proton (pp) and proton-nucleus (pA) collisions have been intensively studied . In fact, these processes are used as a reference for the case of nucleus-nucleus (AB) collisions where a critical energy density may be eventually attained, producing the plasma phase. This should affect various observables and in particular, it may lead to an extra (anomalous) charmonium suppression. There have been indications recently that the plasma has been produced, but it is still too early to draw definite conclusions . These studies will be pursued by the running experiments at RHIC, and future measurements at LHC.

It takes a certain time for a $`c\overline{c}`$ produced in a collision to become a colour singlet state, then it expands to the size of a $`\psi `$ meson (the symbol $`\psi `$ stands from now on for any charmonium meson, i.e. $`J/\psi ,\psi ^{}`$ and $`\chi _c`$). Different models exist at present for the production mechanisms .

Sufficiently fast $`c\overline{c}`$ can traverse the nucleus before a $`\psi `$ meson is fully formed; whereas slow $`c\overline{c}`$ can traverse the nucleus as fully formed $`\psi `$ meson . Experimentally one can pinpoint to these different kinematical regions by measuring charmonium production at positive and negative small Feynman $`x_F`$ and/or making inverse kinematics measurements .

One of the aspects which need to be well understood is the effect of the nuclear medium on the $`c\overline{c}`$ pair evolution. In particular, very little is known on the nuclear medium effects when the formation time of a fully developed charmonium is smaller than the time to traverse the nucleus . Soon, measurements in $`pA`$ and $`AB`$ collisions will be performed by the HERA-B Collaboration  at DESY and the PHENIX Collaboration at RHIC . Moreover the E160 Collaboration at SLAC will perform similar experiments where the charmonia are photoproduced .

In this letter, we focus on this region, i.e. the fragmentation region of the nucleus (negative $`x_F`$ region or slow $`c\overline{c}`$) and study the effects of the nuclear medium on the $`c\overline{c}`$ time evolution and on charmonium suppression in relativistic pA collisions . Because of the kinematical region we are interested in, we assume that the produced $`c\overline{c}`$ becomes quickly a colour singlet state (premeson) and on its way through the nucleus, this premeson expands to a $`\psi `$ meson, while experiencing collisions with the nucleons. To describe this, we use here a quantum-mechanical model where, contrary to previous works , the $`c\overline{c}`$ pair is bound by a realistic potential  and the premeson (and then meson) wave function is expanded on the basis given by the charmonium eigenstates. An imaginary potential depending on the dipole charmonium-nucleon cross section  is included to describe the collisions with the nucleons. This introduces a transition amplitude among the charmonium states. (We assume the premeson wavefunction has a fixed angular momentum quantum number, while it has different radial components.)

We will show that the interaction of the premeson and later on the meson with the nuclear medium produces an interference pattern on the charmonium survival probability. We will show how this affects the charmonium suppression both when the path in the nucleus is varied, by changing the mass number $`A`$, and as a function of $`x_F`$. These effects are particularly important for $`\psi ^{}`$. We compare our results with experimental data obtained at CERN SPS by the Na50  and at Fermilab by the E866/NuSea  Collaborations. We present predictions for the production of $`J/\psi ,\psi ^{},\chi _c`$ in a range of small negative $`x_F`$ values, for different nuclei , for which experimental data will be avalaible soon .

Our $`c\overline{c}`$ pair, produced in a $`pA`$ collision, is described (in its center of mass frame) by the internal Hamiltonian

$$H_0=\frac{p^2}{m_c}+V(r),$$
(1)
where $`m_c`$ is the charm quark mass. The potential $`V(r)`$ which binds $`c`$ and $`\overline{c}`$ is a realistic potential which reproduces well the properties of the charmonium family . The spin-dependent terms are neglected. Only transitions between charmonium states with different $`n`$ quantum numbers are considered.

By solving the static Schrödinger equation with $`H_0`$ Eq.(1), one gets the eigenenergies $`E_{n,\mathrm{}}`$ and the corresponding wave functions $`|n,\mathrm{}`$ (to simplify the notations, we do not write explicitly the dependence on magnetic quantum number $`m`$). The 1S, 2S states correspond in our model to $`J/\psi ,\psi ^{}`$ respectively and 1P to $`\chi _c`$ where the different total spin state of $`\chi _c`$ are not split in energy. (We consider that our $`\chi _c`$ is produced through a two-gluon mechanism).

The time dependent wave function for the $`c\overline{c}`$ in its rest frame is expanded on the basis of eigenstates of $`H_0`$ Eq.(1):

$$|c\overline{c},\mathrm{}(\tau )=\underset{n=0}{\overset{\mathrm{}}{}}c_{n,\mathrm{}}(\tau )e^{iE_{n,\mathrm{}}\tau }|n,\mathrm{}.$$
(2)
($`\mathrm{},c=1`$). In practice, one truncates the sum to a number $`\overline{n}`$ of eigenstates, for a given $`\mathrm{}`$ value. As we will discuss, we have chosen $`\overline{n}`$ large enough so that the results we present are not sensitive to this truncation.

To model the interaction of the $`c\overline{c}`$ with the nuclear medium we add an imaginary part to Eq.(1) :

$$iW=i\frac{\gamma v}{2}\sigma (\stackrel{}{r}_T,\sqrt{s_{\psi N}})\rho (\stackrel{}{b},z)$$
(3)
where $`v`$ is the speed of the nucleus with respect to the $`c\overline{c}`$ frame and $`\sigma `$ is the dipole cross section associated to the interaction of the $`c\overline{c}`$ with a nucleon $`N`$, $`r_T`$ is the transverse distance between $`c`$ and $`\overline{c}`$, $`\sqrt{s_{\psi N}}`$ is the energy in the center of mass of the $`\psi N`$ system and $`\rho `$ is the nuclear density evaluated at the position $`(\stackrel{}{b},z)`$ of the $`c\overline{c}`$ center of mass, $`z`$ being the beam direction. For the dipole cross section we use the parametrization $`\sigma (r_T,\sqrt{s})=\sigma _0(s)(1e^{r_T^2/r_0^2(s)})`$, determined by fitting deep inelastic scattering data and which well reproduces the charmonium photoproduction data . Concerning the nuclear density $`\rho `$, we present results obtained with a Woods-Saxon profile, with parameters chosen such as to well reproduce the nuclear radii.

The time evolution of $`c\overline{c}`$ wavefunction (2) in the nucleus is determined by solving the time dependent Schrödinger equation for $`H=H_0iW`$ Eqs.(1) and (3). This leads to a system of first-order coupled-channel differential equations for the amplitudes $`c_{n,l}(t)`$ :

$$\dot{c}_{n,\mathrm{}}(t)=a\underset{k=1}{\overset{\overline{n}}{}}c_{k,\mathrm{}}(t)e^{i(E_{n,\mathrm{}}E_{k,\mathrm{}})\frac{t}{\gamma }}n,\mathrm{}|\sigma |k,\mathrm{}$$
(4)
with $`a=v\rho (\stackrel{}{b},z)/2`$. The time $`t=\gamma \tau `$ is now in the laboratory frame. We have neglected in this first calculation the higher Fock states that should emerge from the Lorentz boost . We see from (4) that the imaginary potential (3) introduces a transition amplitude among the charmonium eigenstates.

A difficult choice is that of the initial conditions $`c_{n,\mathrm{}}(0)`$ for (4), related to the mechanism of hadroproduction of charmonia, which is still badly known . We use:

$$\varphi _{c\overline{c},\mathrm{}}(0)=c_{\mathrm{}}f_{\mathrm{}}(r_T)e^{\frac{1}{2}\beta ^2(r_T^2+z^2)}$$
(5)
where $`c_{\mathrm{}}`$ is a normalization constant. We describe the conversion of a gluon into a $`c\overline{c}`$ through a gaussian multiplied by the function $`f_{\mathrm{}}(r_T)=r_T^2`$ (or $`\stackrel{}{r}_T`$) to account for the two (one) supplementary gluons necessary to produce $`J/\psi ,\psi ^{}`$ ($`\chi _c`$). The initial wave function depends on one parameter only, $`\beta `$. This parameter has been determined by fixing the ratio $`|c_{2,0}/c_{1,0}|^2`$ to the experimental ratio of $`\psi ^{}`$ over $`J/\psi `$, produced in $`pp`$ collisions at 450 GeV, i.e. $`B_{\psi ^{}\mu \mu }\sigma ^{pp\psi ^{}}/B_{J/\psi \mu \mu }\sigma ^{ppJ/\psi }=(1.60\pm 0.04)\%`$ . where $`B_{\psi \mu \mu }`$ is the branching ratio to dimuon production and $`\sigma ^{pp\psi }`$ the $`pp`$ reaction cross section. This leads to the ratio $`|c_{2,0}/c_{1,0}|^2(0)=0.22\pm 0.03`$ and $`\beta =1.33\pm 0.5`$ GeV (the same $`\beta `$ has been used for the $`\chi _c`$ states).

The number $`\overline{n}`$ of states in (2), to be included in the coupled channel equations (4), has to be large enough that the results do not depend on the truncation. We have checked how the results depend on the inclusion of extra states. Up to 4 channels have been used for the S-states and up to 2 for the P-states. The truncation always affects a little the last state included; whereas the results for the lowest energy states are practically unchanged.

Let us now come to the results, obtained in infinite nuclear matter first, i.e. by taking in (4) a constant nuclear density $`\rho =\rho _0`$. In fig.1 we show the time evolution of the probabilities associated to $`J/\psi `$ and $`\psi ^{}`$ production obtained by solving (4) with the initial conditions (5), both neglecting and including the transition amplitudes between different eigenstates (non-diagonal terms). While in the former case the probabilities decrease exponentially, in the latter they present an oscillation. The interference pattern is more pronounced in the case of $`\psi ^{}`$ whereas for $`J/\psi `$ the oscillations stays very close to the exponential. This effect directly influences the suppression factors in $`pA`$ collisions as we will see. The dominant oscillation frequency in Eq.(4) is $`\omega =(E_{2,\mathrm{}}E_{1,\mathrm{}})/\gamma `$. We have also seen that the deviation from the uncoupled solution becomes stronger the higher $`\gamma `$ is.

Let us now consider a $`c\overline{c}`$ which is produced in a relativistic $`pA`$ collision and evolves according to Eqs.(4). Integrating over all possible paths the ratio $`|c_\psi (\mathrm{})/c_\psi (0)|^2`$, which gives the survival probability of the $`\psi `$ in the nucleus, one gets the $`pA\psi `$ reaction cross section :

$$\sigma ^{pA\psi }=𝑑\stackrel{}{b}𝑑z\rho (\stackrel{}{b},z)\sigma ^{pN\psi }\left|\frac{c_\psi (t(\stackrel{}{b},z))}{c_\psi (0)}\right|^2$$
(6)
where $`t(\stackrel{}{b},z)`$ is the time necessary to a $`c\overline{c}`$ produced at a point $`(\stackrel{}{b},z)`$ to traverse the nucleus. If $`\sigma ^{pA\psi }=A\sigma ^{pN\psi }`$, the suppression factor, defined as

$$S_A^\psi \frac{\sigma ^{pA\psi }}{A\sigma ^{pN\psi }},$$
(7)
becomes $`S_A^\psi =1`$.

Before making predictions, let us compare our calculations to existing data. For $`J/\psi `$, we have to take into account that the dimuons measured in the detector come from the decay of the directly produced $`J/\psi `$’s as well as those produced by the decays of $`\psi ^{}`$’s and $`\chi _c`$’s. We will denote by $`\text{}J/\psi \text{}`$ the total number of produced $`J/\psi `$, i.e. $`S_{pA}^{\text{}J/\psi \text{}}=0.62S_{pA}^{J/\psi }+0.30S_{pA}^{\chi _c}+0.08S_{pA}^\psi ^{}`$  .

Let us first discuss charmonium production as a function of the nuclear mass. The $`c\overline{c}`$ pair travels along different lenghts, spending different time intervals in the nuclei. Therefore, experiments with different $`A`$ explore the time dependence of the $`c_{n,\mathrm{}}`$ amplitudes (Fig.1). We present results on $`\text{}J/\psi \text{}`$ and $`\psi ^{}`$ (Fig.2) in comparison with the experimental values by the NA50 Collaboration , obtained for $`pA`$ collisions with an impinging proton energy of $`450`$GeV. The data are given as branching ratio times cross section divided by the nuclear mass, which differs from the suppression by a constant. (For the comparison, we normalize the data by this constant determined from the average of the ratio experimental/calculated values.) We can see that our results are in good agreement with the experiment both for $`\text{}J/\psi \text{}`$ and $`\psi ^{}`$. In Fig.3 we give predictions for the suppression of $`\chi _c`$ as a function of the nuclear mass, for different values of $`\gamma `$. The largest suppression is observed when $`\gamma `$ is the lowest, corresponding to a longer time spent in the nucleus.

Let us now look at the $`x_F`$ dependence. The suppression factor is often parametrized as $`S_A^\psi (x_F,\sqrt{s_{pp}})=A^{\alpha _\psi (x_F,\sqrt{s_{pp}},A)1}`$ . In fig.4 the $`\alpha `$ values for $`\text{}J/\psi \text{}`$ and $`\psi ^{}`$ are compared to experimental data on $`pA`$ collisions, obtained by the E866/NuSea Collaboration . The impinging proton energy is 800 GeV and the label A is the ratio of the atomic numbers of tungsten over beryllium. In this region of small positive and negative $`x_F`$ values, our calculations are in good agreement with the experimental values for $`\text{}J/\psi \text{}`$, whereas for $`\psi ^{}`$ they overpredict the data. In fact, in the region of small positive $`x_F`$, the $`c\overline{c}`$ can traverse the nucleus before developing into a $`\psi `$ meson. In this kinematical regime, other effects missing in our model, like for example gluon shadowing , can play a role.

Finally, we present predictions for the $`\alpha `$ values (Fig.5) and for the suppression factors (Fig.6) for $`\text{}J/\psi \text{},\psi ^{}`$ and $`\chi _c`$ in the negative $`x_F`$ region, where experimental data from the HERA-B Collaboration  at DESY and the PHENIX Collaboration at RHIC  will soon be available. As we can see, both $`\text{}J/\psi \text{}`$ and $`\chi _c`$ present a rather flat behavior, whereas $`\psi ^{}`$ is very sensitive to the $`x_F`$ values. Besides, the results on $`\psi ^{}`$ also show a stronger variation with the mass of the nucleus. This can be understood from the fact that the effect of the coupling between the charmonium eigenstates is much stronger for $`\psi ^{}`$ than for $`J/\psi `$ (and $`\chi _c`$), as one can see from Fig.1.

Fig.5 also shows that the $`\alpha `$ parametrization of the suppression (often used with $`\alpha `$ independent of $`A`$) is especially not good in the negative $`x_F`$ region, since there is a strong dependence of the $`\alpha `$ values on the nuclei.

In summary, we have studied the production of charmonium ($`J/\psi ,\psi ^{},\chi _c`$) in relativistic proton-nucleus collisions. We focussed on the nuclear effects on the $`c\overline{c}`$ time evolution and on the charmonium production when the time for a $`c\overline{c}`$ to develop to a $`\psi `$ meson is short compared to the time necessary to traverse the nucleus. This kinematical region may be explored experimentally by looking at slow enough $`c\overline{c}`$ and therefore charmonium produced in the negative $`x_F`$ region. The quantum-mechanical model used here includes a realistic potential binding the $`c\overline{c}`$ pair and an imaginary potential, describing the collision of the $`c\overline{c}`$ with the nucleons and depending on the dipole charmonium-nucleon cross section. The imaginary potential introduces transition amplitudes among charmonium eigenstates. We have shown that this produces an interference pattern on the charmonium survival probability, particularly important for $`\psi ^{}`$, and affects the charmonium suppression factors. The results are compared to experimental data from the NA50 and E866/NuSea Collaborations. We have made predictions on the suppression factors of $`J/\psi `$, $`\psi ^{}`$ and $`\chi _c`$, for different nuclei in the fragmentation region of the nucleus, where experiments will be performed soon by the HERA-B Collaboration at DESY and the PHENIX Collaboration at RHIC. We expect that the exploration of the negative $`x_F`$ region will bring new information on the charmonium-nucleon cross section and on the initial conditions which are not well known, especially for $`\chi _c`$. A better understanding of the mechanisms for charmonium production in proton-nucleus collisions in this yet unexplored region will help us in the interpretation of the measurements on relativistic nucleus-nucleus collisions.

We are very grateful to Jörg Hüfner for having interested us in this problem. We thank him for his useful suggestions and comments all along the realization of this work, and for his very careful reading of the manuscript. We thank Boris Kopeliovich and Stephane Peigné for very useful discussions and comments, as well as Yuri Yvanov for his kind help. One of us (C.V.) acknowledge all the people at the Institut für Theoretische Physik of the University of Heidelberg, for the very warm hospitality during her stay.


format: text

text:
# Comment on “Ground State Phase Diagram of a Half-Filled One-Dimensional Extended Hubbard Model”

(August 20, 2003)

In jeckelmann , Jeckelmann argued that the recently discovered bond-order-wave (BOW) phase nakamura ; pinaki ; tf  of the 1D extended Hubbard model does not have a finite extent in the $`(U,V)`$ plane, but exists only on a segment of a first-order spin-density-wave–charge-density-wave (SDW-CDW) phase boundary. We here present quantum Monte Carlo results of higher precision and for larger system sizes than previously pinaki . Using a direct finite-size scaling of the BOW correlations, we reconfirm that the BOW phase does exist a finite distance away from the phase boundary, which hence is a BOW-CDW transition curve. We only address the existence of the BOW phase and focus on a single value, $`U=4`$, of the on-site interaction.

We first determine the critical value $`V_c`$ of the nearest-neighbor interaction. Fig. 1(a) shows the $`V`$-dependence of the charge susceptibility $`\chi _c(q)`$ at the smallest non-zero wave-number, $`q=2\pi /N`$, for different system sizes $`N`$ (for the definition of $`\chi _c`$, see pinaki ). The narrowing of the peak with increasing $`N`$ and the convergence of the height indicate a charge gap for all $`V`$ except at $`V_c`$, i.e., in disagreement with jeckelmann  we find a continuous phase transition at $`U=4`$. The peak position gives $`V_c=2.1602\pm 0.0003`$, which improves considerably on the estimate $`V_c=2.150\pm 0.010`$ reported in jeckelmann .

We next choose $`V=2.10`$, where according to jeckelmann  there should be no long-range BOW. In Fig. 1(b) we show the corresponding correlation function $`C_B(r)`$ at the longest distance, $`r=N/2`$, in periodic chains with up to $`N=512`$ sites. As a function of $`1/N`$ for large $`N`$, the data scale linearly to a value which corresponds to dimerization $`\delta =\sqrt{C_B}/20.053`$ (as defined in jeckelmann ). It is not clear why the density matrix renormalization group (DMRG) calculations of jeckelmann  failed to detect this rather strong BOW order. Dimerization was observed only on the transition curve to the CDW phase jeckelmann , where we instead find coexisting critical BOW and CDW fluctuations. The origin of this discrepancy at $`V_c`$ could be that $`V_c`$ was not determined to sufficient accuracy in jeckelmann .

In summary, our improved calculations do not agree with the phase diagram presented in jeckelmann  but reconfirm the finite extent of the BOW phase pinaki  as first suggested in nakamura . The advantage of the Monte Carlo method we have used pinaki ; sse  is that results can be obtained for large periodic systems, which are better suited for finite-size scaling than the open boundary conditions typically used in DMRG calculations on long chains. With the recently improved simulation algorithm sse  that we have used here, we hope to be able to determine the full extent of the BOW phase in the $`(U,V)`$ plane.

This work is supported by the Academy of Finland, project 26175 (AWS), and by the NSF under Grants No. DMR-99-86948 (PS) and DMR-97-12765 (DKC). The calculations were carried out at the NCSA in Urbana, Illinois, and at the CSC in Helsinki.


text:
# 2-D spectroscopy and modeling of the biconical ionized gas in NGC 4388Based on observations at SAO 6-m telescope (Russia)

## 1 Introduction

The prevailing picture of an AGN structure, the so-called Unified Model, provides for a central supermassive Black-Hole (BH), whose gravitational potential energy is the ultimate source of the AGN luminosity. Matter falling into the BH loses angular momentum through viscous and turbulent processes in an accretion disk, resulting into emission of photons in the ultraviolet, soft-X and hard-X ray wavelength domains (Urry & Padovani (1995)).
According to the Unified Model the central source is surrounded by a dense obscuring thick torus which is able to restrict the emergent radiation from the nucleus to a bipolar cone with an opening angle determined by the torus geometry (Pogge (1989)). Imaging in optical emission lines shows in some cases high ionization regions with a conical and/or biconical morphology. These so-called ionization cones are direct evidence that radiation escapes anisotropically (Mulchaey et al. (1994)).
In this frame we have carried out an observational campaign to study the circumnuclear and extranuclear environment of a spectroscopically selected sample of Seyfert galaxies (Rafanelli et al. (1995)) by means of the integral field spectroscopy. This technique offers much more advantages than the traditional long-slit spectroscopy, because it allows us to observe extended regions of nearby galaxies providing simultaneously recorded spectral information ($`\lambda `$) along two spatial directions ($`\alpha `$ and $`\delta `$).
Here we present results from the observations of the Seyfert-2 galaxy NGC 4388, a nearby (z=0.0084) highly inclined spiral galaxy ($`i`$ 75°), believed located in the core of the Virgo cluster. Indeed observed distortions in the outermost optical isophotes are a clue of the possible interaction between this galaxy and the other members of the cluster.
The IUE spectra of NGC 4388 revealed the presence of an UV source occulted from direct view and producing an ionizing continuum with a steep spectral slope (Ferland & Osterbrock (1986); Kinney et al. (1991); Kinney et al. (1993)).
Narrow-band H$`\alpha `$ +\[N II\] and \[O III\] $`\lambda `$5007 images pointed out the presence of emission line regions with highly ionized gas distributed into two opposite cones: a S-W cone extended to 30″ and a N-E plume extended up to 50″ above the galactic disk. Spectroscopic analysis confirmed the photoionization by the central power-law continuum as the principal excitation mechanism (Pogge (1988); Petitjean & Durret (1993)).

Kinematic data obtained with Fabry-Perot Interferometer (Veilleux et al. 1999b ) revealed non-rotational blueshifted velocities in the extraplanar gas northeast of the nucleus, likely produced by a bipolar outflow.
Using deep narrow-band imaging obtained with the Suprime-Cam mounted at the Subaru Telescope, Yoshida et al. (2002; hereafter YOS02) have recently found a very large emission line region (VEELR) extended up to $``$ 35 kpc northeast of the nucleus. It consists of many gas clouds or filaments, part of which (within 12 kpc) are clearly ionized by the nuclear radiation. They claim that the most plausible origin of this ionized gas is tidal debris due a past interaction with a gas-rich dwarf galaxy.
Moreover at radio wavelengths NGC 4388 shows a two sided, asymmetric structure with a flat spectrum central component, offset by $`\mathrm{\Delta }\alpha `$ =1.2″ and $`\mathrm{\Delta }\delta `$= 2.6″ from the apparent optical nucleus, a diffuse blob to the North and an elongated feature to the South with a compact blob in it (Hummel & Saikia (1991)). Several radio maps have shown that this feature is clearly associated with the S-W ionization cone, suggesting the possibility of an interaction between the radio jet and clouds of the ionized gas (Falcke et al. (1998)). Such interaction generates shock fronts, which cannot be neglected in the interpretation of the spectra.

In this paper we investigate the physical conditions in the extended biconical ionized gas of NGC 4388 through the modeling of the optical emission line and continuum spectra observed in different regions. Both the ionizing radiation from the active nucleus and the shock effects are consistently accounted for. In Sect. 2 the observations and modeling of the S-W cone are presented. In Sect. 3 the N-E cone is considered and the VEELR is modeled on the basis of YOS02 data. The continuum SED is discussed in Sect. 4. Concluding remarks follow in Sect. 5.

## 2 The South-West cone

### 2.1 Observations and data reduction

NGC 4388 was observed in March 1998 at the 6-m telescope of the Special Astrophysical Observatory (Russia) with the MultiPupil Fiber Spectrograph (MPFS), an integral field unit made by an array of $`16\times 15`$ microlenses coupled with a close-packed bundle of optical fibers, which carry the incoming signal to the spectrograph. Each of the 240 elements was covering a region of 1″$`\times `$1″ corresponding to a total field of view of about 2.6$`\times `$2.4 kpc (assuming H0 = 75 km s-1 Mpc-1). Two 600 sec exposures were taken using a TK-1024 CCD, which has a pixel size of 24 $`\mu `$m, and orienting the largest side of the array at P.A. $``$34° (Fig. 1). The grating was chosen to give a spectral range of $``$ 4400–7100 Å, with a resolution of 10 Å and a dispersion of 5.3 Å/px.

The reduction of the data was carried out by means of a software package, developed within our group, and running under IRAF111 IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation environment.
After the usual steps of bias and dark subtraction, all the frames were cross-correlated, in order to compensate for the shifts between them caused by mechanical flexions of the instrument. Then a flat field was chosen as reference image to obtain the geometrical scheme of the spectra arrangement: the position, width and layout along the dispersion direction of every spectrum stored in it were defined and the resulting scheme was subsequently applied to all the images in order to extract for each of them the 240 1-D spectra. The dispersion solution for wavelength calibration was determined by fitting the line positions of He-Ne-Ar lamp exposures with low order Chebyshev polynomial functions. The sky subtraction was performed by averaging the output of 8 fibers devoted to the night sky observation and placed at a distance of 4.5′ from the center of the field of view. Then the data were flux calibrated by observing the spectrophotometric standard star HZ 44: all its flat-fielded, wavelength calibrated and sky subtracted spectra were summed together in order to collect the total incoming flux, to obtain the function necessary to convert counts in physical units and to account for the sensitivity response of the system at different wavelengths. Finally, cosmic rays were removed by combining corresponding spectra of the two exposures of NGC 4388 and a correction for Galactic extinction ($`\mathrm{A}_\mathrm{V}=0.08`$) was applied. Fig. 3 shows some examples of the reduced spectra chosen with different signal-to-noise ratio in the field of view.
The positions, widths and fluxes of all the emission lines detectable 3$`\sigma `$ over the continuum were measured using the IRAF task SPLOT. In case of close blending like H$`\alpha `$ +\[N II\] $`\lambda \lambda `$6548,6583 or \[S II\] $`\lambda \lambda `$6717,6731, better results in performing an interactive multigaussian fit of the total profile could be obtained by using the MIDAS package ALICE, whose input parameters (a first estimate of amplitude, center and sigma of the fitting gaussian) are independently fixable.

### 2.2 Data analysis

In order to get accurate estimate of line fluxes and ratios it is necessary to take into account that emission lines in galaxies are generally affected, and in some cases even dominated, by the underlying starlight contribution. It is possibile to remove the stellar features by subtracting an absorption-line ”template” galaxy spectrum built by means of spectral synthesis techniques (see e.g. Bonatto et al. (1989)), or directly obtained from the observation of early-type spirals. In our case the second approach was used, chosing a spectrum of a S0 galaxy extracted from our database which best matched the stellar features of most of NGC 4388 spectra, and following the suggestions given by Ho et al. (1993). As a whole, the resulting correction was considered satisfactory (an example is given in Fig. 3). The metallic lines (Mg I $`\lambda `$5175, Fe I $`\lambda `$5269 and Na I $`\lambda `$5892) were almost always canceled out after the subtraction. Only in some cases the stellar features of the template have been strengthened to better match the H$`\beta `$ absorption, producing an overcorrection of the metallic absorptions. This can happen since stellar populations cannot be omogeneus everywhere in the field of view. The resulting variations of line intensities have been estimated weak ($`<10`$% for H$`\beta `$ ) in spectra with very bright emission lines, significant ($`1030`$%) in most of the spectra and very important in few spectra (up to 60%).

The emission line fluxes, measured after the starlight contamination removal, were corrected for internal reddening by means of the theoretical Balmer decrement value H$`\alpha `$ /H$`\beta `$ = 2.86 for Case B recombination at electronic temperature $`\mathrm{T}_\mathrm{e}`$=104 K and assuming the reddening law by Miller & Mathews (1972).

Several images of the galaxy within the field of view were reconstructed at different wavelengths, as for instance the map of the visible continuum at 5500 Å. The position of its maximum intensity is located to the N-W of the observed peak of the emission line maps, which we have decided to assume as position of the photoionization source (hereafter nucleus). We do not confirm the large offset between the location of the radio source with flat spectrum and the visible continuum peak indicated in Hummel & Saikia (1991) and references therein: the radio source, the visible continuum and the emission lines have their maxima intensities close to each other (separation $`1`$″). We have discarded the idea to adopt the position of the radio source, which should indicate the real location of the AGN, as nucleus in our field of view since the spatial resolution of our data is only 1 ″/px and an accurate astrometric calibration of the integral field images is very hard to obtain.

The map of the \[O III\] $`\lambda `$5007/H$`\beta `$ emission line ratios revealed part of the S-W ionization cone (Fig. 4). The cone has the axis at P.A. $``$200° and an aperture angle of $``$70° , less than the value measured by Pogge (1988), probably due to our smaller field of view.

The gaseous regions within the ionization cone have been investigated. Their intrinsic FWHM, estimated from the gaussian analysis of the emission lines after having removed the instrumental component ($``$10 Å), range in the interval $``$200–700 $`\mathrm{km}\mathrm{s}^1`$, with a median value of $``$480 $`\mathrm{km}\mathrm{s}^1`$, velocities usually observed in the Narrow Line Regions of Seyfert galaxies. The diagnostic emission line ratios \[N II\] $`\lambda `$6583/H$`\alpha `$ , \[O I\] $`\lambda `$6300/H$`\alpha `$ , \[S II\] $`\lambda \lambda `$6716+6731/H$`\alpha `$ and \[O III\] $`\lambda `$5007/H$`\beta `$ (Veilleux & Osterbrock (1987)), show the typical properties of gas ionized by a non-thermal source (Fig. 5).

Moreover the spatial distribution of the ionization degree within the cone, evaluated by means of the \[O III\] $`\lambda `$5007/H$`\beta `$ ratio as a function of the distance from the nucleus, indicates that the gas is highly ionized even at large distances from the AGN.
In order to understand if the AGN could account for such observed ionization, the H$`\alpha `$ luminosity of each spectrum was calculated from the reddening corrected H$`\alpha `$ fluxes. Then the number of ionizing photons needed to produce such luminosities, $`\mathrm{Q}_{\mathrm{ion}}=7.3\times 10^{11}\mathrm{L}(\mathrm{H}\alpha )`$ photons sec-1 (Kennicutt (1998)), was evaluated and later compared with $`\mathrm{Q}_{\mathrm{nuc}}^{}=\mathrm{Q}_{\mathrm{nuc}}\mathrm{\Omega }/4\pi `$, where $`\mathrm{Q}_{\mathrm{nuc}}`$ is the number of photons coming from the nucleus, diluted by the covering factor $`\mathrm{\Omega }/4\pi `$, and $`\mathrm{\Omega }`$ is the solid angle of a region seen from the nucleus.
By assuming that the nuclear source is an isotropic emitter of radiation occulted along our line of sight, we attempted to estimate $`\mathrm{Q}_{\mathrm{nuc}}`$ by taking the average of the calculated number of photons ionizing the regions surrounding the nucleus, after having removed the effect of the covering factor ($`\mathrm{Q}_{\mathrm{ion}}4\pi /\mathrm{\Omega }`$). A value of $`\mathrm{Q}_{\mathrm{nuc}}4.87\times 10^{52}`$ photons sec-1 was obtained, in good agreement with the 4.98$`\times `$1052 photons sec-1 computed by Colina (1992) using the IUE spectrum of the galaxy and the flux at 912 Å.
The resulting $`\mathrm{Q}_{\mathrm{ion}}/\mathrm{Q}_{\mathrm{nuc}}^{}`$ ratio concentrates around a value of 2.5. This indicates that the active nucleus is the dominant ionizing source for the regions within our field of view, in agreement with what observed in the VO diagrams. In spite of that, the ionization parameter $`\mathrm{U}=\mathrm{Q}_{\mathrm{nuc}}/4\pi \mathrm{cr}_0^2\mathrm{N}_\mathrm{H}`$, calculated within the cone as a function of the distance r0 from the nucleus, ranges between $`10^{5.0}`$ and $`10^{2.8}`$, a result clearly in contrast with the already pointed out high ionization level of almost all the regions. A total hydrogen density $`\mathrm{N}_\mathrm{H}=3\times 10^2\mathrm{cm}^3`$ is assumed, on the basis of the the electron density values estimated over the whole field, using the \[S II\] $`\lambda \lambda `$6717/6731 ratio and assuming T=104 K (except where the sulphur lines were not visible or extremely blended). The problem of the missing data has been solved taking the average value of independent linear interpolations of the \[S II\] ratios along the two orthogonal directions of the spectral array. The resulting Ne distribution presented values $`<`$103 cm-3, with a median of $``$300 cm-3. Of course due to the uncertainties of the \[S II\] emission line strengths, the estimation of the electron density is not very accurate and just a simple decrease by a factor of 10 can account for the ionization parameter in the regions far from the nucleus.

### 2.3 The models

The contribution of the shocks to the line intensities is tested by comparing calculated with observed spectra. The code SUMA (Contini & Viegas 2001, and references therein ), which accounts consistently for the photoionization flux from an external source and for the shocks is adopted.

The general model assumes matter bound gaseous and dusty clouds, which move outwards from the source and emit the line and continuum spectra. A shock front forms on the outer edge of the clouds, while the inner edge is reached by the radiation flux from the source.

The input parameters are : the radiation flux intensity, $`\mathrm{F}_\mathrm{H}`$ (in photons $`\mathrm{cm}^2\mathrm{s}^1\mathrm{eV}^1`$ at 1 Ryd), the spectral index, $`\alpha `$ (= 1.5 for all models), the shock velocity, $`\mathrm{V}_\mathrm{s}`$, the preshock density, $`\mathrm{n}_0`$, the preshock magnetic field, $`\mathrm{B}_0`$ (= $`10^4`$ gauss), the dust-to-gas ratio, d/g, and the geometrical thickness, D, of the emitting clouds. Cosmic relative abundances (Allen (1973)) are adopted.

Composite models lead to a complex distribution of the temperature and the density throughout the clouds. In radiation dominated (RD) models the effect of the flux from the active center (AC) prevails on the shock, while in shock dominated (SD) models the flux is absent ($`\mathrm{F}_\mathrm{H}`$= 0). When shocks and photoionizing fluxes act on the opposite sides of the cloud, two regions of relative high temperature appear near the edges (see Sect. 2.4). Higher temperatures correspond to higher shock velocities (T $``$ $`\mathrm{V}_\mathrm{s}`$2), while the radiation flux cannot heat the gas to $`>`$ 2 - 3 104 K. This temperature leads to high O+2/O oxygen fractional abundance, so the \[OIII\] line strength depends on $`\mathrm{F}_\mathrm{H}`$.

The cooling rate ($``$ n2, where n is the density of the gas) increases with compression downstream of the shock leading to cool gas in the internal region of the cloud. The larger the cool region the strongest the low ionization level lines. The relative sizes of the hot and cool gas regions depend, therefore, on the shock velocity, on the intensity of the flux from the AC, on the preshock density, and on the geometrical thickness of the clouds.

A large grid of models is calculated in order to select the most suitable parameters to the observational evidences for each region in the galaxy. Recall that a first choice is roughly obtained considering that the \[O III\] $`\lambda `$5007 line flux depends on the intensity of the nuclear radiation flux, the \[S II\] $`\lambda \lambda `$6717/6731 line ratio depends on the density, as seen before, and the intensity ratios of high to low ionization lines depend on the geometrical thickness of the clouds and on the shock velocity (see Contini & Viegas 2001). The input parameters are then changed consistently in order to obtain the best fit of all the observed data.

### 2.4 \[O III\]/H$`\beta `$ and \[O I\]/H$`\beta `$

The observations cover the SW cone providing line ratios in many different regions. We focus on \[O III\] $`\lambda \lambda `$4959+5007/H$`\beta `$ and \[O I\] $`\lambda \lambda `$6300+6363/H$`\beta `$ line ratios, because \[O III\] lines are generally strong in Seyfert 2 galaxies and are emitted by gas ionized by the radiation flux from the active nucleus, while \[O I\] lines are emitted from gas at lower temperatures and are generally strong in the presence of shocks.
Moreover, Veilleux et al. (1999a) claim that dynamical processes such as entrainment by AGN-powered radio jets generally have a stronger effect on the kinematics of the highly ionized \[O III\] emitting gas than on those of the low-ionization H$`\alpha `$ emitting material.

Therefore, by modeling the observed spectra it is possible to determine the distribution of the radiation flux intensity and of other parameters, e.g. the shock velocity of the gas, the density, and the geometrical thickness of the clouds in the observed region. The reliability of the models is checked comparing the observed spectral energy distribution (SED) of the continuum with model results (see Sect. 4).

The observed \[O III\]/H$`\beta `$ line ratios overlaid to the continuum of the galaxy at 5500 Å are shown in Fig. 6 in a synthetic representation. The \[O III\]/H$`\beta `$ distribution is not smooth. Particularly, high \[O III\]/H$`\beta `$ ratios appear in some regions far from the active center. This suggest that the clouds in those regions have abnormal intrinsic conditions.

In Fig. 7 \[O III\]/H$`\beta `$ versus \[O I\]/H$`\beta `$ from the observations (black triangles) is compared with model calculations. All the models correspond to $`\mathrm{V}_\mathrm{s}`$=100 $`\mathrm{km}\mathrm{s}^1`$. Higher velocities are less fitting, and velocities lower by a factor $`<`$ 2 lead to similar results. The thick solid line refers to shock dominated models with line ratios increasing with decreasing D and $`\mathrm{n}_0`$. The bulk of the data, however, is scattered throughout a region of the diagram better fitted by radiation dominated models. The thin solid line determines roughly the right edge of this region. The numbers in Fig. 7 refer to Table 1, where the models are described. The dashed lines connect models calculated with different D (increasing from left to right) but all corresponding to the same flux intensity $`\mathrm{F}_\mathrm{H}`$, while the dotted lines connect models calculated with different $`\mathrm{F}_\mathrm{H}`$ (increasing from left to right) and corresponding to the same D. Fig. 7 shows that the trend of both types of models is different from the observed one, i.e. \[O III\]/H$`\beta `$ increasing with \[O I\]/H$`\beta `$ .

The increasing trend of the line ratios can be explained either by a slightly larger D and stronger $`\mathrm{F}_\mathrm{H}`$ (models 1 in Table 1) or by very large clouds (D $``$ 1 pc, models 2 in Table 1), corresponding to a lower density, and a flux intensity not exceeding the characteristic values calculated for most of the observed regions inside the cone. The other observed lines at each position are too few and too uncertain to constrain the models, therefore, we suggest that the clouds reached by a strong nuclear flux (model 1) are located in the central region of the cone, while clouds showing the same high \[OIII\]/H$`\beta `$ ratios, but characterized by a low density and a large geometrical thickness (model 2), are located in some external regions of the cone. Extended clouds of rarefied gas appear in the terminal region of jets and winds.

The distribution of $`\mathrm{F}_\mathrm{H}`$ in the different regions of the cone, as it results from modeling is shown in Fig. 8. The range of the corresponding ionization parameters U is 2 10-5 \- 5 10-4.

To better understand the results obtained above, the distributions of the electron temperature (dotted line) and density (solid line) throughout a cloud are given for different models in Figs. 9 and 10. Each figure is divided by a solid vertical line in two halves. The horizontal axis scale is logarithmic and symmetric to have a comparable view of the two sides of the cloud. The shock front is on the left and the radiation flux from the source reaches the right edge.

The geometrical thickness D is a crucial parameter to low-level and, particularly, to neutral lines, because they are emitted from rather cold gas. Therefore, in the bottom pannels the distribution of the fractional abundance of $`\mathrm{O}^{+0}`$/O (dotted line), $`\mathrm{O}^{+1}`$/O (short-dash line), $`\mathrm{O}^{+2}`$/O (long-dash line), and $`\mathrm{H}^{+1}`$/H (solid line) have been plotted. Fig. 9 corresponds to D= 3 $`10^{16}`$ cm (0.01 pc) and Fig. 10 to D= 4 $`10^{18}`$ cm (1.3 pc). It can be noticed that two zones of gas emit the \[O III\] lines, one from the side of the shock, where collisional ionization prevails, and the other from gas ionized by radiation. The ionization rate due to the photoionization flux from the nuclear source, dominates in a large region of the cloud, particularly in models calculated with $`\mathrm{V}_\mathrm{s}`$=100 $`\mathrm{km}\mathrm{s}^1`$, where the gas density and optical thickness are rather low. A large region of cool gas appears inside clouds corresponding to a large geometrical thickness and low density, leading to relatively strong \[O I\] lines.

### 2.5 \[S II\]/H$`\beta `$

Fig. 11 shows \[S II\] $`\lambda \lambda `$6717+6731/H$`\beta `$ versus \[S II\] $`\lambda \lambda `$6717/6731. The low ionization lines are strongly affected by shocks. For models accounting for the shock a large range of electron densities follows from compression in the downstream region. The lines, calculated integrating through the geometrical thickness of the clouds, account for the different conditions of the gas (see section 2.4).

Some significant models, which are the same as those adopted to explain the \[O III\]/H$`\beta `$ ratios in Fig. 7, are indicated by open circles. They show a good agreement to the observed \[S II\] $`\lambda \lambda `$6717/6731 line ratio between 1 and 1.3 in Fig. 11. The low density model gives \[S II\] $`\lambda \lambda `$6717/6731 $``$ 1.4. Notice, however, that the \[S II\] $`\lambda \lambda `$6717+6731/H$`\beta `$ ratios are higher by a factor $`<`$ 2 than the data for most models, indicating that the S/H relative abundance is lower than cosmic. This is a common characteristic to Seyfert 2 galaxies where sulphur is depleted from the gaseous phase. The models are calculated adopting d/g (by number) = 10-14 \- 10-13, indicating that dust is indeed present in the NGC 4388 cone region.

Reminding that compression increases with $`\mathrm{V}_\mathrm{s}`$, models calculated with higher densities and velocities ($`\mathrm{V}_\mathrm{s}`$=300 $`\mathrm{km}\mathrm{s}^1`$) fit the data corresponding to low \[S II\] $`\lambda \lambda `$6717/6731 ($`<`$1), although their observed errors are large. Finally, models with $`\mathrm{V}_\mathrm{s}`$=500 $`\mathrm{km}\mathrm{s}^1`$ hardly fit any data, suggesting that high velocity clouds are relatively few (see Sect. 4).

The results show that the bulk of the shock velocities is lower than the velocities corresponding to the FWHM of the emission lines (200-700 $`\mathrm{km}\mathrm{s}^1`$), indicating that strong shocks do not form in the observed region of NGC 4388, but rather low shock velocities $`\mathrm{V}_\mathrm{s}`$=100 $`\mathrm{km}\mathrm{s}^1`$ and low densities $`\mathrm{n}_0`$$``$100 $`\mathrm{cm}^3`$ dominate.

## 3 The North-East cone

We now consider the N-E cone and the very extended emission-line region (VEELR) discovered by YOS02 in order to have a larger view of the ionization conditions in NGC 4388. The nuclear region of the galaxy has been successfully investigated by Petitjean & Durret (1993), but only a few emission lines were observed at 21″, 25″, and 38″ from the nucleus.

YOS02 observed the \[O III\] $`\lambda `$5007 and H$`\alpha `$ +\[N II\] $`\lambda \lambda `$6548,6583 fluxes from different clouds in the VEELR, which extends northeastward of the galaxy, very far outside the galaxy disk. The ionization cone is traced by the N-E plume and a high ionization cloud located at 10″-15″ northwest of the nucleus. They claim that there is a sudden change in the excitation state of the gas around the line at P.A. 65°, which coincides with the extrapolation of the northern edge line of the S-W cone. Notice that these authors assume as distance of NGC 4388, the distance of the Virgo cluster (16.7 Mpc), which is different from that one used by us for the S-W cone (33.6 Mpc). Therefore, all the distances should be scaled by a factor of 2.

In Figs. 12, 13, and 14 we compare the results of models calculated with the code SUMA with the YOS02 data. The errorbars refer to the typical uncertainty of 20 % estimated by YOS02. In all the diagrams the fluxes of \[O III\] versus H$`\alpha `$ +\[N II\] (in $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$) are presented. Recall that models are calculated at the emitting nebula, while data are observed at Earth. Therefore, we have multiplied the data by d2/r2, where $`r`$ is the distance of each cloud from the center of the galaxy and $`d`$ is the distance of NGC 4388 from Earth.

All models are calculated adopting a geometrical thickness, D = 3 pc, and $`\mathrm{B}_0`$ = 10-5 gauss, which is similar to the magnetic field in the ISM ($`\mathrm{B}_0`$ = 10-6-10-5 gauss).

Shock dominated models, calculated considering only heating and ionization of the clouds by the shock ($`\mathrm{F}_\mathrm{H}`$=0), are presented in Fig. 12. The models are labeled by the shock velocity in $`\mathrm{km}\mathrm{s}^1`$. Large squares correspond to $`\mathrm{n}_0`$=10 $`\mathrm{cm}^3`$, small squares to $`\mathrm{n}_0`$=2 $`\mathrm{cm}^3`$. Relatively low preshock densities lead to the best agreement with the data, however, the SD models give a poor fit to the observed trend. This confirms YOS02 suggestion that gas excitation depends also on the flux from the active nucleus even at such large distances. We have used, therefore, RD models to fit the data. Recall that RD models are also composed because they account for the shock, even if the effect of the photoionizing flux prevails.

We present in Figs. 13 and 14 radiation dominated models calculated with different velocities and fluxes, adopting $`\mathrm{n}_0`$ = 10 $`\mathrm{cm}^3`$ and 2 $`\mathrm{cm}^3`$, respectively. Higher densities do not fit. In Fig. 13 black triangles indicate models calculated with $`\mathrm{V}_\mathrm{s}`$=50 $`\mathrm{km}\mathrm{s}^1`$, white triangles in the left side of the diagram corresponding to $`\mathrm{V}_\mathrm{s}`$=30 $`\mathrm{km}\mathrm{s}^1`$and at the right to $`\mathrm{V}_\mathrm{s}`$=70-90 $`\mathrm{km}\mathrm{s}^1`$. The models are labeled by the values of log($`\mathrm{F}_\mathrm{H}`$). The maximum \[N II\]+H$`\alpha `$ value is explained by $`\mathrm{V}_\mathrm{s}`$=150 $`\mathrm{km}\mathrm{s}^1`$. Dotted lines join the results of models calculated with $`\mathrm{V}_\mathrm{s}`$= 90 $`\mathrm{km}\mathrm{s}^1`$. Notice that, roughly, all the data can be reached by the models represented by the grid.

In Fig. 14 the models are calculated with $`\mathrm{V}_\mathrm{s}`$=50 $`\mathrm{km}\mathrm{s}^1`$ (black triangles), $`\mathrm{V}_\mathrm{s}`$=70 $`\mathrm{km}\mathrm{s}^1`$ (small white triangles), and $`\mathrm{V}_\mathrm{s}`$=100 $`\mathrm{km}\mathrm{s}^1`$ (large white triangles). The adopted fluxes range between log $`\mathrm{F}_\mathrm{H}`$= 7. and 9. It can be noticed that the calculated trends do not correspond to the observed one, however, some data are well reproduced by the models. In particular, the maximum and minimum observed \[N II\]+H$`\alpha `$ are not explained by these models.

We have chosen the best fitting models from Figs. 13 and 14 and we present in Fig. 15 consistent modeling of the entire dataset. The selected models represented by open circles are described in Table 2.

The fluxes of \[O III\] and H$`\alpha `$ +\[N II\] are given in Fig. 16 as function of distance from the center. The data refer to the fluxes observed at the clouds. Fig. 16 shows that the line fluxes decrease with distance, however, different conditions coexist in the VEELR leading to the large scattering.

The emission line flux from the cloud, identified by YOS02 (see their Fig. 3) as C20, located at the extreme northern edge of the cone, is explained by a low flux (log $`\mathrm{F}_\mathrm{H}`$= 6.5) , $`\mathrm{n}_0`$= 2 $`\mathrm{cm}^3`$, and a low shock velocity (35 $`\mathrm{km}\mathrm{s}^1`$). The largest H$`\alpha `$ \+ \[N II\] flux is emitted by the cloud identified as C26, located in the nebula at the southern edge of the eastern cloud. It corresponds to relatively high $`\mathrm{F}_\mathrm{H}`$ (log $`\mathrm{F}_\mathrm{H}`$=9.5) and high shock velocity ($``$ 150 $`\mathrm{km}\mathrm{s}^1`$). Also cloud C25 is slightly dislocated from the correlating curve in Fig. 15. Therefore, it shows conditions different from those of the other clouds. Cloud C1 shows the lowest distance from center and is modeled by log $`\mathrm{F}_\mathrm{H}`$= 9, $`\mathrm{V}_\mathrm{s}`$= 70 $`\mathrm{km}\mathrm{s}^1`$, and $`\mathrm{n}_0`$ = 10 $`\mathrm{cm}^3`$ (Fig. 11). However, C1 also corresponds to a shock dominated model (Fig. 13) with $`\mathrm{V}_\mathrm{s}`$=80 $`\mathrm{km}\mathrm{s}^1`$and $`\mathrm{n}_0`$=10 $`\mathrm{cm}^3`$. We have not enough data to constrain the model and we conclude that many different conditions coexist in the clouds corresponding to the observations; in fact, the flux from the active nucleus, screened by the intervening matter, cannot reach all of the filaments within the observed clouds. The geometrical thickness adopted in the models is D= 3 pc, while the observed clouds show dimensions of the order of 100 pc. Shock dominated models are radiation bound, so the results should not change even adopting lower geometrical thickness.

To fit the data in Figs. 12, 13 and 14, we multiply the models by factors of 10-5, 10-8 and 10-6.5, respectively. They account for the filling factor. Cloud C20 has a surface area of S = 60$`\times `$ 190 pc2 and the distance from center is r = 14.9 kpc (Yoshida et al. (2002), Table 2). From S/(4 $`\pi `$ r2) $`\times `$ ff = 10-6.5 the filling factor results ff = 0.08. For cloud C26, which corresponds to gas with a higher density, we obtain ff = 1.3 10-4. Regarding cloud C1, ff= 1.4 10-4 results for the radiation dominated filaments and ff= 0.14 for the shock dominated filaments.

YOS02 claim that the inner region (r$`<`$ 12 kpc) of the VEELR may be excited by nuclear ionizing radiation, while the excitation mechanism of the outer region (r$`>`$ 12 kpc) is unclear. Our analysis is not complete because the \[O III\] line fluxes were not observed for r$`>`$ 15.6 kpc. It is found, however, that clouds C19-C22 correspond to lower fluxes and lower velocities, because the nuclear flux from the AC is diluted by distance. The ionization parameter U that results from modeling the VEELR ranges between 10-5 and 1.5 10-3. The shock velocities (50-150 $`\mathrm{km}\mathrm{s}^1`$) are in agreement with the radial velocities in the N-E plume, i.e. 100-200 $`\mathrm{km}\mathrm{s}^1`$(Yoshida et al. (2002)).

## 4 The SED of the continuum

The consistency of the models is checked in Fig. 17 where the SED of the bremsstrahlung radiation and re-radiation by dust from different clouds within NGC 4388 is compared with the data. Therefore, in Fig. 17 two lines correspond to each model, one indicating the bremsstrahlung from the gas and the other emission from dust in the infrared. The observational data (from the NED222NASA/IPAC Extragalactic Database) are taken from Boroson et al. (1983), Scoville et al. (1983), de Vaucouleurs & Longo (1988), Soifer et al. (1989), Moshir et al. (1990), Becker et al. (1991), de Vaucouleurs et al. (1991), Gregory & Condon (1991), Kinney et al. (1993), Fabbiano et al. (1992), Zwicky et al. (1961), Douglas et al. (1996), Dressel & Condon (1978), Spinoglio et al. (1995) and Mould et al. (1980). The two data in the soft X-ray range come from Einstein (Fabbiano et al. 1992), integrated on 0.2-4.0 keV, and from ROSAT (Halderson et al. 2001), integrated on 0.1-2.4 keV.

Indeed, the data from the NED were taken with different precisions and in different epochs. Nevertheless, they are very significant because they correspond to the main contributions to the continuum in the different frequency ranges by 1) bremsstrahlung from the photoionized gas ($`\nu `$ $``$ 5 1013\- 5 1015 Hz), 2) bremsstrahlung from the hot gas in the downstream region of shocked clouds ($`\nu `$ $`>`$ 1016 Hz), 3) emission from the old star population ($`\nu `$ $``$ 5 1013-1015 Hz), 4) dust reradiation in the infrared range ($`\nu `$ $``$ 1012-1014 Hz), and 5) radio emission ($`\nu `$ $`<`$ 1011 Hz).

Besides the line and continuum spectra, the SUMA code calculates reradiation from dust in a consistent way. Dust is heated by the nuclear radiation flux; moreover, gas and dust mutually heat each other by collisions in shock turbulent regimes (Viegas & Contini 1994). Therefore, the frequency corresponding to the dust reradiation peak depends on the shock velocity. At high $`\mathrm{V}_\mathrm{s}`$ the grains will reach high temperatures in the post shock region and dust reradiation will peak in the mid-IR, while at low velocities the peak appears in the far-IR. The peak intensity depends on the dust-to-gas ratio.

Sputtering of the grains, which depends on the shock velocity, is calculated throughout the clouds. Small grains are rapidly sputtered, therefore we adopt silicate grains with an initial radius of 0.2 $`\mu `$m.

Emission in the radio range corresponds to bremsstrahlung from cool gas as well as to synchrotron radiation due to Fermi mechanism at the shock front.

The continuum SED corresponding to single clouds are shown in Fig. 17 in order to explain the role of the different models. Then, the contributions from all the clouds are summed up to fit the data consistently in all the frequency domains.

The models obtained for the S-W cone are used to model the SED of NGC 4388, because the emission in the VEELR is very low and can hardly affect the continuum observed from the entire galaxy. The models presented in Fig. 17 are selected from those which better fit the bulk of the data shown in Fig. 7. We will consider three significant models, regarding the shock velocity, the preshock density, the geometrical thickness of the cloud, and the dust-to-gas ratio. The short-dashed lines refer to a model calculated with $`\mathrm{V}_\mathrm{s}`$=100 $`\mathrm{km}\mathrm{s}^1`$, $`\mathrm{n}_0`$=100 $`\mathrm{km}\mathrm{s}^1`$, D=0.17 pc and d/g=1.5 10-14. Notice that d/g = 10-14 by number corresponds to 4.1 10-4 by mass (Draine & Lee 1984) which is the dust-to-gas ratio in the Galaxy adopting a silicate density of 3 g $`\mathrm{cm}^3`$. Dotted lines represent the SED of the model calculated with D=10 pc, $`\mathrm{n}_0`$=30 $`\mathrm{cm}^3`$, log $`\mathrm{F}_\mathrm{H}`$= 9.8, and d/g=$`10^{14}`$. Considering that the FWHM of the line profiles show velocities up to $``$ 500 $`\mathrm{km}\mathrm{s}^1`$, we have added in Fig. 17 also the continuum of a SD model corresponding to $`\mathrm{V}_\mathrm{s}`$= 500 $`\mathrm{km}\mathrm{s}^1`$, $`\mathrm{n}_0`$=300 $`\mathrm{cm}^3`$, and d/g =1.5 10-13, in order to fit the data at higher frequencies. This relatively high shock velocity leads to a temperature in the immediate postshock region of $``$ 3.7 106 K. Consequently bremsstrahlung radiation peaks in the soft X-ray domain.

Indeed, the analysis of the ROSAT HRI observations by Matt et al (1994) has shown that the 0.1-2.4 keV X-ray emission is dominated by an extended component ($``$ 4.5 pc). It is interpreted as thermal emission from a collisionally heated plasma and a collection of discrete sources. The corresponding temperature of $``$ 0.4 keV ($`4.6\times 10^6`$ K) is in agreement with the temperature of $``$ 4 106 K calculated by $`\mathrm{V}_\mathrm{s}`$=500 $`\mathrm{km}\mathrm{s}^1`$. No correlation of the soft X-ray emission with the optical ionization cones has been found (Matt et al.). This confirms our previous results (Sects. 2 and 3) showing that rather low ($``$ 100 $`\mathrm{km}\mathrm{s}^1`$) shock velocities dominate in the cones.

Fig. 17 shows that the SED of single models calculated with $`\mathrm{V}_\mathrm{s}`$=100 $`\mathrm{km}\mathrm{s}^1`$ nicely fit the data corresponding to the lowest fluxes in the optical-UV range. The data between 1014 and 1015 Hz are, however, contaminated by radiation from old population stars and are explained by a black body model corresponding to 5 103 K (dash-dotted line). This is common to the continua of Seyfert 2 galaxies (Contini & Viegas (2000)).

Recall that dust and gas mutually heat each other across the shock front, therefore, dust re-radiation from high velocity clouds ($`\mathrm{V}_\mathrm{s}`$= 500 $`\mathrm{km}\mathrm{s}^1`$) peaks in the mid-IR (see Contini et al. (2002)). Our results show that these clouds are rather dusty.

The data in the radio range do not follow exactly the slope corresponding to bremsstrahlung from a cool gas. Yet, they do not fit a power-law distribution (thin solid line), which could be explained by synchrotron radiation created by Fermi mechanism at the shock front. Fig. 17 shows that both mechanisms contribute to the radio emission flux.

Finally, the thick solid lines in Fig. 17 represent the weighted sum of the single-cloud models for gas bremsstrahlung and dust reradiation, separately. The line is interrupted at frequencies lower than 1010 Hz, because self-absorption, which increases with wavelength, strongly reduces bremsstrahlung emission in the radio domain.

The model best fit to the data is obtained adopting different weights for different models. The weights account for the relative number of clouds corresponding to the model.

The ratio of the weights used to fit the data by the three models is 1 : 0.2 : 0.03, the maximum weight corresponding to $`\mathrm{V}_\mathrm{s}`$=100 $`\mathrm{km}\mathrm{s}^1`$and D=0.17 pc and the minimum to $`\mathrm{V}_\mathrm{s}`$=500 $`\mathrm{km}\mathrm{s}^1`$. This confirms that velocities of $``$ 100 $`\mathrm{km}\mathrm{s}^1`$and relatively small D ($`<`$ 1 pc) dominate.

## 5 Concluding remarks

The central regions of the Seyfert-2 galaxy NGC 4388 have been observed by means of the integral field spectroscopy technique. The \[O III\] $`\lambda `$5007/H$`\beta `$ map revealed part of the S-W ionization cone, with P.A. $``$200° and aperture angle of $``$70°.
The analysis of the spectra within the cone shows the non-thermal active nature of the emission line regions, whose ionization degree remain high and almost constant even at large distances from the nucleus. With simple energy balance calculations we estimate the number of the ionizing photons produced by the engine and we demonstrate that even if the photoionization is the main physical process responsible for the excitation of the gas inside the cone, it cannot be considered the only one under way, since it cannot explain the observed high ionization far from the source.

We verify our hypothesis by applying to our spectroscopic data the SUMA code, which accounts for the combined effect of the photoionization from an external source and of the shocks. The results of the models, whose reliability is cross-checked by comparing them with the observed SED of the continuum, show that the ionization cone is characterized by a mixture of low density colliding clouds with large (D $``$ 1 pc) and small (D$`<`$0.1 pc) geometrical thickness. Collisions, in fact, cause turbulence and R-T and K-H instabilities which lead to fragmentation of matter. In particular we are able to explain the high \[O III\]/H$`\beta `$ line ratios observed in the outer edges of the S-W cone by the coupled effect of a low density gas in large clouds (D $``$ 1 pc).

The models show that the bulk of the emitting clouds within the S-W cone has shock velocity of $``$ 100 $`\mathrm{km}\mathrm{s}^1`$, definitively lower than those corresponding to the observed FWHM of the line profiles (200 - 700 $`\mathrm{km}\mathrm{s}^1`$). Adopting the hypothesis that the extraplanar material in NGC 4388 may represent the tidal debris from a recent encounter with a small gas-rich galaxy (Pogge (1988); Veilleux et al. 1999b ), we suggest that *head-on-back* shocks accompany the outwards motion of the clouds, which result when an incoming gas catches up with preexistent gaseous clouds moving in the same direction with a different velocity.

To complete the analysis of the biconical structure of NGC 4388 we have modeled the emission line fluxes observed by YOS02 in the VEELR of the N-E cone. We find that the physical conditions are in agreement with those found in the S-W cone, taking into consideration dilution with distance from the AGN. The results lead to log $`\mathrm{F}_\mathrm{H}`$= 8.3 - 10.5 (U = 2 10-5 \- 5 10-4), $`\mathrm{n}_0`$ = 30 - 100 $`\mathrm{cm}^3`$, and $`\mathrm{V}_\mathrm{s}`$$``$ 100 -300 $`\mathrm{km}\mathrm{s}^1`$ in the S-W cone, and log $`\mathrm{F}_\mathrm{H}`$= 6.5 - 9.5 (U = 10-5 \- 1.5 10-3), $`\mathrm{n}_0`$ = 2 - 10 $`\mathrm{cm}^3`$, and $`\mathrm{V}_\mathrm{s}`$= 30 - 150 $`\mathrm{km}\mathrm{s}^1`$ in the VEELR of the N-E cone. We prefer to compare the radiation fluxes instead of the ionization parameter, because U depends on the density of the gas in the edge of the cloud facing the radiation source. So, considering that the density decreases outward the galaxy, dilution is less prominent. The values of U in the nuclear region, calculated by Petitjean & Durret (1993), are much higher, U = 5 10-3 \- 0.014, indicating that different conditions characterize the gas in this region.

Concluding, densities, shock velocities, and the magnetic field are lower in the VEELR, suggesting that beyond the edges of the VEELR the emitting clouds are merging with the ISM.

###### Acknowledgements.

We are very grateful to the referee, D. Alloin, for precious comments which highly improved the presentation of the paper. M. Contini is grateful to the Astrophysikalisches Institut Potsdam and the Department of Astronomy Padova for warm hospitality. We are grateful to Prof. V. Afanasiev and Dr. S. Dodonov of the Special Astrophysical Observatory (Russia) for their decisive contribution in observations and useful discussions about data reduction and analysis. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.


text:
# Effects of resonant single-particle states on pairing correlations

## Abstract

Effects of resonant single-particle (s.p.) states on the pairing correlations are investigated by an exact treatment of the pairing Hamiltonian on the Gamow shell model basis. We introduce the s.p. states with complex energies into the Richardson equations. The solution shows the property that the resonant s.p. states with large widths are less occupied. The importance of many-body correlations between bound and resonant prticle pairs is shown.

One of current topics is concerned with the microscopic structure of nuclei far from the $`\beta `$-stability line. In such nuclei, weak binding of nucleons is expected to cause exotic situations where resonant and continuum single-particle (s.p.) states contribute considerably to nuclear properties even at low energy. The problem has been considered since an early time Migdal , and recent experimental study of drip-line nuclei has stimulated a new interest in it. The resonant and continuum s.p. states could enhance pairing correlations Bennett ; Dobac  which are most important in nuclei once the Hartree-Fock (HF) or shell-model mean field is formed. They could affect binding energies of bound states and nuclear resonant states, and could play important roles in properties of these low-energy states. Recently, theoretical studies Krup ; Sand ; Betan ; Michel  on the contributions from the resonant and continuum states have been reported. Refs. Krup ; Sand  incorporated the effect of the resonant continuum or resonant states in the HF+BCS approximation. Refs. Betan ; Michel  studied the effects of the resonant and continuum s.p. states on one pair state outside a core.

The purpose of this paper is to investigate the effects of the resonant s.p. states on many-body correlations in many particle systems (including correlations between bound and resonant particle pairs), on the shell-model basis. We treat the problem by exactly solving the pairing Hamiltonian including resonant s.p. states. Namely, we extend the Richardson equations Rich1  so as to include s.p. states of complex energies and solve them. The Richardson equations have a defect that the pairing interaction strength cannot be changed for bound and resonant s.p. states. Still, this approach is very useful for getting a general perspective of our subject.

We consider a system of Fermions in which the mean field and residual interaction are approximately represented by a one-body potential and the pairing force. We denote a Fermion by $`c_{\nu \sigma }^{}`$ ($`\sigma =\pm `$ denote time-reversal conjugate states) and a pair operator by $`S_\nu ^{}=c_{\nu +}^{}c_\nu ^{}`$. The eigenstates of our system are classified by the number of unpaired Fermions. The states with no unpaired Fermion are lower in energy than the others. We pay our attention only to such pair-correlated states composed of the basis states $`S_\mu ^{}S_\nu ^{}\mathrm{}|0`$. In this subspace, the pairing Hamiltonian is written as $`H=_\nu 2\epsilon _\nu N_\nu G_{\mu \nu }S_\mu ^{}S_\nu `$, where $`\epsilon _\nu `$ denotes the s.p. energies and $`G`$ is the pairing force strength. The operator $`N_\nu =S_\nu ^{}S_\nu `$ satisfies the relation $`[N_\nu ,S_\nu ^{}]=S_\nu ^{}`$ from the Pauli principle $`c_{\nu \sigma }^{}c_{\nu \sigma }^{}=0`$, and counts the pair number in the level $`\nu `$. Richardson Rich1  proved that a pair-correlated eigenstate is expressed as $`_{k=1}^N\{_\nu S_\nu ^{}/(2\epsilon _\nu z_k)\}|0`$, using the roots $`z_k`$ of the coupled equations

$$\underset{\nu }{}\frac{1}{2\epsilon _\nu z_k}=\frac{1}{G}+\underset{k^{}(k)}{}\frac{2}{z_k^{}z_k}.$$
(1)
Here, the number of the equations (1) is the number of pairs $`N`$. For an eigenstate, the total energy and occupation probability of the s.p. level $`\nu `$ are given by

$`E=H={\displaystyle \underset{k=1}{\overset{N}{}}}z_k,`$ (2)
$`v_\nu ^2=N_\nu =G^2{\displaystyle \frac{d}{dG}}{\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle \frac{1}{2\epsilon _\nu z_k}}.`$ (3)
Although Richardson Rich1  supposed $`\epsilon _\nu `$ to be real number, his proof is not relevant to whether $`\epsilon _\nu `$ are real or complex. We can therefore introduce resonant s.p. states with complex energies to our pairing problem by considering a Woods-Saxon potential. In compensation for this, however, the Hamiltonian stops being Hermitian. Even when the s.p. energies $`\epsilon _\nu `$ are real, some of the roots $`z_k`$ are complex conjugate (the total energy of the system is real). If any of $`\epsilon _\nu `$ are complex in Eqs. (1), all the roots $`z_k`$ are complex but not complex conjugate. If the resonant s.p. states participate in the pairing correlations, therefore, the total energies of the pair-correlated states including bound states become complex in general.

Let us first investigate basic properties of a general system which has resonant s.p. states. We consider a model that has a dozen s.p. levels with equal spacing $`e`$, where the lower eight levels are bound states and the upper four levels have the same resonant width $`\gamma `$,

$$\epsilon _\nu =\{\begin{array}{cc}(\nu 8.5)e\hfill & \text{for }\nu =1,2,\mathrm{}8\hfill \\ (\nu 8.5)ei(\gamma /2)e\hfill & \text{for }\nu =9,10,11,12.\hfill \end{array}$$
(4)
We use the energy unit $`e=1`$ (also for $`G`$ and $`\gamma `$) in the following pages.

When the pair number $`N`$ is small, the Fermi surface lies deeply below the resonant s.p. levels. The resonant states can be neglected as long as the pairing force is not very strong. This is the case of ordinary treatment for the nuclear pairing problem. Calculated results for the $`N=4`$ system are summarized as follows. The system gets additional energy gain when the resonant levels are taken into account and the energy gain increases more rapidly than that in case of no resonant levels as the pairing force becomes stronger. The total energy becomes complex, namely, eigenstates including the ground state have finite resonant width when the resonant s.p. states contribute to the pairing correlations, as mentioned above. Interestingly, the ground state has pair decay width larger than low-energy excited states.

We are interested in the situation of weak binding of nuclei where the Fermi surface lies near the resonant s.p. levels. Fig. 1 shows the effects of the resonant s.p. levels on the $`N=8`$ system with $`\gamma /2=0.5`$, in which the Fermi surface is between the bound and resonant levels. If we do not take account of the resonant levels, there is no configuration mixing like the closed shell ($`v_\nu ^2=1`$ for $`\nu 8`$) and hence the energy gain due to the pairing correlations increases monotonously depending on $`G`$. Fig. 1 indicates more rapid increase of the energy gain when the resonant s.p. states are taken into account. The imaginary part of the total energy shown in the lower graph of Fig. 1 displays interesting behavior. When $`G=0`$, the ground state is the closed-shell state with no decay width and the second pair-correlated state is a resonant state of the Fermion pair $`S_9^{}`$ with the half width $`\mathrm{\Gamma }/2=\gamma =1`$. As the pairing correlations become stronger, the “width” of the ground state increases and after a critical strength of $`G`$ it exceeds the width of the second state which is originated in a resonant state.

How does the imaginary part of the ground-state energy Im($`E`$) depend on the parameter $`\gamma `$? The upper graph of Fig. 2 illustrates Im($`E`$) as a function of $`\gamma /2`$, in case of $`N=8`$ and $`G=0.3`$. The value of Im($`E`$) is proportional to $`\gamma `$ in a small $`\gamma `$ region, but reaches the top at $`\gamma /22`$ after becoming a gentle slope around $`\gamma /21`$. To see the structure change, we calculated the occupation probabilities $`v_\nu ^2`$ of respective s.p. states for the ground state. The real parts of $`v_\nu ^2`$ for the resonant s.p. states are illustrated in the lower graph of Fig. 2. This figure indicates that the occupation probabilities Re($`v_\nu ^2`$) cross the zero point one after another. Especially, Re($`v_9^2`$) of the lowest resonant s.p. level $`\nu =9`$ decreases drastically as $`\gamma `$ increases, and becomes negative at $`\gamma /21.15`$. The occupation probability must be positive to keep the physical meaning. If one of Re($`v_\nu ^2`$) is negative, the state is unphysical. We should, therefore, stop our treatment in the unphysical situation for $`\gamma /21.15`$. This model calculation suggests that the resonant s.p. states with too large widths should be removed from the space of the pairing correlations.

It should be noted that Figs. 1 and 2 give information about how the quantities vary when the pairing interactions are weakened for the resonant s.p. states or widths $`\gamma `$ are narrowed.

When $`N=9`$, the last pair occupies the resonant s.p. levels in our model. This situation is similar to the system of a core plus one pair considered in Refs. Betan ; Michel , which discussed important roles of the resonant and continuum s.p. states in binding of the last one pair. In Fig. 3, we show the total energy of the $`N=9`$ system (ground state) as a function of $`G`$. The real part Re($`E`$) is measured from the base level -64 that is the unperturbed core energy $`\mathrm{\Sigma }_{\nu =1}^82(\nu 8.5)`$. If we regard the $`N=8`$ inner part as a core after Refs. Betan ; Michel , the system is reduced to one pair moving in the 4 resonant s.p. levels. Its energy shown by the dash curve (a) decreases slowly as $`G`$ increases and becomes negative. The negative energy means binding (confinement) of the last pair in the system, according to Ref. Michel . The imaginary part of energy Im($`E`$) remains being 1 in the truncated model (a) which has only 4 resonant s.p. levels with $`\gamma /2=0.5`$.

If we take the 8 bound s.p. levels into account, the ground-state energy is very much lowered as shown by the solid curve (b) in Fig. 3. The energy crosses the zero point at much smaller $`G`$ ($`0.11`$) as compared with (a). Obviously, the inclusion of the bound s.p. states strengthens the correlation energy. This, however, does not mean the binding of the last pair, because the 8 pair system without the last pair also gets large binding energy. The energy difference between the 9 and 8 pair systems, which is equal to the pair decay $`Q`$-value (we denote it by $`Q_{\mathrm{pair}}`$), is plotted by the dot curve (b) in Fig. 3. The curve indicates that the $`N=9`$ system cannot bind the last pair in our model assuming the same potential depth for $`N=8`$ and $`N=9`$. We can say that the correlations of the last pair in the unbound s.p. states with the inner part are very important but whether the last pair is bound or not in the system may depend on complicated conditions. (It is expected that the $`2N+1=17`$ system is more loosely bound due to the blocking effect than the $`2N=18`$ system.) The imaginary part of energy Im($`E`$) increases also in the $`N=9`$ system as $`G`$ increases, in our model without the continuum s.p. states.

The previous work Betan ; Michel  pointed out important roles of the continuum s.p. states in the core plus one pair system. The calculations Betan ; Michel  showed that the continuum states contribute to the binding of the last pair in the system without bound s.p. states like (a). If the continuum s.p. states are correctly treated, therefore, we can expect that the continuum states may contribute to the binding of the last pair in the system including bound s.p. states like (b) and that the imaginary part of energy Im($`E`$) decreases reasonably.

Now let us consider sixteen neutrons of 24O as an example of realistic system. We neglect the continuum s.p. states, because we do not have a reliable treatment of them in our method. We generated the single-neutron basis for 24O by a Woods-Saxon potential with the radius $`R_0=1.2A^{1/3}`$ fm, the surface diffuseness $`a=0.65`$ fm, the potential depth $`V_0=58`$ MeV and the strength of the spin-orbit term $`U_{so}=15`$ MeV Vertse . The s.p. energies obtained are listed in the second column of Table I. The energies $`\epsilon _a`$ are similar to those obtained with the Cwoik code Cwoik , except for the very low energies of the $`2p_{3/2}`$ and $`2p_{1/2}`$ orbits. The resonant widths are large for $`2p_{3/2}`$ and $`2p_{1/2}`$. Using the s.p. energies, we solved the Richardson equations for the lowest three pair-correlated states of 24O, by varying the pairing force strength $`G`$. Though we omit figures, calculated result shows the increase of the pairing correlation energy as $`G`$ increases, which resembles the small $`G`$ region of Fig. 1. The general features shown by the model (4) are confirmed in 24O.

Table I lists the real parts of the occupation probabilities $`v_a^2`$ and pairing gap $`\mathrm{\Delta }=\mathrm{\Sigma }_a\mathrm{\Delta }_a=\mathrm{\Sigma }_aG(j_a+1/2)v_a\sqrt{1v_a^2}`$, for the ground state. Here, we set $`G=0.61`$ MeV so that $`\mathrm{\Delta }`$ is nearly the standard value $`12/\sqrt{A}`$ MeV when $`\gamma _a=0`$. Table I indicates that the widths $`\gamma _a`$ disturb the pairing correlations and pairs hesitate in jumping to the resonant orbits with large widths. The pairing gap $`\mathrm{\Delta }`$ is 2.45 MeV when neglecting $`\gamma _a`$, while it is 2.09 MeV when we take account of $`\gamma _a`$. This feature agrees with that of the HF+BCS calculation Sand . It is important to notice that the resonant orbits with large widths participate little in the pairing correlations. This property is sharply revealed if the widths of $`2p_{3/2}`$ and $`2p_{1/2}`$ are larger. When using another parameters for the Woods-Saxon potential, we had negative values of Re($`v_a^2`$) for $`2p_{3/2}`$ and $`2p_{1/2}`$, which corresponds to the unphysical situation pointed out in Fig. 2.

Kruppa et al. Krup  argued that only resonant s.p. states with narrow widths should be included, to keep physical meaning. If we exclude the orbits $`2p_{3/2}`$ and $`2p_{1/2}`$ following them, the results become very different as shown in the seventh column of Table I. Excluding $`2p_{3/2}`$ and $`2p_{1/2}`$ changes the occupation probabilities Re($`v_a^2`$) very much. The pairing gap $`\mathrm{\Delta }`$ decreases from 2.09 MeV to 1.27 MeV. (It is interesting to note that the values of Re($`v_a^2`$) approach those of the fifth column including $`2p_{3/2}`$ and $`2p_{1/2}`$, if $`G`$ is enlarged so as to give the same gap $`\mathrm{\Delta }2.09`$ MeV.) The significant reduction means that the resonant s.p. states with large widths, if they are included, are little occupied but still contribute to the pairing correlations. If we neglect the last resonant s.p. orbit $`1f_{7/2}`$ further, the pairing correlations are much weakened as shown in the last column of Table I ($`\mathrm{\Delta }=0.73`$ MeV). On the ground-state energy tabulated in the bottom line of Table I, whether including or excluding the resonant s.p. states $`2p_{3/2}`$ and $`2p_{1/2}`$ is more decisive than the effect of the resonant widths $`\gamma _a`$. Exclusion of $`2p_{3/2}`$ and $`2p_{1/2}`$ reduces the energy gain about 1.5 MeV (from 189.29 to 187.77), while the resonant widths $`\gamma _a`$ reduces the energy gain 0.66 MeV (from 189.95 to 189.29). The imaginary part of the ground-state energy Im($`E`$) becomes very small -0.09 MeV from -1.12 MeV, if the resonant orbits $`2p_{3/2}`$ and $`2p_{1/2}`$ with large widths are excluded, which complies with the physical demand Krup .

In conclusion, we have investigated the effects of the resonant s.p. states on the pairing correlations by an exact treatment of the pairing Hamiltonian. The present study shows that many-body correlations between bound and resonant particle pairs are important and that resonant s.p. states with too large widths should not be included in the space of the pairing correlations. We have not completed the study of the contribution from the continuum s.p. states. The interactions related to the resonant s.p. states and especially to the continuum s.p. states must be appropriately evaluated. Our approach using the Richardson equations is not adequate to this problem. It is necessary to use different pairing interactions in the next step. The conclusion that resonant s.p. states with very short life time should be excluded from the space of the pairing correlations must be investigated further in calculations which deal with the continuum s.p. states correctly.


text:
# Quantitative Spectroscopy of Supergiants

line: formation – stars: abundances – stars: evolution – stars: fundamental parameters – supergiants – galaxies: abundances

Dissertation Summary

Electronic mail: norbert@ifa.hawaii.edu

Thesis work conducted at: Ludwig-Maximilians-Universität München

Ph.D. Thesis directed by: R.P. Kudritzki/R. Bender; Ph.D. Degree awarded: June 2002

Received

Blue supergiants of spectral types B and A are the visually brightest stars in spiral and irregular galaxies, with their most luminous members (at $`M_V`$$``$$``$10) outshining entire dwarf galaxies. This characteristic allows us to use them as probes to study the Local Universe in great detail. In principle, already the existing large telescopes and instrumentation facilitate quantitative spectroscopy of these objects as far as the Virgo and Fornax clusters of galaxies. Beyond their challenging stellar atmospheres and opportunities for testing sophisticated non-LTE physics they offer numerous applications to modern astrophysics. Quantitative spectroscopy of supergiants will contribute to improve our understanding of massive star evolution. Galactic abundance gradients and abundance patterns, as will be obtained from studies of large ensembles of supergiants in our own and other galaxies, will foster the understanding of galactochemical evolution. Finally, they are promising independent indicators for calibrating the extragalactic distance scale, by application of the wind momentum–luminosity and the flux-weighted gravity–luminosity relationships (Kudritzki et al. 1999; Kudritzki, Bresolin & Przybilla 2003).

In view of this large potential, the objective of this thesis is to improve the status of quantitative spectroscopy of BA-type supergiants and to provide first applications on a sample of Galactic and extragalactic targets. It is shown that among the model atmospheres available at present, the best suited for analyses of supergiants are line-blanketed classical LTE models. An investigation of the impact of various parameters such as helium abundance and line blanketing on the atmospheric structure shows that for the most luminous objects an accurate treatment of these parameters is essential for a quantitative analysis, whereas the less luminous supergiants react less sensitive. Spectrum synthesis is used to model the line spectra. It is the only technique capable of providing analyses of spectra of different qualities from low to high resolution, and able to cope with heavy line blending at a broad range of signal-to-noise ratios. The stellar parameters are determined from purely spectroscopic indicators, from temperature and gravity sensitive ionization equilibria and the Balmer line wings. Elemental abundances are derived by modelling individual spectral features. Several ten-thousand spectral lines from 28 chemical species are included in the line formation, allowing almost the entire observed spectra to be reproduced. Non-LTE effects become important in blue supergiants, where a strong radiation field at low particle densities favours deviations from LTE. Comprehensive model atoms are therefore constructed for C i/ii, N i/ii, O i and Mg i/ii in order to determine non-LTE level populations (Przybilla, Butler & Kudritzki 2001; Przybilla & Butler 2001; Przybilla et al. 2000, 2001) Highly accurate radiative and collisional atomic data recently determined for astrophysical and fusion research using the $`R`$-matrix method in the close-coupling approximation are incorporated. In addition, model atoms for H, He i, O ii, S ii/iii, Ti ii and Fe ii are adopted from the literature, the atomic data being updated to more modern values in some cases. Thus, an improved treatment for the main elements of astrophysical interest is achieved.

Extensive testing of the atomic data is performed for the nearby bright main sequence standard Vega (A0 V), at well-determined stellar parameters and atmospheric structure. A high-resolution and low-noise spectrum with large wavelength coverage from the visual to the near-IR is used for this purpose. Further tests are performed for the Galactic supergiants $`\eta `$ Leo, HD 111613, HD 92207 and $`\beta `$ Ori, with similar high-quality spectra in order to study the non-LTE effects across the parameter space. Accurate and consistent stellar parameters are derived for these. Non-LTE ionization equilibria of several elements – typically N i/ii, O i/ii, Mg i/ii and S ii/iii – agree simultaneously, provided a realistic treatment of line-blocking is used. These parameters also constitute important input data for further studies of the stellar winds of these objects. Accounting for non-LTE reduces the random errors and removes systematic trends in the analyses. In particular, the improved treatment of electron collisions largely removes long-standing discrepancies in analyses of lines from different spin systems of a given ion. The computed non-LTE line profiles fit the observations well for the different species at a given elemental abundance. In the parameter range covered, all lines from He i, C i/ii, N i/ii, O i/ii and S ii/iii are significantly strengthened by non-LTE effects; Mg ii remains almost unaffected, except for the strongest lines; non-LTE weakening is found for the lines of Mg i, Ti ii and Fe ii in supergiants. The nature of the non-LTE effects is investigated: non-LTE strengthening generally occurs in conjunction with a strong overpopulation of metastable energy levels, while non-LTE weakening is due to overionization of minor ionic species. In extreme cases, as for several strong lines of N i and O i, non-LTE abundance corrections up to a factor of 50 to 100 are found, whereas typical mean non-LTE abundance corrections for the diagnostic lines are within a factor of up to 3. Estimates of the systematic uncertainties in the non-LTE abundance analysis of CNO and Mg are provided. Accounting for these and random errors, it is shown that absolute elemental abundances can be derived with accuracies of $``$0.1 to 0.25 dex (1$`\sigma `$-uncertainties), depending on the element, in contrast to $``$0.2–0.3 dex (only random contribution) typically achieved in previous studies. The statistical significance of the analyses is also largely improved due to the large wavelength coverage of the spectra, with many lines per element being available. In addition, hitherto unaccounted effects on metal line strengths are found from the consistent treatment of microturbulence in the non-LTE computations and line formation.

The abundance analysis for Vega confirmes its status as a mild $`\lambda `$ Bootis star: the light elements, CNO, show an underabundance of $``$0.25 dex when compared to the solar composition, while the heavier elements are depleted by $``$0.55 dex. All four Galactic supergiants have metallicities close to solar. The non-LTE abundances for individual heavier elements in a star cluster around a mean offset to the solar composition, whereas in a classical LTE analysis misleading ‘abundance patterns’ are seen. In particular, LTE analyses tend to overestimate the abundances of the $`\alpha `$–process elements and to underestimate the Iron Group abundances; this systematic effect strengthens with increasing luminosity. It might be suspected that other elements with an atomic structure comparable to the species investigated are susceptible to similar non-LTE mechanisms; these need to be quantified in future studies. The increasing scatter of individual abundances around a mean value with increasing stellar effective temperature and luminosity can be interpreted in terms of neglected effects in the atmospheric modelling. Non-LTE effects on the atmospheric structure become more pronounced at higher temperature, while sphericity and hydrodynamical outflow velocity fields are noticeable only for very luminous objects.

A similar analysis is performed on high-resolution spectra of supergiants in nearby Local Group galaxies, although at lower signal-to-noise ratios. VLT/UVES and Keck/HIRES spectra are available for this (Venn et al. 2000, 2001). One early A-type supergiant near the centre of the dwarf irregular galaxy NGC 6822 is found to be metal-poor by $``$0.55 dex, confirming similar values from the literature, both from stellar and H ii region studies. The non-LTE overionization of Ti ii shows a distinctive strengthening at this low metallicity. Two objects of early-A spectral type are analysed in M 31, both located at a galactocentric distance of $``$12 kpc. The objects have abundances compatible with and slightly above solar, also in good agreement with measurements from close-by H ii regions in M 31. The non-LTE effects follow the same trends as in the Galactic counterparts.

While indications are found that the abundances of the heavier elements in the supergiants cluster around a mean value, distinctive abundance patterns for the light elements are derived, with enriched He and N, depleted C, and O being compatible with the heavier element abundances. The predictions for chemical mixing from recent stellar evolution models, accounting for mass-loss and rotation, can be verified: good agreement is found for the He enrichment and the N/C ratios in the Galactic supergiants, whereas the observed N/O ratios are lower than those from the models. With the presently achieved accuracy in the abundance determination, different stages of stellar evolution – first crossing of the Hertzsprung-Russell diagram from the main sequence to the red vs. a blue-loop scenario with first dredge-up abundances – can be distinguished with confidence. From the Galactic sample, $`\eta `$ Leo is apparently undergoing a blue-loop, while the other stars have directly evolved from the main sequence. In the NGC 6822 and M 31 supergiants similar He enrichments and unchanged O abundances are derived. However, as the diagnostic lines of the main indicators for mixing, C and N, are not covered by the available spectra, further conclusions on their evolutionary status cannot be drawn at present.

Low-resolution spectra of supergiants in the Sculptor Group spiral galaxy NGC 300 and in NGC 3621 in the field, at distances of 2.0 and 6.6 Mpc – well beyond the Local Group –, were recently obtained with FORS1 on the VLT (Bresolin et al. 2001, 2002). In order to perform quantitative analyses of these stars, the applicability of the spectrum synthesis technique is tested on the most luminous Galactic sample star, HD 92207, its spectrum being artificially degraded to FORS1 resolution. It is shown that the stellar parameters can be determined with sufficient accuracy from the spectral classification and the Balmer line strengths to constrain the metallicity to $`\pm `$0.2 dex. Also, individual abundances for a few key elements of astrophysics can be determined at this resolution. Two early A-type supergiants in NGC 3621 and two similar objects in NGC 300 are studied, each among the most luminous stars in these galaxies. The objects have metallicities of slightly sub-solar to $``$0.2 $`\times `$ solar, typically, in concordance with expectation from their position in their host galaxies. The stellar metallicities agree with literature data on abundance gradients for these galaxies, as derived from H ii regions; the reddening of the objects is compatible with being due to the Galactic foreground. Analyses of stars at such distances are performed for the first time.