text:
# Codes for optically thick and hot photoionized media - Radiative transfer and new developments

## 1. Introduction

The average continuum observed in Active Galactic Nuclei (AGN) is shown in Fig. 1 from the optical to X-ray range, as obtained from a composite optical-UV spectrum of Francis et al. (1991) and Zheng et al. (1997), and from an average (slightly modified) soft X-ray spectrum of Laor et al. (1997). It is widely admitted since several years that the optical-UV part of this continuum is produced by a warm (temperature of the order of $`10^510^6`$ K), optically thick or effectively thick medium, most likely an accretion disc. Moreover, one deduces from the spectral distribution in the hard X-ray range, in particular the 30 keV hump observed in many Seyfert 1 galaxies, and from the presence of other features like the iron K line around 7 keV, that this medium is irradiated by a hard X-ray source, which is partly reprocessed and reemitted (improperly named “reflected”). So one sees a combination of the primary and of the reflected spectra (Ross & Fabian 1993, and many subsequent works). In another interpretation the UV-X continuum is produced by clouds from a disrupted disk (Collin-Souffrin et al. 1996, Czerny & Dumont 1998), and in this case one would also see the spectrum emitted by the non illuminated surface of the clouds. In both interpretations the UV-soft X spectrum is emitted by a dense and thick shell (density higher than 1012 cm-3, Thomson thickness higher than unity), irradiated by an X-ray continuum. Similar conditions are met in X-ray binary stars.

It is necessary to compute the structure of such an irradiated medium to predict the observed spectrum. Therefore one has to build “intermediate”codes, between those designed for planetary nebulae or Narrow or Broad Line Regions in AGN, and those designed for stellar atmospheres. These codes have to be valid for thick or semi-infinite dense media, eventually hot, irradiated by a non thermal continuum extending in the hard X-rays, to cover the whole range of situations encompassed in AGN and in binary stars.

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

We have undertaken to build a code which satisfies these requirements, to compute the emission spectrum produced by irradiated thick and hot media like those commonly assumed in AGN close to the black hole, in a wide range of photon energy. Precisely we have built several interconnected codes, which allow more flexibility. The ensemble is far from being perfect and still contains several approximations which restrict its use, but we improve it gradually.

TITAN is designed to study the structure of a warm or hot thick photoionized gas, and to compute its emission - reflection - transmission spectrum from the infrared up to about 20 keV. It solves the energy balance, the ionization and the statistical equilibria, the transfer equations, in a plane-parallel geometry, for the lines and the continuum. Then, given the thermal and ionization stratification, the computation of the emitted spectrum from 1 keV to a few hundreds keV is performed with NOAR which uses a Monte-Carlo method taking into account direct and inverse Compton scattering (it allows also to study various geometries). In a previous paper (Dumont, Abrassart & Collin 2000) TITAN and its interconnection with NOAR were described. We recall here the main characteristics of the code, and we describe recent improvements as well as new results, focussing only on the aspects which are not treated in a standard way.

We briefly summarize below the physical processes (Section 2), the transfer method and the iteration procedure (Section 3). The influence of the physical parameters, of the atomic data, and of the density distribution, are discussed in Section 4. The coupling of TITAN and NOAR is briefly described in Section 5.

## 2. Generalities about TITAN

### 2.1. Physical processes

In TITAN, the physical state of the gas (temperature, ionic abundances and level populations of all ionic species), is computed at each depth, assuming stationary state, i.e.: local balance between ionizations and recombinations, local balance between excitations and deexcitations, local energy balance, and finally total energy balance (equality between inward and outward fluxes).

The ionization equilibrium equations include radiative ionizations by continuum and line photons, collisional ionizations and recombinations, radiative and dielectronic recombinations, charge transfer by H and He atoms, the Auger effects, and ionizations by high energy electrons arising from ionizations by X-ray photons. The emission-absorption mechanisms of the continuum include free-free and free-bound processes, two-photon process, and Thomson electron scattering. Inverse processes (except autoionization for dielectronic recombinations, and recombinations onto excited levels for some ions, see next section) are computed through the equations of detailed balance. Induced processes are taken into account. Energy balance equations include the same processes and Compton heating/cooling.

### 2.2. Atomic data

Hydrogen and hydrogen-like ions are treated as 5 or 6-level atoms. In a first step, in order to save computation time, non hydrogen-like ions were treated with a rough approximation: interlocking between excited levels was neglected and the populations of excited levels were computed separately using a two-level approximation. This approximation does not predict correctly the details of the line spectrum, including the resonance lines themselves. Moreover ionizations from excited levels were not taken into account, though recombinations onto excited states were taken into account in the ionization equilibrium. This approach is is not correct close to LTE, and even far from LTE, if ones wants to predict accurately detailed spectral features (cf. Section 4.2). However it gives a correct overall spectral distribution in this last case. Such are the upper layers of an irradiated atmosphere where at the same time the fractional abundances of heavy elements in high ionization states are important, and the ground level of these ions is strongly overpopulated with respect to LTE. But it is not the case in the deep layers of a disk atmosphere, emitting the UV part of the spectrum, where the density is high, and the influence of the external irradiation is small.

This is why helium-like ions treated as complete 8-level atoms and lithium-like ions treated as complete 5-level atoms have been recently implemented in the code. All processes are thus taken into account for each level. This treatment is valid as far as the upper levels are in LTE with the continuum. To handle intermediate cases, we replace the missing upper levels which are very important especially for helium-like ions by adding an additional recombination rate, not compensated by ionization.

In the future we plane to add several other levels, and to sum the contributions of the higher levels to handle at the same time cases close to and cases far from LTE.

Photoionization cross sections are fitted from Topbase (cf. Cunto et al 1993). Cross sections of neutral and once ionized heavy elements are not well represented, but they do not have any influence on these hot media.

The gas composition include 10 elements (H, He, C, N, O, Ne, Mg, Si, S, Fe), and all their ionic species.

The slab is divided in about 300 layers with variable thicknesses. The code can work with a prescription for the density or for the pressure, such as a constant pressure or a pressure corresponding to hydrostatic equilibrium.

### 2.3. Illumination

The external illumination can be characterized by an “ionization parameter”, which we choose here as $`\xi =4\pi F_{\mathrm{inc}}/n_\mathrm{H}`$, $`F_{\mathrm{inc}}`$ being the frequency integrated flux incident on one side of the slab and $`n_\mathrm{H}`$ the hydrogen number density at this surface. Note that there are several other definitions of the ionization parameter, for instance $`U`$, or $`\mathrm{\Xi }`$, but they correspond basically to the same concept, namely a parameter which determines the ionization state. It is possible to take also into account a flux incident on the back side of the slab.

Presently we mimic a semi-isotropic illumination, but with the ALI transfer method which is implemented (cf. below) it will be possible to take into account the angular distribution of the incident radiation.

## 3. Radiation transfer

We are concerned by media with a continuum optical thickness larger than unity in a range of frequencies rich in emission or in absorption lines, and having a very inhomogeneous structure. Moreover in the hot photoionized layers, the absorption coefficient is weak at all wavelengths and Thomson diffusion is dominant. On the contrary absorption dominates and must be carefully treated to compute the spectrum emitted by the non illuminated cold part of the medium. Radiation transfer requires therefore some attention.

The transfer equation is

$$\mu \frac{dI}{dz}=(\kappa +\sigma )I+\sigma J+\eta ,$$
(1)
where $`z`$ is the distance to the illuminated surface, $`\mu `$ is the cosine of the angle between the normal and the light ray, $`\kappa `$ is the absorption coefficient (for the continuum it is due to photoionizations and free-free transitions), $`\sigma `$ is the diffusion coefficient (here Thomson scattering), $`\eta `$ is the emissivity (due to radiative recombinations and free-free emission), and $`J`$ is the mean intensity, equal to $`I_\nu 𝑑\mu /2`$ (all these quantities are local, $`I`$ depends on the direction and on the frequency $`\nu `$, and J, $`\kappa `$ and $`\eta `$ depend on the frequency).

Introducing the source function $`S`$ and the optical depth $`\tau `$, Eq. 1 can be written:

$$\mu \frac{dI}{d\tau }=I+S$$
(2)
with

$$S=\frac{\sigma J+\eta }{\sigma +\kappa },$$
(3)
which can be separated between J-dependent and J-independent components:

$$S=(1ϵ)J+ϵB.$$
(4)
The formal solution is:

$$J(\tau )=\frac{1}{2}𝑑\mu S(t)e^{(t\tau )/\mu }𝑑t/\mu .$$
(5)
### 3.1. Transfer of the continuum: the ALI method

In the previous version of the code (Dumont et al. 2000) we used the same method for line and continuum transfer, described in Section 3.2. Presently the transfer of the continuum is performed with the Accelerated Lambda Iteration method (ALI), which was set up with the help of subroutines kindly provided to us by F. Paletou (cf. Auer & Paletou 1994, and Kunasz & Auer 1988).

In this method $`J(\tau )`$ is written with the help of an operator $`𝚲`$:

$$J(\tau )=𝚲[S(t)].$$
(6)
ALI uses an operator splitting, described first by Cannon (1973), then by Olson et al. (1986):

$$𝚲=𝚲^{}+(𝚲𝚲^{})\mathrm{and}S=S^{}+\delta S$$
(7)
where $`𝚲^{}`$ is an approximate $`𝚲`$ operator chosen to facilitate the inversion, and $`S^{}`$ has been calculated in the previous iteration. One obtains:

$$\delta S=[1(1ϵ)𝚲^{}]^1[SS^{}].$$
(8)
We use here a parabolic interpolation in the evaluation of the integral giving the intensity from the source function, and an acceleration type Ng (in this method, at every fourth iteration a least squares extrapolation is made in order to reduce the residuals, cf. Ng 1974).

To perform this computation, one needs to know the opacities and emissivities at each frequency as functions of $`z`$, hence the temperature and the populations at each depth, so an iteration procedure is required.

* For each layer, starting from the illuminated side, the ionization and thermal balance equations are solved by iteration; the temperature, the opacity and the emissivity are computed;
* when the back side of the cloud is reached, the transfer is solved with the ALI method, and $`J_\nu (z)`$ is computed;
* the whole calculation is repeated until convergence with the new values of $`J_\nu (z)`$. It is stopped when the energy balance is achieved for the whole slab (i.e. when the flux entering on both sides of the slab is equal to the flux coming out from both sides). This method is almost as slow as the previous one, albeit more secure.

Fig. 2 shows that the difference in the spectrum computed with the previous transfer method and with ALI is very small. It is certainly due to the fact that the continuum optical thickness is never very large (at most 104). This is not the case at the center of intense lines, which are typically 100 times thicker. We are presently implementing ALI for the line transfer and we expect to find more important differences for a few lines like OVIII L$`\alpha `$ which are not completely converged even after 1000 iterations (cf. Dumont et al. 2000).

### 3.2. Line transfer

Since ALI is presently not implemented for the line transfer, it is treated like in the previous version with a simple method based on the Eddington two-stream approximation. In this approximation the transfer equations are written:

$`{\displaystyle \frac{1}{\sqrt{3}}}{\displaystyle \frac{dI_\nu ^+}{dz}}`$ $`=`$ $`(\kappa _\nu +{\displaystyle \frac{\sigma }{2}})I_\nu ^++{\displaystyle \frac{\sigma }{2}}I_\nu ^{}+ϵ_\nu `$ (9)
$`{\displaystyle \frac{1}{\sqrt{3}}}{\displaystyle \frac{dI_\nu ^{}}{dz}}`$ $`=`$ $`(\kappa _\nu +{\displaystyle \frac{\sigma }{2}})I_\nu ^{}+{\displaystyle \frac{\sigma }{2}}I_\nu ^++ϵ_\nu .`$
The mean intensity $`J_\nu `$ is equal to $`(I_\nu ^++I_\nu ^{})/2`$, and the flux $`F_\nu `$ defined as usual by $`I_\nu \mathrm{cos}\theta d\omega `$ is equal to $`(I_\nu ^+I_\nu ^{})2\pi /\sqrt{3}`$. The “reflected” flux is equal to $`I_\nu ^{}(0)2\pi /\sqrt{3}`$, and the outward flux to $`I_\nu ^+(H)2\pi /\sqrt{3}`$. The optical depth and the total optical depth are defined as:

$$\tau _\nu (z)=_0^z\sqrt{3}(\kappa _\nu +\sigma )𝑑z^{}\mathrm{and}T_\nu =\tau _\nu (H)$$
(10)
where $`H`$ is the total thickness of the slab. For $`I_\nu ^+(z)`$ the formal solution of the transfer equations between $`z\delta z`$ and $`z`$ is:

$$I_\nu ^+(z)=I_\nu ^+(z\delta z)e^{\delta \tau _\nu }+e^{\tau _\nu }_{\tau _\nu \delta \tau _\nu }^{\tau _\nu }S_\nu (t)e^{+t}𝑑t.$$
(11)
We approximate the source function $`S_\nu (t)`$ by a constant in the interval $`z,z+dz`$. It leads to:

$$I_\nu ^+(z)=I_\nu ^+(z\delta z)e^{\delta \tau _\nu }+\frac{1e^{\delta \tau _\nu }}{2}[S_\nu (z\delta z)+S_\nu (z)]$$
(12)
and:

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

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

The procedure is then the following:

* for each layer, after computation of the temperature and the source function, $`I_\nu ^+(z)`$ is transferred through the layer according to Eq. 12, while $`S_\nu (z)`$ and $`\tau _\nu (z)`$ are buffered for each layer and each frequency;
* when the back side of the cloud is reached, new values of $`I_\nu ^{}(z)`$ are calculated, starting from the back side where $`I_\nu ^{}(H)`$ is given by Eq. 12, and using Eq. 13 with the previous values of $`S_\nu (z)`$ and $`\tau _\nu (z)`$ (actually we use the values given by several previous iterations to accelerate the procedure with a Ng method).

This method requires a large number of iterations to get a complete convergence of the model, and one cannot be sure to have reached the correct value. It is why we are now changing for the ALI method.

We assume presently complete redistribution in the lines, and partial redistribution in intense lines is mimicked by assuming complete redistribution in a pure Doppler profile. It will be possible to treat partial redistribution with the ALI method.

A line profile is represented by a Voigt function at several frequencies distributed symmetrically around the line center. The frequencies in the profile are chosen to take into account the fact that the Doppler width may decrease by one order of magnitude inside the slab, owing to large temperature variations. Moreover, to solve the statistical equilibrium equations one needs to know the total line intensity weighted by the profile $`J_\nu \varphi (\nu )𝑑\nu `$, while to compute the thermal and ionization equilibria one simply needs the line intensity integrated over frequencies $`J_\nu 𝑑\nu `$. For intense lines, the first integral is dominated by the Doppler core, and the second by the wings. The set of frequencies chosen to describe the line must therefore cover at least two orders of magnitude. The integrals over the line profiles are achieved with a 15-point Gauss-Legendre quadrature. We also assume that the lines do not overlap, which is valid as far as multiplets are treated globally, and because the number of lines taken into account in TITAN is relatively small (presently 580, including subordinate lines in the optical-UV range).

As shown by Fig. 3, the profiles of $`I_\nu ^{}`$ and $`I_\nu ^+`$ for OVIII Ly$`\alpha `$ vary considerably from the surface to the deepest layers.

## 4. Some results of TITAN

We shall call “reference model” a model with the following characteristics:

* it is a plan-parallel slab with a constant density, equal to 1012 cm-3,
* its column density equals 1026 cm-2 (corresponding to a Thomson thickness $`\tau _{\mathrm{es}}`$ of 80),
* it is illuminated on one side by an incident power law continuum proportional to $`\nu ^1`$ from 0.1 eV to 100 keV,
* the ionization parameter is $`\xi =10^3`$ erg cm s-1,
* there is no illumination on the other side,
* the abundances are (in number): H: 1, He: 0.085, C: 3.3 10-4, N: 9.1 10-5, O: 6.6 10-4, Ne: 8.3 10-5, Mg: 2.6 10-5, Si: 3.3 10-5, S: 1.6 10-5, Fe: 3.2 10-5.

In several examples given below, we also use a smaller value of the ionization parameter, $`\xi =10^2`$ (corresponding to a value of the usual ionization parameter $`U`$ equal to 0.88), allowing an easier comparison with computations made with CLOUDY or with XSTAR.

Figs. 4 and 5 show the results corresponding to different values of the ionization parameter, the other parameters being the same as in the reference model. These results were obtained with the previous version of the code. Indeed each model takes a long time to run, and all models have still not been recalculated with the ALI method. However, we have checked in several cases that the results obtained with the new version differ very slightly from the previous ones.

Fig. 4 displays the temperature versus the depth. Close to the surface the temperature is high and reaches the Compton temperature for the highest value of $`\xi `$. The extension of this region (called the “hot skin”) increases with the value of $`\xi `$. Its optical thickness reaches $`\tau _{\mathrm{es}}10`$ for $`\xi =10^4`$. As we shall see later, this is quite different from the case of constant pressure, where the thickness of this zone is always limited to $`\tau _{\mathrm{es}}1`$. This hot skin is almost absent for low values of $`\xi `$. Below the hot skin the temperature decreases smoothly, again at contradistinction with the case of a constant pressure. This decrease is followed by a region of a quasi constant temperature, which is in a state of “quasi-LTE”, i.e. each process important for the energy equilibrium (lines, free-bound from and onto ground levels, free-free) is exactly balanced, due to the high optical thickness and relatively high density. However the region is not in LTE, in the sense that the fractional ionic abundances do not satisfy Saha equations. We shall discuss this issue in more detail in the next section.

On Fig. 5 are displayed the “reflected” spectrum emitted by the illuminated side of the slab, and the “outward” spectrum emitted by the back side. Note that the outward spectrum is distinct from the transmitted one, as it includes also emission (actually our transfer method does not allow to separate emission and transmission, contrary to approximate methods using for instance “on the spot” approximation). The outward spectrum lies only in the optical near UV range, owing to the low temperature of the back side (cf. Fig. 4). This is due to the large column density of the model which does not allow any transmission in the X-ray range. Note the strong edges in emission (including a Balmer edge), and the cut off of the spectrum at a few tens of eV, due to absorption by HeII and by low ionization heavy element species. The reflected spectrum shows a characteristic shape with a smooth “powerlaw” soft X-ray excess, and displays many lines and ionization edges whose intensity decreases with increasing values of $`\xi `$, owing to the increasing influence of pure reflection when the ionization state increases. Note that the Lyman discontinuity is in emission for low $`\xi `$ and in absorption for high $`\xi `$.

Fig. 6 illustrates the influence of the column density on the reflected and on the outward spectra for $`\xi =100`$, the other parameters being the same as in the reference model. First we see that the reflected spectrum “saturates” above 1025 cm-2, because all the available photons have been absorbed and no longer succeed in escaping from the illuminated side. Second, the optical-UV part of the outward spectrum saturates above 1025 cm-2, but not the X-ray part, which is due to transmission. As the column density decreases, the slab becomes more transparent and less emissive, so the outward spectrum becomes basically only the transmitted spectrum, while the reflected spectrum decreases in intensity. Note that the “saturation” of the reflected spectrum would occur at a larger value of the column density for a larger $`\xi `$, as the hot reflecting skin becomes thicker. It is thus necessary to take into account deeper layers, exceeding a Thomson thickness $`\tau _{\mathrm{es}}`$ of few units, to compute correctly the reflected spectrum in the UV range for $`\xi 10^3`$, when dealing with constant density slabs. This has always been overlooked in the past.

### 4.1. Influence of the excited levels

In the first version of the code, non hydrogen-like ions were treated with a rough approximation: the populations of the excited levels were computed separately using a two-level approximation, i.e. interlocking between excited levels was neglected, as well as ionizations from these levels. Moreover only a very small number of levels (at most four) were taken into account. However recombinations onto excited states were taken into account in the ionization equilibrium, and the transfer of these photons was treated in an approximative way, according to an escape probability formula proposed by Canfield & Ricchiazzi (1980). Each recombination was also assumed to produce a resonant photon after cascades.

Hydrogen and hydrogen-like ions were - and are - treated as 6-level atoms. Levels 2s and 2p are treated separately, while full l-mixing is assumed for higher levels. All processes including collisional and radiative ionizations and recombinations are taken into account for each level (cf. Mihalas 1978). In the previous version, recombinations onto levels higher than five were not taken into account, which amounted assuming that the higher levels were in LTE with the continuum. It was true for hydrogen in the conditions in which this code is presently used, but not for other hydrogen-like ions. We have therefore recently added recombinations onto higher levels, and we intend to implement “fictive” levels to take ionizations from upper levels into account as it is done for instance in CLOUDY or XSTAR. Level populations are obtained as usual by matrix inversion.

The rough treatment of non hydrogen-like ions could have several consequences.

* First it should not predict correctly the details of the line spectrum, since it neglects subordinate lines and even important resonance lines.
* Also it can lead to a bad estimation of the energy losses and consequently to a change of the temperature. The recombination continuum corresponds generally to a smaller absorption coefficient than the neglected lines, and therefore can escape more easily from the deep layers. It would correspond to an overestimation of the losses and an underestimation of the temperature. But other parameters, such as the ionization state, play also a crucial role in the energy balance, so the net effect is difficult to predict. Fortunately in these models cooling due to lines is not predominant except in some regions.
* It can lead to important errors in the fractional ion abundances. This effect is “self-cumulating”, as an increased proportion of a given ion leads to a stronger ionization edge, which in turns inhibits further photoionizations. So again the net effect is not predictable.
* Finally recombinations onto excited levels are not balanced by ionizations from these levels. It is not important in conditions prevailing in a photoionized slab far from equilibrium (like for instance in the hot skin of an irradiated accretion disk), as the lowest atom levels are strongly overpopulated with respect to LTE. But this treatment is not valid in the innermost layers of a disk atmosphere, where the density is very high so complete LTE is almost reached.

We have recently introduced interlocked multi-level computations for He-like and Li-like ions, in order to obtain a better description of the spectrum, and to treat cases closer to LTE. However we are for the moment restricted to 8 levels for He-like ions and 5 for Li-like ones, which is not sufficient close to LTE. We take into account recombinations onto the levels higher than 8 for He-like ions (respectively 5 for Li-like), adding them as a contribution to the recombination of the highest level. This treatment is similar to the old one, except that it takes into account many subordinate lines and ionizations from excited levels, and it is valid as far as LTE is not reached for the highest levels and ionizations are dominated by the ground level. Like for hydrogen-like ions we plane to add several other levels, and to mimic the contributions of the higher levels.

We mention also that radiative recombinations are computed differently according to the relative values of $`kT`$ to the photon energy bin, in order to get an accurate frequency dependence.

Figs. 7 to 10 illustrate the influence of the new treatment with respect to the old one, for the reference model (except that the thickness is only $`\tau _{\mathrm{es}}=25`$ instead of 80; we have already mentioned that above this value nothing varies anymore for $`\xi =10^3`$, cf. Fig. 4). In this model the level populations and the ionization degrees are far from LTE, even in the deepest layers. It is quite comparable to a classical photoionized model, where excited levels are populated by radiative recombinations and depopulated by deexcitations onto lower levels. Note that these results are obtained with the new transfer treatment.

Fig. 7 shows the new temperature profile. Except at the surface $`T`$ is increased with respect to the previous computation. The change is quite important, reaching a factor two in the transition layer. Fig. 8 shows the change of fractional abundances of Hydrogen, Helium, and Iron. It is relatively important in the hot skin for HeI, but the abundance of this ion is completely negligible in this region, so it should not induce any change in the emitted spectrum around the ionization edge, i.e. 20 eV. The result of a computation with full He-like and Li-like atoms, but without additional recombinations on the highest level, is also shown. Except precisely for HeI in the hot skin, the results with and without the additional term are quite similar.

Figs. 9 and 10 display the reflected spectrum for this model, computed with the old and the new atomic treatment, with a low resolution in the whole energy range (Fig. 9), and with a higher resolution in a smallest but interesting energy range (Fig. 10).

The spectral distribution does not differ much for the two treatments. This is not the case of the detailed features, as it can be seen on Fig. 10. The flux of the main lines of H-like ions in the reflected spectrum are almost not changed (which is expected since the treatment is the same). For the He-like and Li-like ions, the total flux for all lines is multiplied by a factor 3 or 4. To illustrate this effect, Fig. 11 displays the fluxes of the OVI, OVII, and OVIII lines. Even the resonant lines differ in both treatments, as the energy is distributed among a much larger number of lines: for instance the OVIII line at 665 eV decreases by one order of magnitude with the new treatment. The net increase of line flux is simply due to the redistribution of continuum recombination photons into line photons through cascades, an effect which was only partly taken into account in the previous treatment.

Besides, we have added the result of a computation performed with the complete treatment of He-Li-like ions, when suppressing the recombination rates added on the highest level (cf. Figs. 9 and 10). This suppression does not affect at all the spectral distribution, and only slightly the detailed features of the spectrum.

The next further step is to introduce forbidden lines, which could be important for the cooling of the medium in some places.

Note that all the results given in the following sections are obtained without Li-like and He-like complete atom models.

### 4.2. Influence of the density distribution

#### Constant pressure

:

Up to now, we have presented results only for slabs of constant density. It is a peculiar case, as the pressure is thus varying according to the density and the temperature. When the pressure is imposed and not the density, one is faced to the well-known problem of thermal instability described first in the framework of AGN by Krolik, McKee & Tarter (1981). If the pressure is imposed, three (or even more) solutions to the Gain=Loss equation can exist, and one of these solutions is unstable. The gas has therefore to shift either to the Compton dominated hot solution at a temperature close to the Compton value, or to the cold solution dominated by atomic processes, at a few $`10^4`$ K.

Krolik et al. were considering optically thin media with a uniform temperature, whose value is given by the external irradiating spectrum. Here we are considering a stationary case, where the temperature changes across the slab because the mean intensity $`J_\nu `$ is modified as a function of depth. Since $`J_\nu `$ is absorbed in the soft X-ray band, its spectral distribution becomes more dominated by UV as the depth increases, and such a spectrum does not lead to the S-shape curve characteristic of thermal instability (cf. Krolik et al. 1981). However, even if there is no instability, the temperature decreases very rapidly with increasing depth, as the constancy of the pressure leads to an increase of the density, therefore to an increase of the cooling, and to a further decrease of temperature (note that the temperature is not strictly inversely proportional to the density, as radiation pressure generally exceeds gas pressure).

This is illustrated on Fig. 12 which shows the Gain-Loss function versus the temperature for different layers. The model is a slab with a constant pressure, a density at the surface equal to 1011 cm-3, an ionization parameter equal to 1000. The illuminating spectrum is a power law $`F_\nu \nu ^{0.7}`$, chosen because such a flat spectrum is well known to correspond to an “S-shape” curve with a strong instability. However we see that as the temperature decreases with the depth, there are never multiple solutions of the Gain=Loss equation. Note that it is by no means a general result. In particular, a medium with a more rapid variation of the pressure with the depth, like given by hydrostatic equilibrium, could likely lead to a thermal instability, as the transition to the cold region occurs in a region where the spectrum is less absorbed than in a constant density medium (cf. Nayakshin, Kazanas & Kallman 2000). Also if the UV band is absorbed before the soft X-ray one, which is the case when the ionization parameter is small, it should lead to a thermal instability.

Fig. 13 displays the temperature versus depth for slabs of constant density and slabs of constant pressure, with different ionization parameters. Note that it would be better to speak of “ionizing flux” instead of ionizing parameter in the constant pressure case. Actually the comparison is performed with the same incident fluxes. As expected, and in contrast with the constant density case where there is a smooth decrease of the temperature with increasing depth, the constant pressure case gives a sharp transition between the hot and dilute, and the cold and dense, parts of the slab. Moreover the thickness of the “hot skin” is smaller for the same flux in the constant pressure case, in particular for high illumination. For instance, for $`\xi =10^3`$ the Thomson thickness of the hot skin is about unity in the constant density case, and only 0.4 in the constant pressure case. All these effects have strong consequences on the emission spectrum.

For constant pressure the spectrum displays a narrow UV bump, while the constant density the UV bump is more extended in frequency, as it corresponds to intermediate temperatures which do not exist in the constant pressure case (cf. Fig. 14). The difference between the two cases is particularly apparent for high illumination, as expected. It is interesting to note that the Lyman edge is quite weak in the constant pressure case, since the UV spectrum is closer to a blackbody. This is a point that constant density models have difficulties to account for in AGN.

The distribution of fractional ionic abundances is also completely modified. As an illustration, Fig. 15 displays the fractional abundances of iron versus depth. The most remarkable difference is the presence at $`\tau _{\mathrm{es}}1`$ of low ionization species, such as FeIV or FeV, which are completely absent in the constant density case. One notes also that the extension of highly ionized species is smaller than in the constant density case, as expected. For instance FeXIV is present up to $`\tau _{\mathrm{es}}=0.4`$ for constant pressure, instead of $`1.2`$ for constant density. As a consequence, the Iron K lines produced are different from the constant density, as we shall see later.

#### Hydrostatic equilibrium

:

In order to describe an irradiated accretion disk, we have introduced in collaboration with A. Rosanska and B. Czerny a variation of the density given by the hydrostatic equilibrium law. TITAN treats then only the layer at the surface of the disc where the diffusion approximation does not hold. In this case there is a flux incident on the back side of the slab, equal to the flux generated by viscous dissipation inside the disk. Iterations are needed with the hydrostatic equilibrium equations (for more details see Różańska et al. 2001).

Fig. 16 shows the results for the following case: $`M_{\mathrm{BH}}=10^8`$ M, $`\dot{m}=0.03`$, $`R=10R_{\mathrm{sch}}`$, and an X-ray irradiating powerlaw spectrum, with a spectral index 0.9 in flux, and an X-ray flux $`F_X=3.5\times 10^{13}`$, i.e. equal to half the viscous flux. Note that the temperature decreases abruptly like in the constant pressure case, and that the hot skin is very thin. The temperature reaches that given by viscous dissipation at about $`\tau _{\mathrm{es}}=0.4`$. The reflected spectrum displays a prominent blue bump with no soft X-ray excess, and a strong absorption in the UV range due to low ionization species in the deep layers of the atmosphere. It is important to note that such a spectrum differs strongly from that of an irradiated disc in which the density is assumed to be constant (Ross & Fabian 1993, and subsequent works, cf. Różańska et al. 2001 for a detailed discussion). Our results are similar to those of Nayakshin, Kazanas & Kallman (2000).

## 5. Coupling of TITAN and NOAR

### 5.1. Monte Carlo code NOAR

In order to compute the effect of Compton scattering on the high energy part of the reprocessed radiation and to determine precisely the Compton heating-cooling rate in the energy balance equation, one needs another numerical approach. In this aim Abrassart (1998) developed a Monte Carlo code, NOAR (cf. Dumont et al. 2000). The asset of such an approach is also that it enables to investigate an arbitrary geometry and to determine the angular dependence of the observed spectra. Moreover, it allows to easily extract time variability information.

NOAR takes into account direct and inverse Compton scattering, according to the method proposed by Pozdniakov, Sobol & Sunyaev (1983) and by Gorecki & Wilczewski (1984). Given all the fractional abundances and the temperature in each layer provided by TITAN, it computes the absorption cross sections in each layer. Free-free absorption is taken into account, as well as recombination continua of hydrogen and helium like ions (the spectral distribution of the recombination continuum is determined by the local temperature). The proper yields to account for the competition with the Auger effect are used for the fluorescence lines. Fluorescence of Iron XVII to XXIII is suppressed by resonant trapping (but not in TITAN). For the pseudo fluorescence of Iron XXV and XXVI, a simplifying prescription is adopted which includes the two most probable outcomes of a K shell photoionization, that is: direct recombination on ground level or L-K transition. Other line emission is not included.

One use of NOAR is to provide TITAN with the local Compton gains and losses in each layer. This is necessary, because Compton heating-cooling is dominated by energy losses of photons $`>`$ 20 keV for high values of $`\xi `$. In Dumont et al. (2000) it was shown that given an incident spectrum the dependence of the Compton heating-cooling over $`\xi `$ ratio on the Thomson depth depends only slightly on $`\xi `$ except for high values of $`\xi `$. The Compton heating-cooling rate obtained with NOAR is fitted analytically as a function of $`z`$, which is transferred to TITAN. The spectra computed with TITAN are then completed. (cf. Fig. 17, and also Figs. 14 and 16).

### 5.2. Some results of TITAN and NOAR

TITAN and NOAR lead to similar energy distributions of the reflected spectra in the 1-20 keV range. The reflected spectrum above 10 keV exhibits a “Compton hump” depending on the high energy cut-off of the incident spectrum, which is here 100 keV, and on the ionization parameter. The spectrum in the higher energy range is used for the energy balance.

Comptonisation leads also to line broadening, and when it is important, to smearing of the ionization edges. Direct Compton downscatterings produce a red wing, while inverse Compton upscatterings produce a blue wing. The effect can be comparable to relativistic broadening near the black hole, and can account partly for the broad profiles of the iron line in AGN (Abrassart & Dumont 1998).

Fig. 18 displays (with a relatively high resolution, in order to show detailed features) the spectrum around 7 keV in the constant density case, for different values of $`\xi `$, the other parameters being as in the reference model. The feature near 7 keV is made of a mixture of iron edges and of several iron lines:

\- for $`\xi `$ equal to a few hundreds erg cm s-1, the iron line is dominated by low ionization stages, smaller than iron XVII; the absorption edge is proeminent above 7 keV;

\- for $`\xi `$ equal to $`10^3`$ erg cm s-1, FeXXV and FeXXVI lines are apparent, modified by Compton scattering. The lines are significantly broadened, the broadening is asymmetric, the profile is skewed towards the red. The absorption edge is still important, but extends from 8 to 10 keV, due both to Comptonization and to the influence of several edges;

\- for $`\xi `$ equal to $`10^4`$ erg cm s-1, the line has a large red Compton wing and a weak blue wing, and it is weak. The absorption edge is completely erased.

Fig. 19 displays the iron line complex for the constant pressure case, to be compared to the constant density case. Both a “cold line” at 6.4 keV and FeXXV and XXVI lines are present in the constant pressure case, while only highly ionized species are visible in the constant density case. For constant pressure the cold layers are indeed closer from the surface than in the constant density case, so the photons emitted by these layers are able to leave the medium.

## 6. Foreseen improvements

The coupling of TITAN and NOAR allows to compute the structure and the emission of hot Compton thick irradiated media, by solving consistently the energy balance of the medium. Moreover it allows also to take into account Compton broadening of the lines and of the edges in the X-ray range.

In the previous description of the code (Dumont et al. 2000), we had already recalled the importance of the returning flux, which is often neglected, even for relatively low column densities, and we had compared the line transfer treatment for the lines with the escape probability approximation generally used in these problems. Here we have in addition stressed the influence of the pressure law, by showing the strong differences existing both in the structure and in the emitted spectrum, between a constant density medium, a constant pressure medium, or a medium in hydrostatic equilibrium. We have also shown that the atom models are of fundamental importance as they determine the ionization state, the energy balance, the temperature, and in fine the detailed spectral features which will be soon accessible through X-ray spectroscopy.

All along this paper we have mentioned what improvements we intend to bring to the code.

The Accelerated Lambda Iteration method is used now for the continuum. We will soon implement it for the lines, which are the most difficult to converge. It will in particular allow to treat partial redistribution in the lines.

In collaboration with S. Coupe and M.C. Artru (Coupe et al. 2001) we have begun to take into account subordinate lines for all ions by solving a complete multi-level atom for Li-like and He-like ions, after H-like ions which were already introduced in the previous version of the code. In the future we plan to implement other iso-electronic sequences, and to add to the present ones several other levels. In particular we will try to get a better representation of the highest levels for computations close to LTE. We will take into account forbidden lines, which could be important for the cooling of the medium in some places. We also intend to add to the line transfer an option allowing the UV and soft X-ray photons to escape from the medium through Compton scattering.

Though there are still many improvements to perform, we presently use these codes to model AGN spectra in the UV and X-ray range, in conditions where they are valid.

## References

Abrassart A. & Dumont A.-M. 1998, Proceedings of the First XMM workshop on Science with XMM”, M. Dahlem (ed.)

Abrassart A. 1999, Proceedings of the 32nd Cospar Scientific Assembly, Nagoya, Japan, Advances in Space Research, Vol 25 issue/3-4, p 465

Auer H.L., Paletou F. 1994, A&A 284, 675

Canfield R.C., Ricchiazzi P.J. 1980, ApJ 239, 1036

Cannon C.J. 1973, J. Quant. Spectroscop. Radiat.Transfer 13, 627

Cunto W., Mendoza C., Ochsenbein F., Zeippen C.J. 1993, A&A 275, L5

Coupe S., Dumont A-M., Artru, M-C., Collin S. 2001, in preparation

Collin-Souffrin S., Czerny B., Dumont, A. M. & Zycki P.T. 1996, A&A 314, 393

Czerny B. & Dumont A.-M. 1998, A&A, 338, 386

Dumont A.-M., Abrassart A., Collin S. 2000, A&A 357, 823

Francis P.J., et al. 1991, ApJ 373, 465

Ferland G.J., Korista T., Verner D.A., Ferguson J.W., Kingdom J.B., Verner E.M., 1998, PASP 110, 761

Gorecki A., Wilczewski W. 1984, Acta Astronomica 34, 141

Krolik, J.H., McKee, C.F., Tarter, C. B., 1981, ApJ 249, 422

Kunasz P. & Auer H. L. 1988, J. Quant. Spectroscop. Radiat. Transfer 39, 67

Laor, A., Fiore, F., Elvis, M., Wilkes, B.J., McDowell, J.C. 1997, ApJ 477, 93

Ng K. C. 1974, J. Chem. Phys. 61, 2680

Olson G.L., Auer L.H., Buchler J.R. 1986, J. Quant. Spectroscop. Radiat. Transfer 35, 431

Pozdniakov L.A., Sobol I.M., Sunyaev R.A. 1983, Astrophysics and Space Physics Reviews 2, 189

Ross R.R., Fabian A.C., 1993, MNRAS 261, 74

Różańska A., Czerny B., Dumont A-M. Collin, S. 2001, A&nnnA, submitted

Zheng W., Kriss G.A., Telfer R.C., Grimes J.P., Davidsen A.F. 1997, ApJ 475, 469


format: text

text:
# Comparison of Dirac and microscopic Schrödinger optical model potentials for nucleon-nucleus scattering

## Abstract

Detailed comparisons are made between results of calculations for intermediate energy nucleon–208Pb scattering using optical potentials obtained from Dirac phenomenology and a microscopic Schrödinger model. The two approaches yield quite similar results for all of the observables with each projectile, including spin observables, suggesting the two models are physically equivalent. A new phenomenon in the analyzing power is confirmed in the results for 100 MeV scattering.

preprint: LA-UR-01-1496

The physics of the interaction of a nucleon with the nucleus has traditionally been represented by the optical potential. Once the optical potential is specified, and thus the scattering matrix, all observables may be calculated. For intermediate energies there have been a number of methods of calculating the optical potential, the most successful of which have been the Dirac phenomenology  and a recent microscopic, coordinate space, Schrödinger model . Both approaches have had success in not only predicting differential cross sections, but spin and integral observables as well . Note that these models are fundamentally different. The phenomenological Dirac model is dependent on the fitting of data to determine the parameters in the assumed potentials. The microscopic Schrödinger model is derived from the bare nucleon-nucleon ($`NN`$) potential to obtain an effective interaction in-medium that when folded with a suitable representation of the target density gives the optical potential. While a Schrödinger-like potential may be derived from the phenomenological Dirac one, there is no a priori reason to assume that that nonrelativistic potential resembles in any way the potential obtained from the microscopic model.

Yet no real comparison has been made of the results from the Dirac and microscopic Schrödinger approaches for both projectiles simultaneously. An earlier comparison of the Dirac and phenomenological Schrödinger models showed problems inherent in the latter . The phenomenological Schrödinger model usually assumes a local Woods-Saxon form of the optical potential, whose parameters are determined by fitting a complete set of data. Such an approach generally fails to reproduce the spin observables due to the assumptions made in defining the spin-orbit potential and in the neglect of nonlocality in the potential, arising microscopically from the antisymmetrization of the wave function of the projectile and bound state nucleon in the initial state. The purpose of this Letter is to compare in detail the Dirac phenomenology and microscopic Schrödinger models for both proton and neutron scattering in the energy region in which they have had success. The case we have chosen is that of nucleon elastic scattering from 208Pb at 100–300 MeV. Previous work has been done on this system at this energy range in the Dirac approach , from which a striking phenomenon was identified in the analyzing power at 100 MeV: the proton and neutron analyzing powers are exactly out of phase with each other and the one is nearly the inverse of the other about their mean.

The Dirac phenomenology begins with spherically symmetric complex Lorentz isoscalar-scalar and isoscalar-vector potentials together with the Coulomb potential leading to a Walecka-like scattering model  with a large attractive scalar potential and an almost as large repulsive vector potential. A second-order reduction of the Dirac equation then leads to a Schrödinger-equivalent equation with physically correct effective central and spin-orbit potentials by which the observables are accurately predicted. The spin-orbit term and a Coulomb correction term (accounting for part of the difference between proton and neutron in-medium projectile motion) appear naturally. Projectile and target isospin dependences are treated by introducing corresponding spherically symmetric complex Lorentz isovector-scalar and isovector-vector potentials in a relativistic generalization of the standard Lane model . The resulting Dirac equation is suitable for simultaneous analyses of proton-nucleus and neutron-nucleus scattering data up to several GeV, provided that the parameters of the potential have been determined by least-squares adjustment to existing experimental data . Such was done by Kozack and Madland  for nucleon–208Pb scattering, using energy-independent symmetrized Woods-Saxon form factors, and those results are reproduced herein.

The microscopic Schrödinger approach  does not assume any phenomenological form for the potential nor is it derived from any phenomenological potential. It begins instead with the $`g`$ matrices of the $`NN`$ potential; those $`g`$ matrices are solutions of the Brueckner-Bethe-Goldstone equation in infinite nuclear matter. The Bonn-B $`NN`$ potential  was chosen as the starting point for the calculations presented herein. A local density approximation is used to map the infinite matter solution to the nucleus in question by which an effective $`NN`$ interaction is defined in medium. When folded with an appropriate ground state density of the target, the microscopic optical potential is obtained naturally, incorporating Pauli blocking and density dependences. The (coordinate space) potential contains both direct and exchange parts, with the exchange terms arising from the antisymmetrization of the projectile and bound state nucleon wave functions, and so the potential is fully nonlocal. There are no parameters in the model which must be adjusted after the fact to achieve a reasonable fit: all results are obtained from a single predictive calculation. We use the code DWBA98  to calculate the optical potential and observables. Success has been achieved predicting the observables from proton-nucleus scattering for a number of nuclei across a range of energies . For the present case, the ground state density for 208Pb was obtained from a Skyrme-Hartree-Fock calculation by Brown . As the $`g`$ matrix is defined for all two-body spin and isospin channels, the isospin of the projectile selects the correct components of the matrix to define the appropriate optical potential for that projectile. As such the same $`g`$ matrix defines the optical potentials for both proton and neutron scattering in a natural way. The Coulomb interaction is, of course, also included in the calculations of proton scattering.

In Fig. 1, we compare the differential cross sections for 100 \[1(a)\], 200 \[1(b)\] and 300 MeV \[1(c)\] proton and neutron scattering from 208Pb from the Dirac and microscopic Schrödinger models. (Hereafter, for simplicity, we will use the term Schrödinger model as referring to the microscopic Schrödinger model.) The results for proton and neutron scattering calculated from the Schrödinger model are shown by the solid lines, while those from the Dirac model are portrayed by the dashed lines. We first consider the results for 100 MeV scattering. The cross sections for proton scattering as obtained from both models are very similar in shape and magnitude, and agree quite well with the 98 MeV data of Schwandt et al. . While the similarity between the two models is also observed for neutron scattering, the level of agreement worsens above $`20^{}`$, although both models predict a lack of structure in this angular region compared to the proton scattering results. The results of the Schrödinger model calculations reproduce the 200 MeV data of Hutcheon et al.  well, but not so well for the 300 MeV case. This energy is currently the upper limit for the nonrelativistic microscopic formalism using a bare $`NN`$ potential. The underlying $`NN`$ potential must account for particle production at 300 MeV and higher as each nucleon resonance is excited before agreement in the nucleon-nucleus scattering may be achieved. This may account for the disagreement in magnitude between the two models in the neutron scattering at 300 MeV. Over all energies for both projectiles, however, the agreement in the calculated differential elastic scattering from the two approaches is surprisingly good.

As already stated, Kozack and Madland  observed a striking phenomenon in the analyzing powers for 100 MeV nucleon–208Pb elastic scattering. In that work, they noticed that the proton and neutron scattering analyzing powers were out of phase, and the one almost the inverse of the other about their mean, with the neutrons exhibiting a high polarization above $`20^{}`$ (Fig. 3 of ). (A detailed investigation of this phenomenon is outside the scope of this Letter, but work in progress indicates the dominant element to be the presence or absence of a strong Coulomb field.) We compare the results for the analyzing powers at 100, 200, and 300 MeV from the Schrödinger and Dirac models in Fig. 2. In Fig. 2(a), one observes immediately that the same phenomenon as observed by Kozack and Madland is reproduced by the Schrödinger model, although the peaks are more exaggerated than in the Dirac results. Note that both models reproduce the 98 MeV data of Schwandt et al. . The phenomenon is explicitly illustrated in Fig. 3, which displays the spin observables for 100 MeV scattering. Both the Dirac and Schrödinger models predict the effect. A significant difference between the 100 MeV proton and neutron scattering spin rotations is also observed with both models, though not as dramatic as that between the analyzing powers. Their similarity again gives confidence in the results obtained from both models and we encourage measurements of these spin observables, especially the analyzing powers, to see if this is a real effect. Such a measurement would also serve to assess the models, especially for neutron scattering. We note here that neither model has yet explicitly included Mott-Schwinger scattering, but that would be the next step especially if contrary experimental evidence appears. The 200 and 300 MeV results are less striking in the proton and neutron analyzing power differences and beyond $`15^{}`$ the proton and neutron analyzing powers are similar for both energies. Significant differences are observed only in the forward angle scattering, especially at 300 MeV, where in both models the minimum at $`8^{}`$ in the neutron scattering analyzing power is missing in that for proton scattering. Again, the agreement between the two models, especially in this spin observable, is much better than had been expected. Note that the 200 MeV proton scattering data were used in the global least-squares adjustment to determine the Dirac optical potential at that energy, and so the level of agreement between those data and the Dirac result is to be expected. The Schrödinger 200 MeV proton result is a prediction.

The dramatic difference in the analyzing power at 100 MeV is also seen in the spin rotation function, displayed in Fig. 4. As with the other two observables, the relative features between the proton and neutron scattering results at all three energies are common to both models. However, the variations in the neutron scattering results from the two models are more distinct at 100 MeV, and less so as one increases the energy. At 200 MeV, the agreement between the Dirac model and the data of Hutcheon et al.  again is quite good. While that level of agreement is not obtained with the Schrödinger model that calculation reproduces most of the features in the data. As with the analyzing powers at 200 and 300 MeV, the differences between the proton and neutron spin rotation functions from both models significantly differ only below $`20^{}`$. Thus, the two model approaches are in good agreement for both spin observables for both projectiles in their similarities as well as in their differences.

Fundamentally, our Dirac and Schrödinger approaches are different. That the two models agree so well for the differential cross sections as well as both spin observables for both projectiles not only gives confidence in the calculated results at the three energies considered but also suggests, rather strongly, that the models are physically equivalent over this energy range. Namely, the microscopic Schrödinger model incorporates all the dominant medium modifications in the optical potential without significant approximation by using a realistic ground state density yielding a reasonable specification of all terms in the optical potential, whereas the Dirac approach provides for a natural specification of such terms. The comparison obtained in this work is exemplary of the dilemma in judging the relative merits of relativistic vs. nonrelativistic approaches in intermediate energy proton-nucleus scattering analyses . As Ref.  speculates, the answer throughout the energy regime may lie in QCD-based models of nuclear scattering systems. Such theoretical models and concomitant experiments will require simultaneous treatment of proton and neutron scattering in order to be complete; the latter’s experimental database is currently sadly lacking. Finally, the question remains: is the difference between the proton and neutron analyzing powers at 100 MeV real? Only experiment will provide the final answer to that question.

This work was supported by the United States Department of Energy Contract no. W-7405-ENG-36. It is dedicated to the memory of Richard Kozack.


text:
# 1 Singlet fermion contribution to muon anomalous magnetic moment.

Muon anomalous magnetic moment in models with singlet fermions
Subhash Rajpoot

Department of Physics and Astronomy

California State University

Long Beach, CA 90840, USA

## Abstract

It is shown that a minimal extension of the standard model with leptons that transform as singlets under the $`SU(2)`$ symmetry of weak interactions can explain the recently reported derivation in the measured value of the anomalous magnetic moment of the muon from that expected in the standard model.

Recently the Muon $`(g2)`$ Collaboration  has reported a measurement of the anomalous magnetic moment of the positively charged muon,

$$a_\mu ^{exp}=11659202(14)(6)\pm 16510^{11}.$$
(1)
The value expected in the standard model is

$$a_\mu ^{SM}=11659159.6(6.7)10^{11}$$
(2)
There seems to emerged a discrepancy of 2.6 $`\sigma `$ deviation between theory and the world averaged experimental value

$$a_\mu ^{exp}a_\mu ^{SM}=426\pm 16510^{11}.$$
(3)
Moreover, at 90$`\%`$ confidence level, the measurement indicates a relatively large effect due to ”new physics” beyond the standard model. Since the value tends to lie in the ”positive” rather than in the ”negative” direction, it would appear that the discrepancy can only accommodate a very selective type of new physics. To the contrary, the discrepancy allows for several possible scenarios  to be entertained; supersymmetry , scalar leptoquarks , substructure, muon substructure , new gauge bosons, exotic fermions  and models for neutrino masses. Here, a minimal extension of the standard model is proposed which contains fermions transform as singlets under the $`SU(2)`$ of electroweak interactions.

One way to motivate the extension of the standard model with singlet fermions is the generalized see-saw mechanism , widely discussed in the context of neutrino masses. The generalized see-saw mechanism explains in similar terms the smallness of the masses of the fermions belonging to at least the first two generations. Here, to keep matters simple, we propose to extend the second family with singlet charged leptons $`M_L`$ and $`M_R`$ that are vectorlike , i.e., they transform as

$$M_L(1,1,2),M_R(1,1,2),$$
(4)
under the standard model gauge group $`SU(3)\times SU(2)\times U(1)`$. The model still contains the conventional quarks and leptons and a single doublet of Higgs and is free from triangle anomalies. The singlet fermions may have a bare mass term since no symmetry forbids it.

The coupling of the fermions to the muons are given by

$$L=Y_{\mu M}(\overline{\nu }_\mu \overline{\mu })_Li\sigma _2\mathrm{\Phi }^{}M_R+m_{\mu M}\overline{M_L}\mu _R+\text{(H.C)},$$
(5)
where $`Y_{\mu M}`$ and $`m_{\mu M}`$ refer to Yukawa couplings. Without loss of generality we take $`Y_{\mu M}`$ equal to $`m_{\mu M}`$. After spontaneous symmetry breaking, there is left over the usual Higgs boson $`h^0`$ with coupling $`Y_{\mu M}`$ to the $`\mu `$ and $`M`$.

The contribution of the singlet heavy charged lepton $`M`$ to $`a_\mu `$ is due to its couplings to the neutral Higgs boson and the muon as shown in Fig.1. This contribution is given by

$$a_\mu =\frac{Y_{\mu M}^2m_\mu ^2}{8\pi ^2}_0^1\frac{x^2(1x)+x^2m_M/m_\mu }{x(x1)m_\mu ^2+xm_M^2+(1x)m_h^2}𝑑x.$$
(6)
Here $`a_\mu `$ is defined as the coefficient of the term

$$ie\frac{\sigma _{\mu \nu }q^\nu }{4m_\mu }$$
(7)
in the effective Lagrangian. There are additional contributions with standard model gauge boson exchange but they are known to be small. In evaluating $`a_\mu `$, we take the Higgs boson to be light with mass $`m_H=100`$ GeV, which is compatible with the present bounds from phenomenology. In Fig.2, $`a_\mu `$ is plotted as a function of $`m_M`$ varying between 0.1 TeV and 5 TeV for three different values of $`Y_{\mu M}`$. At 90$`\%`$ confidence level $`a_\mu `$ lies between

$$21.5\times 10^{10}a_\mu 63.7\times 10^{10}.$$
(8)
As can bee seen from Fig.2 these bounds select the value of $`Y_{\mu M}=0.1`$ and the following bounds on the mass of the singly charged lepton

$$3\text{TeV}m_M1\text{TeV}.$$
(9)
The present limits on the mass of singly charged heavy lepton are $`m_M95`$ GeV . Thus, with $`Y_{\mu M}0.1`$ and $`m_M`$ of order 1 TeV, the anomaly can easily be explained. It is to be noted that the addition of the singlet lepton leads to flavor changing processes such as $`\mu e\gamma `$. However, the predicted rate is model dependent and even if we make the conservative estimate of taking $`Y_{\mu M}=Y_{eM}`$, the rate is well below the bounds deduced from present experiments.

To conclude, if the reported anomaly  in the theoretical and experimental values of the muon anomalous magnetic moment of the muon withstands the test of time , then a simple viable explanation can be the standard model with an extended second generation of conventional fermions with weak $`SU(2)`$ singlets charged leptons. The singlet leptons required in the scheme lie within the reach of the near future upcoming accelerators.


text:
# The Formation of Disks in Elliptical Galaxies

## 1 Introduction

The merger hypothesis (Toomre & Toomre, 1972) assumes that elliptical galaxies form by major mergers of spiral galaxies. This is supported by observations of nearby merging systems that seem to evolve into systems with properties comparable to elliptical galaxies (Hibbard & Yun, 1999). The gravitational interaction during the merger process also triggers massive star formation and, probably, the formation of massive central black holes and AGNs (Sanders & Mirabel, 1996; Rigopoulou et al., 1999). On the theoretical side, numerous numerical simulations have shown that mergers lead to spheroidal systems with global properties comparable to elliptical galaxies (Barnes & Hernquist, 1992; Barnes, 1998; Bekki & Shioya, 1997; Naab, Burkert & Hernquist, 1999). During the merging epoch the galaxies are far from dynamical equilibrium. Gas, if added to the simulations, condenses into stars by massive bursts of star formation during the merging process or effectively looses angular momentum and falls into the central parts of the merger remnant, probably feeding a central black hole (Mihos & Hernquist, 1996; Springel, 2000; Bekki, 2000). When the merger remnants have settled into dynamical equilibrium, phase mixing and violent relaxation have erased most information about the initial conditions. However, as violent relaxation is incomplete, the stellar phase space distribution of elliptical galaxies should still show fine structure which affects the isophotal shapes and velocity distributions. This fine structure provides fundamental insight into the formation epoch of elliptical galaxies.

Numerical simulations of collisionless mergers of disk galaxies have confirmed this conclusion. The LOSVDs of merger remnants and the orbit structure of their stars contain information about the initial disk orientations (Heyl, Hernquist & Spergel, 1996). The isophotal shapes and the dynamics of the merger remnants are primarily determined by their initial mass ratios (Naab, Burkert & Hernquist, 1999; Bendo & Barnes, 2000).

In this letter we show that the LOSVDs of most elliptical galaxies contain evidence for a second disk component with roughly 10% to 20% the luminosity of the spheroid and a large scale radius which if is order the effective radius of the spheroid. This disk might have formed by late gas infall from gas-rich progenitors, followed by star formation. Observationally, photometric and kinematic investigations of elliptical galaxies have already indicated the existence of large stellar disks of the same order of magnitude as predicted here, especially in disky, rotating ellipticals (Rix & White, 1990; Scorza et al., 1998). We demonstrate that this component is required in almost all ellipticals, independent of whether they are disky or boxy in order to bring the observations in agreement with the merger scenario of the formation of early-type galaxies.

## 2 The merger models

We use the method described by Hernquist (1993) to construct disk galaxies in dynamical equilibrium. The system of units is: gravitational constant G=1, exponential scale length of the larger disk in the merger h=1 and mass of the larger disk $`M_d=1`$. The disks are exponential with an additional spherical, non-rotating bulge with mass $`M_b=1/3`$, a Hernquist density profile (Hernquist, 1990) and a scale length $`r_b=0.2h`$. The system lives in a spherical pseudo-isothermal halo with a mass $`M_d=5.8`$, cut-off radius $`r_c=10h`$ and core radius $`\gamma =1h`$.

The N-body simulations were performed using a direct summation code with the new special purpose hardware GRAPE-5 (Kawai et al., 2000). The 1:1 merger was calculated adopting in total 400000 particles with each galaxy consisting of 60000 disk particles, 20000 disk particles and 120000 halo particles. For the 3:1 merger the parameters of the more massive galaxy were as described above. The low-mass galaxy contained 1/3 the mass and number of particles of the larger galaxy, with a disk scale length of $`h=\sqrt{1/3}`$, as expected from the Tully-Fisher relation.

For both mass ratios, the galaxies approached each other on nearly parabolic orbits with an initial separation of 30 length units and a pericenter distance of 2 length units. In this letter we focus on one corotating and one counterrotating geometry. We have found that these two impact parameters are representative for a much larger set of simulations that have been performed and that will be discussed in a subsequent paper. The inclinations of the corotating disk models (called model A for the equal mass merger and model C for the 3:1 merger) relative to the orbital plane are $`t_1=30^{}`$ and $`t_2=30^{}`$ with arguments of pericenter of $`\omega _1=30^{}`$ and $`\omega _2=30^{}`$, respectively. The counterrotating disks (called model B for 1:1 and model D for 3:1) have $`t_1=30^{}`$ and $`t_2=150^{}`$, $`\omega _1=30^{}`$ and $`\omega _2=30^{}`$, respectively. In all simulations the merger remnants were allowed to settle into dynamical equilibrium for approximately 10 dynamical timescales after the merger was complete. Then their equilibrium state was analysed.

## 3 LOSVD analysis of the merger remnants

To measure the line-of-sight velocity distribution of a merger remnant we shifted the densest region of every two-dimensional projection to the origin. Then we placed a slit with a width of 0.2 unit lengths along the apparent long axis of each projected remnant. Thereafter we binned all particles falling within each grid cell in velocity along the line-of-sight. The grid spacing was chosen to be 0.15 times the projected half-mass radius. The width of the velocity bins was set to a value of 0.2 for line-of-sight velocities $`v_{\mathrm{los}}`$ in the range $`4v_{\mathrm{los}}4`$. This results in 80 velocity bins over the whole velocity interval. Using the binned velocity data we constructed line-of-sight velocity profiles (LOSVD) for each bin along the grid. Subsequently we parametrized deviations from the Gaussian shape of the velocity profile using Gauss-Hermite basis functions (van der Marel & Franx, 1993; Gerhard, 1993; Bendo & Barnes, 2000). The kinematic parameters of each profile ($`\sigma _{\mathrm{fit}}`$, $`v_{\mathrm{fit}}`$, $`H_3`$, $`H_4`$) were then determined by least squares fitting. The large number of simulated stellar particles ($`>`$ 100000) guaranteed that at least 2000 particles fall within each slit inside one effective radius.

## 4 Comparison with observations

The straight lines in Figure 1 indicate the observed local correlation between $`H_3`$ and $`v/\sigma `$ for a sample of elliptical galaxies in low density environments, published by Bender, Saglia & Gerhard (1994). In all cases $`H_3`$ and $`v/\sigma `$ have opposite signs. The data indicate that all elliptical galaxies have LOSVDs with steep prograde wings and broad retrograde wings. Mehlert et al. (2000) have shown that this result also holds for cluster ellipticals in Coma.

The local correlation between $`H_3`$ and $`v/\sigma `$ for the simulated corotating (A,C) and counterrotating (B,D) 1:1 and 3:1 mergers, respectively, is shown by the dots in Figure 1. Every remnant is analysed as seen from 50 random viewing angles. For the corotating mergers the correlation between $`H_3`$ and $`v/\sigma `$ is almost opposite to the observed one. This is also reflected in a positive effective value for $`H_3`$ (left plot in Figure 2) for almost all projection angles, where $`<H_3>`$ is defined as the mean value between 0.25 and 0.75 effective radii (Bender, Saglia & Gerhard, 1994). Here, a positive $`<H_3>`$ corresponds to an anticorrelation between $`H_3`$ and $`v/\sigma `$, such that $`H_3>0`$ for $`v/\sigma <0`$ and vice versa. We find that $`<H_3>`$ versus $`v/\sigma _0`$ does not follow the observed correlation, indicated by the solid and dashed lines. The profiles have broad prograde wings and narrow retrograde wings. The counterrotating 3:1 merger (D) in Fig. 1 shows a large spread around zero but also does not agree with observations. The only exception is the equal mass merger of counter-rotating disks (B) which leads to a very anisotropic elliptical with no signature of rotation.

Note, that the four simulations which we focus on in this paper have been chosen because they are representative for a much larger set of simulations with different orbital geometries. These four models have been recalculated with high resolution for better statistics.

We conclude that collisionless mergers of disk galaxies in general fail to explain the detailed kinematics of all observed elliptical galaxies (also massive, boxy ones) that exhibit a significant amount of rotation ($`v/\sigma 0.2`$) inside one $`r_{\mathrm{eff}}`$. Further evidence for a possible failure of the collisionless merger picture comes from Cretton et al. (2001) who showed that the kinematical properties of very faint and fast rotating, disky elliptical galaxies can not be explained by collisionless mergers of disk galaxies.

## 5 Theoretical evidence for disks in elliptical galaxies?

One plausible explanation for line-of-sight velocity distributions with negative $`H_3`$ is the superposition of a spheroidal body with a disk-like component (Bender, Saglia & Gerhard, 1994). With a simple experiment one can test if the wrong correlations found in our remnants result from a lack of such a component. We artificially added a thin (scale height $`0.05r_{\mathrm{eff}}`$), cold (velocity dispersion perpendicular to the disk $`\sigma _z=0`$), stellar disk with an exponential surface density profile. The disk was placed in the plane defined by the long and intermediate axes of the main stellar body, rotating in the same direction as the main stellar body. The particles of the disk were assumed to move around the center of the galaxy in centrifugal equilibrium with the gravitational potential arising from the total enclosed mass. No additional random motion was added. Under these simple assumptions only two parameters remain free: the total mass of the disk $`M_d`$ in units of the total luminous mass of the remnant and its scale length $`r_d`$ in units of the projected half-light radius $`r_{\mathrm{eff}}`$. The results for the corotating 3:1 remnant (C) are summarized in the Figures 3 and 4. This case is representative for all our unequal mass merger remnants. In general, we find that disks with small masses or radii do not change the LOSVDs of the stellar component. Disks with masses and radii in the region indicated by the black dots in Figure 4, on the other hand, lead to a significant change in the resulting line profile. The prograde wings become steeper than the retrograde ones in very good agreement with observations. This effect is shown in Figure 3 for a disk with 15% of the spheroid mass and $`r_d=1.25r_{\mathrm{eff}}`$. The influence on the $`<H_3>`$ versus $`v/\sigma _0`$ correlation and its dependence on projection effects is shown on the right hand side of Figure 2. Now the values fall in the observed regime.

If the additional disk becomes too massive, the absolute values for $`H_3`$ are larger than observed and the surface brightness profiles change from de Vaucouleurs to exponential profiles which is again not in agreement with observations. We therefore conclude that the existence of an additional stellar disk component with 10% to 20% the luminosity of the spheroid and a scale length of the order of $`r_{\mathrm{eff}}`$ can explain the observed correlation between $`H_3`$ and $`v/\sigma `$ in elliptical galaxies.

## 6 Discussion and Conclusions

Our simulations indicate that models of pure collisionless mergers in general fail to reproduce LOSVDs of observed elliptical galaxies. All simulated profiles locally have the wrong sign of $`H_3`$ with respect to $`(v/\sigma )`$, compared with the observations. A similar result has been reported by Bendo & Barnes (2000). They, however, did not investigate projection effects. A detailed investigation of local properties of the LOSVDs and projection effects shows that the presence of an additional stellar disk component added to the stellar body after the merger is complete could solve the problem.

One possible origin for such a disk is gas that must have been present during the formation epoch of massive ellipticals and that formed a disk after the major merger event was completed. In order to fit the observed profiles, the amount of gas that settled into the equatorial plane must have been significant. It is puzzling why the gas did not turn into stars prior to its infall into the equatorial plane and how the gas could keep its angular momentum required for the disk to have a scale length equal to the effective radius of the spheroid.

Barnes & Hernquist (1996) do find gaseous disk-like components in their equal mass merger models. Tidal torques lead however to efficient angular momentum loss in the gaseous component resulting in gas infall to the center. Only 20% of the initial gas mass (which was 10% of the initial disk mass) and by this less than 2% of the stellar mass settled into an extended disk-like component. In this case the disk would not be massive enough to change the LOSVD of the system significantly. Equal mass mergers of very gas rich disks (more than 50% gas) could, however, lead to the formation of more massive disk-like components in the end.

First simulations by Naab, & Burkert (2001) have shown that for unequal mass mergers of gas rich galaxies which lead to disky, fast rotating ellipticals due to the stronger centrifugal support, a large fraction of the gas settles into extended disks after the merger is complete, if star formation is suppressed. This seems to be an attractive scenario to explain the existence of large disks in disky elliptical galaxies. The origin of large stellar disks in massive boxy ellipticals with a small but significant amount of rotation that form presumably from equal-mass mergers is, however, still unclear.

We thank Ralf Bender, Hans-Walter Rix and Nicola Cretton for helpful discussions.