text:
# Unitarity and geometrical aspects at low 𝑥

## Introduction

Experimental data obtained at HERA  clearly demonstrated rising behaviour of the structure function $`F_2(x,Q^2)`$ at small $`x`$ which is translated to the rising dependence of the total cross–section $`\sigma _{\gamma ^{}p}^{tot}(W^2,Q^2)`$ on center of mass energy $`W^2`$. This effect is consistent with various dependencies on $`W^2`$ and has been treated in different ways, e. g. as a manifestation of hard BFKL Pomeron , a confirmation of the DGLAP evolution in the perturbative QCD , a transient phenomena, i.e. preasymptotic effects  or as a true asymptotical dependence of the off–mass–shell scattering amplitude . This list is far from being complete and other interpretations can be found, e.g. in the review papers .

It is worth to note here that the essential point in the study of low-$`x`$ dynamics is that the space-time structure of the scattering at small values of $`x`$ involves large distances $`l1/mx`$ on the light–cone  and the region $`x0`$ is sensitive to the nonperturbative contributions. The deep–inelastic scattering in this region turns out to be a coherent process where diffraction plays a major role and nonperturbative models such as Regge or vector dominance model can be competitive with perturbative QCD successfully applied for description of the experimental data.

The strong experimentally observed rise of $`\sigma _{\gamma ^{}p}^{tot}(W^2,Q^2)`$ when

$$\sigma _{\gamma ^{}p}^{tot}(W^2,Q^2)(W^2)^{\lambda (Q^2)}$$
(1)
with $`\lambda (Q^2)`$ rising with $`Q^2`$ from about $`0.1`$ to about $`0.4`$ was considered to a somewhat extent as a surprising fact on the grounds of our knowledge of the energy dependence of total cross–sections in hadronic interactions, where $`\lambda 0.1`$. The above comparison between photon–induced and hadron–induced interactions is quite legitimate since the photon is demonstrating its hadronlike nature for a long time. The apparent difference between the hadron and virtual photon total cross–section behaviours however has no fundamental meaning since there is no Froissart–Martin bound in the case off–shell particles . Only under some additional assumptions this bound can be applied .

This problem was addressed in  on the basis of unitarity for off–shell scattering starting from the eikonal representation for the scattering amplitude. It was argued that the observed rise of $`\sigma _{\gamma ^{}p}^{tot}(W^2,Q^2)`$ can be considered as a true asymptotic behaviour and extension of the eikonal representation for off-shell particles does not provide limitations for $`\sigma _{\gamma ^{}p}^{tot}(W^2,Q^2)`$ at large $`W^2`$. It was claimed that HERA data can be described by the classical universal (with $`Q^2`$–independent intercepts) Regge trajectories.

In the present paper we treat similar problems on the basis of the off–shell $`U`$–matrix approach to the amplitude unitarization. It is shown that the unitary representation for off–shell particles and the respective extension of the chiral quark model for the $`U`$–matrix can lead to Eq. (1), where the exponent $`\lambda (Q^2)`$ is related to the $`Q^2`$–dependent interaction radius of a virtual (off–shell) constituent quark.

It is to be stressed here the importance of the effective interaction radius concept . The study of the effective interaction radius dependence on the scattering variables seemed very useful for understanding of the relevant dynamics of high energy hadronic reactions . Now on it is widely known that the respective geometrical considerations about interaction provide a deep insight in hadron dynamics and deep–inelastic scattering (cf. ).

## 1 Off–shell scattering in the $`U`$–matrix method

The extension of the $`U`$–matrix unitarization for the off-shell scattering was considered in . It was supposed that the virtual photon fluctuates into the virtual vector meson states which afterwards interact with a hadron. We considered a single effective vector meson field and used the notions $`F^{}(s,t,Q^2)`$, $`F^{}(s,t,Q^2)`$ and $`F(s,t)`$ for the amplitudes when both initial and final mesons are off mass shell, only initial meson is off mass shell and both mesons are on mass shell, correspondingly. The virtualities of initial and final vector mesons were chosen to be equal $`Q^2`$ since we need the amplitude of the forward virtual Compton scattering.

The equations for the amplitudes $`F^{}`$ and $`F^{}`$ have the same structure as the equation for the on–shell amplitude $`F`$ but relate the different amplitudes. In the impact parameter representation ($`s4m^2`$) they can be written as follows

$`F^{}(s,b,Q^2)`$ $`=`$ $`U^{}(s,b,Q^2)+iU^{}(s,b,Q^2)F^{}(s,b,Q^2)`$
$`F^{}(s,b,Q^2)`$ $`=`$ $`U^{}(s,b,Q^2)+iU^{}(s,b,Q^2)F(s,b).`$ (2)
The solutions are

$`F^{}(s,b,Q^2)`$ $`=`$ $`{\displaystyle \frac{U^{}(s,b,Q^2)}{1iU(s,b)}},`$ (3)
$`F^{}(s,b,Q^2)`$ $`=`$ $`{\displaystyle \frac{U^{}(s,b,Q^2)}{1iU(s,b)}}+i{\displaystyle \frac{[U^{}(s,b,Q^2)]^2U^{}(s,b,Q^2)U(s,b)}{1iU(s,b)}}.`$ (4)
We also assumed the following relation

$$[U^{}(s,b,Q^2)]^2U^{}(s,b,Q^2)U(s,b)=0.$$
(5)
which can be identically fulfilled if the following factorization occurs:

$`U^{}(s,b,Q^2)`$ $`=`$ $`\omega ^2(s,b,Q^2)U(s,b)`$
$`U^{}(s,b,Q^2)`$ $`=`$ $`\omega (s,b,Q^2)U(s,b).`$ (6)
Such factorization is valid, e. g. in the Regge model with factorizable residues and the $`Q^2`$–independent trajectory. It is also valid in the off–shell extension of the chiral quark model for the $`U`$–matrix. We consider the latter further in detail.

Thus, we have for the amplitudes $`F^{}`$ and $`F^{}`$

$`F^{}(s,b,Q^2)`$ $`=`$ $`{\displaystyle \frac{U^{}(s,b,Q^2)}{1iU(s,b)}}=\omega (s,b,Q^2){\displaystyle \frac{U(s,b)}{1iU(s,b)}}`$ (7)
$`F^{}(s,b,Q^2)`$ $`=`$ $`{\displaystyle \frac{U^{}(s,b,Q^2)}{1iU(s,b)}}=\omega ^2(s,b,Q^2){\displaystyle \frac{U(s,b)}{1iU(s,b)}}`$ (8)
and unitarity provides inequalities

$`|F^{}(s,b,Q^2)|`$ $``$ $`|\omega (s,b,Q^2)|,`$
$`|F^{}(s,b,Q^2)|`$ $``$ $`|\omega ^2(s,b,Q^2)|.`$ (9)
To discuss the asymptotical behaviour of $`\sigma _{\gamma ^{}p}^{tot}`$ we consider off–shell extension of the model for hadron scattering , which is based on the ideas of chiral quark models. Valence quarks located in the central part of a hadron are supposed to scatter in a quasi-independent way by the effective field. In accordance with the quasi-independence of valence quarks we represent the basic dynamical quantity in the form of product:

$$U(s,b)=\underset{i=1}{\overset{n_{h_1}}{}}f_{Q_i}(s,b)\underset{j=1}{\overset{n_{h_2}}{}}f_{Q_j}(s,b)$$
(10)
in the impact parameter representation, $`N=n_{h_1}+n_{h_2}`$ is the total number of constituent quarks in the initial hadrons. Factors $`f_Q(s,b)`$ correspond to the individual quark scattering amplitude smeared over transverse position of $`Q`$ inside hadron $`h`$ and over fraction of longitudinal momentum of the initial hadron carried by quark $`Q`$. Factorization Eq. (10) reflects the coherence in the valence quark scattering and may be considered as an effective implementation of constituent quarks’ confinement. This mechanism resembles Landshoff mechanism of quark–quark independent scattering . However, in this case we refer not to pair interaction of valence quarks belonging to the colliding hadrons, but rather to Hartree–Fock approximation for the constituent quark scattering in the mean field.

The picture of hadron structure implies that the overlapping of hadron structures and interaction of the condensates occur at the first stage of collision. Due to excitation of the condensates, the quasiparticles, i.e. massive quarks arise. These quarks play role of scatterers. To estimate number of such quarks one could assume that part of hadron energy carried by the outer condensate clouds is being released in the overlap region to generate massive quarks. Then their number can be estimated by the quantity:

$$\stackrel{~}{N}(s,b)\frac{(1k_Q)\sqrt{s}}{m_Q}D_c^{h_1}D_c^{h_2},$$
(11)
where $`m_Q`$ – constituent quark mass, $`k_Q`$ – hadron energy fraction carried by the constituent valence quarks. Function $`D_c^h`$ describes condensate distribution inside the hadron $`h`$, and $`b`$ is an impact parameter of the colliding hadrons $`h_1`$ and $`h_2`$. Thus, $`\stackrel{~}{N}(s,b)`$ quarks appear in addition to $`N`$ valence quarks. Those quarks are transient ones: they are transformed back into the condensates of the final hadrons in elastic scattering. It should be noted that we use subscript $`Q`$ to refer the constituent quark $`Q`$ and the same letter $`Q`$ is used to denote a virtuality $`Q^2`$. However, they enter formulas in a way excluding confusion.

The amplitude $`f_Q(s,b)`$ describes elastic scattering $`QQ`$ of a single valence on-shell quark $`Q`$ in the effective field and we use for the function $`f_Q(s,b)`$ the following expression

$$f_Q(s,b)=[\stackrel{~}{N}(s,b)+(N1)]V_Q(b),$$
(12)
where $`V_Q(b)`$ has a simple form

$$V_Q(b)g\mathrm{exp}(m_Qb/\xi ),$$
(13)
which corresponds to the quark interaction radius

$$r_Q=\xi /m_Q.$$
(14)
This picture can be extended for the case when the hadron $`h_2`$ (vector meson in our case) is off mass shell. The off–shell $`U`$–matrix, i.e. $`U^{}(s,b,Q^2)`$ should be then presented as the product

$$U^{}(s,b,Q^2)=\underset{i=1}{\overset{n_{h_1}}{}}f_{Q_i}(s,b)\underset{j=1}{\overset{n_{h_2^{}}}{}}f_{Q_j^{}}(s,b,Q^2)$$
(15)
and the function $`f_Q^{}(s,b,Q^2)`$ is to be written as

$$f_Q^{}(s,b,Q^2)=[\stackrel{~}{N}(s,b)+(N1)]V_Q^{}(b,Q^2).$$
(16)
The notion $`f_Q^{}(s,b,Q^2)`$ stands for the the smeared amplitude which describes elastic scattering $`Q^{}Q^{}`$ of a single valence constituent off-shell quark $`Q^{}`$ in the effective field. In the above equation

$$V_Q^{}(b,Q^2)g(Q^2)\mathrm{exp}(m_Qb/\xi (Q^2))$$
(17)
and it corresponds to the virtual constituent quark interaction radius

$$r_Q^{}=\xi (Q^2)/m_Q.$$
(18)
Under this we mean constituent quark composing effective virtual meson.

The $`b`$–dependence of $`\stackrel{~}{N}(s,b)`$ is weak compared to the $`b`$–dependence of $`V_Q`$  and therefore we have taken this function to be independent on the impact parameter $`b`$.

Dependence on virtuality $`Q^2`$ comes through dependence of the intensity of the virtual constituent quark interaction $`g(Q^2)`$ and the parameter $`\xi (Q^2)`$, which determines the quark interaction radius (in the on-shell limit $`g(Q^2)g`$ and $`\xi (Q^2)\xi `$).

Following these considerations, the explicit dependencies of the generalized reaction matrices $`U^{}`$ and $`U^{}`$ on $`s`$, $`b`$ and $`Q^2`$ can easily be written in the form of Eq. (6) with

$$\omega (s,b,Q^2)=\frac{f_Q^{}(s,b,Q^2)}{f_Q(s,b)}.$$
(19)
Note that Eqs. (5) and (6) imply that

$$f_{Q^{}Q}(s,b,Q^2)=[f_Q^{}(s,b,Q^2)f_Q(s,b)]^{1/2}.$$

We consider high–energy limit and for simplicity suppose that all constituent quarks have equal masses and parameters $`g`$ and $`\xi `$ as well as $`g(Q^2)`$ and $`\xi (Q^2)`$ do not depend on quark flavor. We also consider for simplicity pure imaginary amplitudes. Then we have for the functions $`U`$, $`U^{}`$ and $`U^{}`$ the following expressions

$`U(s,b)=ig^N\left({\displaystyle \frac{s}{m_Q^2}}\right)^{N/2}\mathrm{exp}\left[{\displaystyle \frac{m_QNb}{\xi }}\right]`$ (20)
$`U^{}(s,b,Q^2)=\omega (b,Q^2)U(s,b)`$ (21)
$`U^{}(s,b,Q^2)=\omega ^2(b,Q^2)U(s,b)`$ (22)
where the function $`\omega `$ is an energy-independent one and has the following dependence on $`b`$ and $`Q^2`$

$$\omega (b,Q^2)=\frac{g(Q^2)}{g}\mathrm{exp}\left[\frac{m_Qb}{\overline{\xi }(Q^2)}\right]$$
(23)
with

$$\overline{\xi }(Q^2)=\frac{\xi \xi (Q^2)}{\xi \xi (Q^2)}.$$
(24)
It is clear that for the on–shell particles $`\omega 1`$ and using Eqs. (3) and (4) we will arrive at large $`W^2`$ to the result obtained in

$$\sigma _{\gamma p}^{tot}(W^2)\frac{\xi ^2}{m_Q^2}\mathrm{ln}^2\frac{W^2}{m_Q^2},$$
(25)
where the usual for deep–inelastic scattering notation $`W^2`$ instead of $`s`$ is used. Similar result is valid also for the off mass shell particles when the interaction radius of virtual quark does not depend on $`Q^2`$ and is equal to the interaction radius of the on–shell quark, i.e. $`\xi (Q^2)\xi `$. The behaviour of the total cross–section at large $`W^2`$

$$\sigma _{\gamma ^{}p}^{tot}(W^2)\left[\frac{g(Q^2)\xi }{gm_Q}\right]^2\mathrm{ln}^2\frac{W^2}{m_Q^2},$$
(26)
corresponds to the result obtained in .

For the off–shell case (and $`\xi (Q^2)>\xi `$) the situation is different and we have at large $`W^2`$

$$\sigma _{\gamma ^{}p}^{tot}(W^2,Q^2)G(Q^2)\left(\frac{W^2}{m_Q^2}\right)^{\lambda (Q^2)}\mathrm{ln}\frac{W^2}{m_Q^2},$$
(27)
where

$$\lambda (Q^2)=\frac{\xi (Q^2)\xi }{\xi (Q^2)}.$$
(28)
We shall further concentrate on this as we currently think the most interesting case.

However, it should be noted that for $`\xi (Q^2)<\xi `$

$$\sigma _{\gamma ^{}p}^{tot}(W^2)\left[\frac{g(Q^2)\xi }{g\lambda (Q^2)m_Q}\right]^2,$$

i. e. asymptotically cross–section would be energy-independent and this option cannot be excluded in principle, since we are dealing with the limit $`W^2\mathrm{}`$. The last scenario would mean that the experimentally observed rise of $`\sigma _{\gamma ^{}p}^{tot}`$ is transient preasymptotic phenomena . It can be realized when we replace in the formula for the interaction radius of the on–shell constituent quark $`r_Q=\xi /m_Q`$ the mass $`m_Q`$ by $`m_Q^{}=\sqrt{m_Q^2+Q^2}`$ in order to obtain the interaction radius of the off-shell constituent quark and write it down as $`r_Q^{}=\xi /m_Q^{}`$, or equivalently replace $`\xi (Q^2)`$ for $`\xi (Q^2)=\xi m_Q/\sqrt{m_Q^2+Q^2}`$.

All the above expressions for $`\sigma _{\gamma ^{}p}^{tot}(W^2)`$ can be rewritten as the corresponding dependencies of $`F_2(x,Q^2)`$ at small $`x`$ according to the relation $`F_2(x,Q^2)=\frac{Q^2}{4\pi ^2\alpha }\sigma _{\gamma ^{}p}^{tot}(W^2)`$ where $`x=Q^2/W^2`$. In particular, the Eq. (27) will appear in the form

$$F_2(x,Q^2)\stackrel{~}{G}(Q^2)\left(\frac{1}{x}\right)^{\lambda (Q^2)}\mathrm{ln}(1/x),$$
(29)
## 2 Phenomenological implications

It is interesting that the value and $`Q^2`$ dependence of the exponent $`\lambda (Q^2)`$ is related to the interaction radius of the virtual constituent quark. The value of parameter $`\xi `$ in the model is determined by the slope of the differential cross–section of elastic scattering at large $`t`$ , i. e.

$$\frac{d\sigma }{dt}\mathrm{exp}\left[\frac{2\pi \xi }{m_QN}\sqrt{t}\right].$$
(30)
and from $`pp`$-experimental data $`\xi =22.5`$. Uncertainty is related to the ambiguity in the constituent quark mass value. Using for simplicity $`\xi =2`$ and the data for $`\lambda (Q^2)`$ obtained at HERA  we calculated “experimental” $`Q^2`$–dependence of the function $`\xi (Q^2)`$:

$$\xi (Q^2)=\frac{\xi }{1\lambda (Q^2)}.$$
(31)
Results are represented in Figure 1. It is clear that experiment leads to $`\xi (Q^2)`$ rising with $`Q^2`$. This rise is slow and consistent with $`\mathrm{ln}Q^2`$ extrapolation. Thus, assuming this dependence to be kept at higher $`Q^2`$ and using Eq. (28), we predict saturation in the $`Q^2`$–dependence of $`\lambda (Q^2)`$, i.e. at large $`Q^2`$ the flattening will take place.

## Conclusion

We considered limitations unitarity provides for the $`\gamma ^{}p`$–total cross-sections and geometrical effects in the model dependence of $`\sigma _{\gamma ^{}p}^{tot}`$. In particular, it was shown that the constituent quark’s interaction radius dependence on $`Q^2`$ can lead to a nontrivial, asymptotical result: the behaviour of $`\sigma _{\gamma ^{}p}^{tot}`$ will be $`\sigma _{\gamma ^{}p}^{tot}(W^2)^{\lambda (Q^2)}`$, where $`\lambda (Q^2)`$ will be saturated at large values of $`Q^2`$. This result is valid when the interaction radius of the virtual constituent quark is rising with virtuality $`Q^2`$. The reason for such rise might be of a dynamical nature and it could originate from the emission of the additional $`q\overline{q}`$–pairs in the nonperturbative structure of a constituent quark. In the present approach constituent quark consists of a current quark and the cloud of quark–antiquark pairs of the different flavors . Available experimental data are consistent with the dependence $`\xi (Q^2)\mathrm{ln}Q^2`$.

We would also like to note that unitarity is the most essential property the asymptotical behaviour of the cross–section $`\sigma _{\gamma ^{}p}^{tot}`$ depends on, and the above power energy dependence is not the only possible one. Unitarity transforms strong energy dependence of the “Born” amplitude $`(W^2)^{N/2}`$ into $`\mathrm{ln}^2W^2`$ dependence for the on–shell particles or when interaction radius of the virtual constituent quark does not depend on virtuality and equal to the interaction radius of the on–shell constituent quark. Unitarity transforms this strong energy behaviour into the one, which depends on energy mildly $`(W^2)^{\lambda (Q^2)}`$ when $`\xi (Q^2)>\xi `$. It can even lead to the asymptotically constant cross–section $`\sigma _{\gamma ^{}p}^{tot}`$ when $`\xi (Q^2)<\xi `$.

The available experimental data for the structure functions at low values of $`x`$ continue to demonstrate the rising total cross-section of $`\gamma ^{}p`$–interactions and therefore we currently consider rising with virtuality interaction radius of a constituent quark as a most relevant option, however, it does not mean that the other possibilities are principally excluded.

## Acknowledgements

We would like to thank V. A. Petrov for the stimulating discussions. This work was supported in part by RFBR under Grant 99-02-17995.


format: text

text:
# Untitled Document

Nature of traps responsible for failure of MOS devices

V.A. Gritsenko, P.M. Lenahana, Yu.N. Morokovb, and Yu.N. Novikov

Institute of Semiconductor Physics, Novosibirsk, 630090, Russia
a The Pennsylvanian State University, University Park, Pennsylvania 16802
b Institute of Computational Technologies, Novosibirsk, 630090, Russia

## Abstract

A failure of chips in a huge amount of modern electronic devices such as computers, TV, radio etc is connected as a rule with the undesirable capturing of charge (electrons and holes) by traps in thin insulating film of silicon oxide (SiO2) in transistors. It leads to a breakdown of transistors or to a destructive change of their characteristics. It is suggested that silicon oxide will be replaced in the next generation of nanoscale devices by silicon oxynitride (SiOxNy) . Therefore, it is very important to understand the nature of traps in this material. We discuss this nature using the quantum-chemical simulation.

More than 90 percent of electronic devices are based on silicon. More than 85% of silicon devices are based on Metal-Oxide-Semiconductor (MOS) transistors. The MOS transistor consists of silicon substrate, source, drain, and gate (Fig. 1). Gate is separated from silicon by SiO2. The traps in the thin film of SiO2 can capture electrons or holes from a channel. This is one of the main source of the failure of electronic devices. A key unsolved problem with SiO2 as a gate dielectric is a hole capture in the Si-Si defects near the Si/SiO2 interface and the following device degradation .

Silicon oxide, which has been used as the gate dielectric more than thirty years, will be replaced in next generation of nanoscale MOS devices by silicon oxynitride . A gate oxynitride fabrication with high-temperature annealing of thermal SiO2 in ammonia results in a hole traps removal from the Si/SiO2 interface  (Fig. 2). However, such oxynitride contains undesirable large density of electron traps. To improve and to control the dielectric properties it is important to understand which structural defects are the dominant traps in oxide and oxynitride. In the present paper we consider, using the quantum-chemical simulation, the nature of such traps in SiOxNy and mechanism of creation and removal of traps.

The defects and traps are absent in ”ideal” oxynitride containing only the Si-O and Si-N bonds. The oxygen and nitrogen atoms in amorphous silicon oxynitride as in SiO2 and Si3N4 (silicon nitride) are bonded with two and three silicon atoms, respectively . The Si atoms in these materials are bonded with four other atoms (O or N). The traps in SiOxNy are related with defects in atomic network.

The main defects, identified in SiOxNy, are the next: the silicon-silicon bond ($``$Si-Si$``$), the three-fold coordinated silicon atom $``$Si\*, the two-fold coordinated nitrogen atom $``$Si2N\*, and the non-bridged oxygen atom $``$SiO\*. Here ”\*” marks one unpaired electron. The main defects related with impurities are the hydrogen bonds: $``$SiOH, $``$Si2NH, and $``$SiH.

To understand the nature of traps in SiOxNy we considered the capturing properties of the main defects in SiOxNy. We have studied the electronic structure of the clusters (small fragments of material) in different charge states corresponding to the capture of electrons or holes by traps. The dangling bonds at the cluster surface were saturated with hydrogen atoms. Atomic relaxation of the atoms in the first coordination sphere of defects was considered in the simulation. The calculations were maid with the quantum-chemical method MINDO/3 with the parameters used in Refs. .

We used the Si18N17H35 and Si20N28H36 clusters for simulation of the Si3N4 bulk electronic structure. The clusters with different numbers of oxygen and nitrogen atoms in the second coordination sphere were considered to study the effect of the chemical composition on the capturing properties of the defects in SiOxNy. We estimated the energy gain for the electron or hole capture from the bulk of material as the differences between the total energies of the bulk and defect clusters in different charge states.

Simulation of the neutral $``$Si2N\* defect (nitrogen dangling bond)  shows that unpaired electron is localized for all considered clusters in a nitrogen 2pπ orbital, which is oriented normally to the Si2N plane. This result agrees with the previous theoretical simulation of this defect  and with the ESR experiments  for the $``$Si2N\* defect in Si3N4. The positive energy gain for electron capture by the $``$Si2N\* defect in Si3N4 is about 0.8 eV. It means that the defect can capture an electron. The calculated results show that the $``$Si2N\* defect in silicon nitride and in oxynitride with a high concentration of nitrogen can capture an electron and cannot capture a hole .

The analogue of the $``$Si2N\* defect in Si3N4 is the $``$SiO\* defect in SiO2. It was shown in Ref.  that the $``$SiO\* defect is the electron trap in SiO2. Our calculations for the $``$SiO\* defect in SiOxNy and SiO2 confirm this result. Thus, the obtained results together with experiments show that the $``$Si2N\* and $``$SiO\* defects are the electron traps in gate dielectric.

The $``$Si2N\* defect creation in SiOxNy is related with a breaking of the $``$Si2N-H bond according to the reaction

$``$Si2NH $``$ $``$Si2N\* + H.

We associate the electron traps in oxynitride with these $``$ Si2N\* defects. In order to remove the electron traps from gate oxynitride in silicon devices it is used the reoxidation of oxynitride. The removal of the electron traps at oxynitride reoxidation is explained by replacing of the =N\* atom with oxygen atom according to the reaction

$``$Si2N\* + O $``$ $``$Si-O-Si$``$ \+ N.

With this reaction we can understand the chemical nature of one of the key technological step in the MOSFET device fabrication - the electron traps removal at the gate silicon oxynitride reoxidation (Fig. 2).

The calculations have shown that the main extrinsic defects, related with the hydrogen bonds, $``$SiOH, $``$Si2NH, and $``$SiH did not capture electrons and holes.

Our simulation shows that the Si-Si bonds in Si3N4 and SiO2 can be both the electron and hole traps. We have obtained also that the $``$Si\* defect in Si3N4 and SiO2 can capture an electron. However, the ESR experiments do not show the existence of this defect in Si3N4 and SiO2 .

The recent density functional calculations of the $``$Si\* and $``$ Si2N\* defects in silicon nitride  confirm that both defects can capture an electron.

Thus, the calculations of the electronic structure of the main intrinsic defects in SiOxNy in different charge states support the identification of the $``$Si2N\* and $``$SiO\* defects with the electron traps. The hole traps are related with the Si-Si bonds. This conclusion is in agreement with the results of experiments . It was suggested earlier by Robertson and Powell  that the negatively charged nitrogen defect $``$Si2N\* can be a hole trap. Our results agree with it and show also that the electron localization by neutral defect $``$Si2N\* leads to the ESR signal disappearance. This effect was experimentally observed earlier . The results of calculations allow to explain the mechanism of creation and removal of the intrinsic defects and to understand the chemistry of the main technological steps in the gate silicon oxynitride fabrication. With these results the nature of the traps responsible for the electron and hole capture in gate oxynitride and so the main source of the silicon devices failure become more clear.

References

Electronic address: grits@isp.nsc.ru
1. E.P. Gusev, H.C. Lu, E.L. Garfunkel, T. Gustafsson, and M.L. Green, IBM J. Research and Develop. 43 (1999) 1.
2. P.M. Lenahan and P.V. Dressendorfer, J. Appl. Phys. 55 (1984) 3495.
3. V.A. Gritsenko, J.B. Xu, Y.H. Ng, R.W.M. Kwok, and I.H.Wilson, Phys. Rev. Lett. 81 (1998) 1054.
4. I.A. Britov, V.A. Gritsenko, and Yu.N. Romaschenko, Sov. Phys. JETP 62 (1985) 321.
5. V.A. Gritsenko, Yu.N. Morokov, and Yu.N. Novikov, Appl. Surf. Sci. 113/114 (1997) 417.
6. Yu.N. Morokov, V.A. Gritsenko, Yu.N. Novikov, and H. Wong, Microelectronic Engineering 48 (1999) 175.
7. W.L. Warren, J. Robertson, and J. Kanicki, Appl. Phys. Lett. 63 (1993) 2685.
8. L. Martin-Moreno, E. Martinez, J.A. Verges, and F. Yndurain, Phys. Rev. B35 (1987) 9683.
9. W.L. Warren, P.M. Lenahan, and S.E. Curry, Phys. Rev. Lett. 65 (1990) 207.
10. J.T. Yount and P.M. Lenahan, J. Non-Cryst. Solids 164-166 (1993) 1069.
11. A.H. Edwars and W.B. Fowler, in Structure and Bonding in Noncrystalline Solids, G.E. Walrafen and A.G. Revesz, Eds. (Plenum Press, New York and London, 1986) 139-155.
12. V.A. Gritsenko and A. D. Milov, JETP Lett. 64 (1996) 531.
13. G. Pacchioni and D. Erbetta, Phys. Rev. B61 (2000) 15005.
14. V.A. Gritsenko, I.P. Petrenko, and S.N. Svitasheva, Appl. Phys. Lett. 72 (1998) 462.
15. J. Robertson and M.J. Powell, Appl. Phys. Lett. 44 (1984) 415.
16. I.A. Chaiyasena, P.M. Lenahan, and G.J. Dunn, Appl. Phys. Lett. 58 (1991) 2141.


text:
# Emission from dust in galaxies: metallicity dependence

## 1 Introduction

When we want to know the entire evolutionary history of a galaxy, its color and metallicity are important diagnostic quantities (e.g., Tinsley tinsley72 (1972)). The evolution of the color and the metallicity results from the superposition of successive star formation during the lifetime of the galaxy. Thus, if we wish to reveal the evolution of galaxies, we always need to estimate the star formation rate (SFR) on a galaxy-wide scale. There are many methods to estimate the SFR of galaxies from observational quantities (Kennicutt 1998a ). In this paper, we are especially interested in the conversion formula from infrared (IR) luminosity of the galaxies to their present SFR. The IR luminosity originates from dust that is heated by stellar radiation field. Kennicutt (1998b ) estimated SFR from IR luminosity for starburst galaxies. Buat & Xu (buat96 (1996)) adopted an empirical approach by utilizing the observed relation between the ultraviolet (UV) luminosity and the IR luminosity.

Recently, Inoue et al. (2000a , hereafter IHK00) have derived a theoretical conversion formula from IR luminosity to the SFR by developing a standard model of H ii regions by Petrosian et al. (petrosian72 (1972)). The present paper examines how metallicity (or dust-to-gas ratio) affects the conversion factor of IHK00. It is possible because the effect of the metallicity is well parameterized in their formula. Here, we present it as

$`{\displaystyle \frac{\mathrm{SFR}}{M_{}\mathrm{yr}^1}}={\displaystyle \frac{3.3\times 10^{10}(1\eta )}{0.40.2f+0.6ϵ}}{\displaystyle \frac{L_{\mathrm{IR}}}{L_{}}},`$ (1)
where $`f`$ is the fraction of ionizing photons absorbed by hydrogen, $`ϵ`$ is the efficiency of dust absorption for nonionizing photons, $`L_{\mathrm{IR}}`$ is the observed luminosity of dust in the wavelength range of 8–1000 $`\mu `$m, and $`\eta `$ is the cirrus fraction of $`L_{\mathrm{IR}}`$. In Eq. (1), $`f`$ and $`ϵ`$ depend on the dust-to-gas ratio111In this paper, we focus our main attention on only $`f`$ and $`ϵ`$, although we should keep in mind that the ($`1\eta `$) may depend on the dust-to-gas ratio. The effect of ($`1\eta `$) on SFR is discussed in Appendix A.. According to Hirashita (1999a, b), the dust-to-gas ratio is expressed as a function of metallicity (see also Dwek dwek98 (1998); Lisenfeld & Ferrara lisenfeld98 (1998)). Therefore, we find easily that the relation between $`L_{\mathrm{IR}}`$ and SFR depends on metallicity via $`f`$ and $`ϵ`$.

In this paper, we quantitatively examine the importance of the metallicity to our conversion law, the conclusion being that the metallicity dependence is not always negligible when we compare SFRs of very young galaxies and those of present galaxies. First, we examine the formula by IHK00 in Sect. 2, focusing on the effect dust-to-gas ratio. Then, in Sect. 3, we present the metallicity dependence of the IHK00’s formula. In Sect. 4, we extend our discussion to the cosmic star formation history. Finally, we summarize our conclusions in Sect. 5 .

## 2 Dependence on dust-to-gas ratio

In this section, we derive the metallicity dependence of the conversion formula (Eq. (1)). It is convenient to define the factor $`C_{\mathrm{IR}}`$ as

$`\mathrm{SFR}`$ $`=`$ $`C_{\mathrm{IR}}L_{\mathrm{IR}},`$ (2)
$`C_{\mathrm{IR}}`$ $``$ $`{\displaystyle \frac{3.3\times 10^{10}(1\eta )}{0.40.2f+0.6ϵ}}[M_{}\mathrm{yr}^1L_{}^1].`$ (3)
We focus on the dependence of $`C_{\mathrm{IR}}`$ on metallicity. The dependence of $`C_{\mathrm{IR}}`$ on dust-to-gas ratio is included through $`f`$ and $`ϵ`$ as will be described in Sects. 2.1 and 2.2. We will check our formula observationally in Sect. 2.6. Then, using the relation between the dust-to-gas ratio and the metallicity, we will obtain the metallicity dependence of $`C_{\mathrm{IR}}`$ in the next section. We assume Case B (Baker & Menzel baker38 (1938)), where the optical depths for the Lyman series are large enough.

### 2.1 Dependence of $`f`$ on dust-to-gas ratio

The dependence of $`f`$ (Eq. (3)) on dust-to-gas ratio is obtained from Spitzer (spitzer78 (1978), hereafter S78). We define $`\tau _{\mathrm{S},\mathrm{d}}`$ and $`\tau _{\mathrm{S},\mathrm{H}}`$ as the optical depths of the dust and of the neutral hydrogen atoms, respectively, for the Lyman-continuum photons over a path length equal to the Strömgren radius, $`r_\mathrm{S}`$ (see S78 or Eq. (6) for the definition of $`r_\mathrm{S}`$). Since the radius of a dust-free H ii region is estimated to be $`r_\mathrm{S}`$, it is useful to normalize the length scale by $`r_\mathrm{S}`$. We assume spherical symmetry and spatial uniformity of H ii regions in this paper for the simplicity.

First, we calculate $`\tau _{\mathrm{S},\mathrm{d}}/\tau _{\mathrm{S},\mathrm{H}}`$ for an H ii region. Here, both $`\tau _{\mathrm{S},\mathrm{d}}`$ and $`\tau _{\mathrm{S},\mathrm{H}}`$ are approximated with the optical thickness at the Lyman limit (912 Å). Then we obtain

$`{\displaystyle \frac{\tau _{\mathrm{S},\mathrm{d}}}{\tau _{\mathrm{S},\mathrm{H}}}}={\displaystyle \frac{\tau _{\mathrm{S},\mathrm{d}}}{n_\mathrm{H}sr_\mathrm{S}}}=12{\displaystyle \frac{E_{BV}}{N_\mathrm{H}s}},`$ (4)
where the dust extinction at 912 Å, $`A_{912}`$, is taken to be $`13E_{BV}`$ mag according to the Galactic extinction curve, and $`\tau _{\mathrm{S},\mathrm{d}}=A_{912}/(2.5\mathrm{log}_{10}e)`$. Moreover, $`n_\mathrm{H}`$ and $`s=6.30\times 10^{18}\mathrm{cm}^2`$ denote the number density of the hydrogen and the absorption cross section for a hydrogen atom in the $`n=1`$ level (Eq. (5.6) of S78), respectively, and the column density $`N_\mathrm{H}`$ is defined as $`N_\mathrm{H}n_\mathrm{H}r_\mathrm{S}`$.

If the physical properties of grains such as the extinction curve are unchanged, $`E_{BV}/N_\mathrm{H}`$ is proportional to the dust-to-gas mass ratio $`𝒟`$. According to S78, $`𝒟=6\times 10^3`$ when $`N_\mathrm{H}/E_{BV}=5.9\times 10^{21}\mathrm{mag}^1\mathrm{cm}^2`$ (Sects. 7.2 and 7.3 of S78). Here, we apply these values to H ii regions. That is, we assume that the dust-to-gas ratio in the H ii region is the same as that in the mean value in the interstellar space. Then, Eq. (4) reduces to

$`{\displaystyle \frac{\tau _{\mathrm{S},\mathrm{d}}}{\tau _{\mathrm{S},\mathrm{H}}}}={\displaystyle \frac{1}{3100}}\left({\displaystyle \frac{𝒟}{6\times 10^3}}\right).`$ (5)
This is the same as Eq. (5.23) in S78 but an explicit expression of the dependence on $`𝒟`$.

Next, we calculate $`\tau _{\mathrm{S},\mathrm{H}}=n_\mathrm{H}sr_\mathrm{S}`$. The Strömgren radius is determined as

$`r_\mathrm{S}=1.4\left({\displaystyle \frac{N_\mathrm{u}}{10^{48}\mathrm{s}^1}}\right)^{1/3}\left({\displaystyle \frac{n_\mathrm{H}}{10^2\mathrm{cm}^3}}\right)^{2/3}\mathrm{pc},`$ (6)
where $`N_\mathrm{u}`$ represents the number of ionizing photons emitted from central stars per second. In this equation, we have assumed that the temperature of the H ii region is 8000 K (i.e., the recombination coefficient to the $`n=2`$ level is $`\alpha ^{(2)}=3.09\times 10^{13}\mathrm{cm}^3\mathrm{s}^1`$; see Eq. (5.14) and Table 5.2 of S78) and that $`n_\mathrm{e}=n_\mathrm{H}`$, where $`n_\mathrm{e}`$ is the number density of electrons, is satisfied within the Strömgren radius. Using the estimation of the Strömgren radius above, we obtain

$`\tau _{\mathrm{S},\mathrm{H}}=2.7\times 10^3\left({\displaystyle \frac{n_\mathrm{H}}{10^2\mathrm{cm}^3}}\right)^{1/3}\left({\displaystyle \frac{N_\mathrm{u}}{10^{48}\mathrm{s}^1}}\right)^{1/3},`$ (7)
Combining this with Eq. (5), $`\tau _{\mathrm{S},\mathrm{d}}`$ is estimated as

$`\tau _{\mathrm{S},\mathrm{d}}`$ $`=`$ $`0.87\left({\displaystyle \frac{𝒟}{6\times 10^3}}\right)\left({\displaystyle \frac{n_\mathrm{H}}{10^2\mathrm{cm}^3}}\right)^{1/3}`$ (8)
$`\times \left({\displaystyle \frac{N_\mathrm{u}}{10^{48}\mathrm{s}^1}}\right)^{1/3}.`$
We note that $`\tau _{\mathrm{S},\mathrm{d}}`$ becomes large as $`N_\mathrm{u}`$ increases. This is because the probability that the dust grains absorb the photons inside an H ii region increases as the size of the region becomes larger.

Next, we estimate the fraction of the ionizing photons absorbed by dust grains. Due to the grain absorption, the size of an H ii region is smaller than $`r_\mathrm{S}`$. In other words, if we define $`y_\mathrm{i}`$ as the ionization radius normalized by $`r_\mathrm{S}`$, $`y_\mathrm{i}<1`$. (Without dust grains, $`y_\mathrm{i}=1`$.) A useful relation between $`\tau _{\mathrm{S},\mathrm{d}}`$ and $`y_\mathrm{i}`$ is given in Table 5.4 of S78, where $`y_\mathrm{i}`$ is determined from the following expression (Eq. (5.29) of S78):

$`3{\displaystyle _0^{y_\mathrm{i}}}y^2e^{y\tau _{\mathrm{S},\mathrm{d}}}dy=1.`$ (9)
Using $`y_\mathrm{i}`$, the fraction of the ionizing photons absorbed by hydrogen, $`f`$, is expressed as

$`f=y_\mathrm{i}^3.`$ (10)
This is the same as $`f`$ in Eq. (1). In Fig. 1a, we show $`f`$ as a function of $`\tau _{\mathrm{S},\mathrm{d}}`$.

When we consider the dependence of $`f`$ on the dust-to-gas ratio, a large ambiguity exists: The number of ionizing photons per H ii region is unknown, since the typical number and mass function of OB stars in an H ii region is difficult to determine exactly. Fortunately, for the purpose of finding the dependence on $`𝒟`$, this is resolved by calibrating the “Galactic” $`f`$ with the value of Orion Nebula (Petrosian et al. 1972). According to them, $`f=0.26`$ in the Nebula. They also commented that the value explains the IR emission from H ii regions. Adopting $`f=0.26`$ as the typical value of the Galaxy (see also Inoue et al. 2000b ), we obtain $`\tau _{\mathrm{S},\mathrm{d}}=2.7`$ for the typical Galactic H ii regions (we consider $`𝒟6\times 10^3`$ for such regions) from Fig. 1a. Hence we write

$`\tau _{\mathrm{S},\mathrm{d}}=2.7\left({\displaystyle \frac{𝒟}{6\times 10^3}}\right).`$ (11)
This is consistent with Eq. (8) if we assume $`n_\mathrm{H}=10^2\mathrm{cm}^3`$ and $`N_\mathrm{u}3.0\times 10^{49}\mathrm{s}^1`$. Thus, the net effects of mass function and number of OB stars are included in the value of $`N_\mathrm{u}`$ by adopting Eq. (11). This simplicity is meaningful for our motivation to find the dependence of $`C_{\mathrm{IR}}`$ on $`𝒟`$. Combining Eq. (11) with Eq. (9), we obtain $`f`$ as a function of $`𝒟`$ as shown in Fig. 1b.

### 2.2 Dependence of $`ϵ`$ on dust-to-gas ratio

In IHK00, $`ϵ`$ in Eq. (3) is defined by

$`ϵ1e^{\tau _{\mathrm{nonion}}},`$ (12)
where $`\tau _{\mathrm{nonion}}`$ is the mean optical depth of dust for nonionizing photons. That is, $`ϵ`$ represents the efficiency of the dust absorption of nonionizing photons. Inoue et al. (2000c ) estimated $`ϵ=0.9`$ ($`\tau _{\mathrm{nonion}}=2.3`$) from the averaged visual extinction of Usui et al. (usui98 (1998))’s sample ($`A_V1`$ mag) and the Galactic extinction curve between 1000 Å and 4000 Å by Savage & Mathis (savage79 (1979)). This wavelength range is fit for our purpose, since most of the nonionizing photons from OB stars are emitted in the wavelength range much shorter than 4000 Å. The efficient extinction in this short wavelength makes $`ϵ`$ large (i.e., nearly unity).

It is obvious that $`\tau _{\mathrm{nonion}}`$ depends on the dust-to-gas ratio. If the column density of gas contributing to the absorption of nonionizing photons is fixed, $`\tau _{\mathrm{nonion}}`$ is proportional to the dust-to-gas ratio. Since we are interested in the dependence of $`C_{\mathrm{IR}}`$ on $`𝒟`$, we simply adopt $`\tau _{\mathrm{nonion}}𝒟`$. Here, we determine the numerical value of $`\tau _{\mathrm{nonion}}`$ as

$`\tau _{\mathrm{nonion}}=2.3\left({\displaystyle \frac{𝒟}{6\times 10^3}}\right),`$ (13)
so that $`\tau _{\mathrm{nonion}}`$ becomes 2.3 for the Galactic dust-to-gas ratio. By combining Eqs. (12) and (13), we obtain $`ϵ`$ as a function of $`𝒟`$ as shown in Fig. 2.

### 2.3 Cirrus Fraction

The fraction of the cirrus component, $`\eta `$, remains to be determined. In this paper, an empirical value of $`\eta `$ is simply adopted. According to Lonsdale Persson & Helou (lonsdale87 (1987)), $`\eta 0.5`$–0.7 for their sample spiral galaxies. We use $`\eta =0.5`$ as adopted in IHK00 (the value averaged for Usui et al. 1998’s sample) for the current estimate for our SFR or $`C_{\mathrm{IR}}`$. As a first step, we adopt a constant $`\eta `$, since a theoretically reasonable determination of it is very difficult. The variation of the cirrus fraction by the change of $`𝒟`$ will be briefly considered in Appendix A. We note that the range of $`C_{\mathrm{IR}}`$ becomes larger if we consider the metallicity dependence of the cirrus fraction (Appendix A).

### 2.4 Absorption of Lyman-$`\alpha `$ photons by dust

Since IHK00’s derivation of Eq. (1) is based on Case B, Ly$`\alpha `$ photons are assumed to be easily trapped in an H ii region (Osterbrock osterbrock89 (1989)). Thus, during the resonant scatterings in an H ii region, all the Ly$`\alpha `$ photons are assumed to be absorbed by grains in IHK00, which this paper is based on.

If the dust-to-gas ratio is significantly smaller than the Galactic value, the dust grains might not efficiently absorb the Ly$`\alpha `$ photons in H ii regions. However, H i envelopes on a galactic scale generally exists around H ii regions. A Ly$`\alpha `$ photon is absorbed in the H i envelopes soon after it escapes from an H ii region (Osterbrock osterbrock61 (1961)). Hence, even if there is only a small amount of dust, we expect a sufficient chance for dust to absorb the Ly$`\alpha `$ photons in the H i envelope. We check this point below.

First of all, we define a path length of the Ly$`\alpha `$ photons as $`l_{\mathrm{path}}`$. It is estimated as being $`l_{\mathrm{path}}\tau _{\mathrm{Ly}\alpha }^2l_{\mathrm{Ly}\alpha }`$, where $`\tau _{\mathrm{Ly}\alpha }`$ means the optical depth for Ly$`\alpha `$ photons and $`l_{\mathrm{Ly}\alpha }`$ is the mean free path of the Ly$`\alpha `$ photons. The square to $`\tau _{\mathrm{Ly}\alpha }`$ means that the resonant scattering is assumed to be a random-walk process of photons. Next, we define an optical depth of dust grains for Ly$`\alpha `$ photons over the length of $`l_{\mathrm{path}}`$ as $`\tau _{\mathrm{dust}}`$. It is estimated as $`\tau _{\mathrm{dust}}\pi a^2l_{\mathrm{path}}n_{\mathrm{dust}}`$, where $`a`$ is the size of grains and $`n_{\mathrm{dust}}`$ is their number density. Here, we will discuss whether the Ly$`\alpha `$ photons have a chance to escape from galaxies. We are interested in a star-forming region surrounded by H i gas envelopes. The scale length of the H i envelope, $`L`$, may be estimated to be about 100 pc, which corresponds to the thickness of the disk of spiral galaxies. In such a case, $`\tau _{\mathrm{Ly}\alpha }`$ is estimated to be about $`L/l_{\mathrm{Ly}\alpha }`$. Then, we find $`\tau _{\mathrm{dust}}\pi a^2n_{\mathrm{dust}}L^2/l_{\mathrm{Ly}\alpha }`$. Using $`𝒟`$, it is expressed as

$`\tau _{\mathrm{dust}}{\displaystyle \frac{\pi a^2n_\mathrm{H}𝒟L^2m_\mathrm{p}}{l_{\mathrm{Ly}\alpha }m_{\mathrm{dust}}}},`$ (14)
where $`m_\mathrm{p}`$ is the proton mass and $`m_{\mathrm{dust}}`$ is the dust mass for an assumed spherical dust with mass density of 3 g cm-3. Adopting $`n_\mathrm{H}=0.1`$ cm-3 for a diffuse H i medium, $`a=0.1`$ $`\mu `$m, $`𝒟=6\times 10^5`$ (0.01 times smaller than the Galactic value), $`L=100`$ pc, and the resonant-scattering cross section of Ly$`\alpha `$ photons $`2\times 10^{13}`$ cm2 (i.e., $`l_{\mathrm{Ly}\alpha }=5\times 10^{13}`$ cm), we find $`\tau _{\mathrm{dust}}5\times 10^2`$. Thus, $`\tau _{\mathrm{dust}}`$ is much larger than unity. This means that most of the Ly$`\alpha `$ photons are absorbed even if the dust-to-gas ratio is as small as 0.01 times the Galactic value.

### 2.5 Dependence of $`C_{\mathrm{IR}}`$ on dust-to-gas ratio

In the above subsections, we have expressed $`f`$ and $`ϵ`$ as a function of the dust-to-gas ratio ($`𝒟`$), while $`\eta `$ is treated as a constant ($`\eta =0.5`$). We also assume a typical star-forming region whose mean density of gas and production rate of ionizing photons are about $`10^2`$ cm-3 and $`3.0\times 10^{49}`$ s-1, respectively. Then, we can express $`C_{\mathrm{IR}}`$ defined in Eq. (3) as a function of the dust-to-gas ratio. In Fig. 3, we present $`C_{\mathrm{IR}}`$ as a function of $`𝒟`$. From this figure, we find that the coefficient of the conversion from IR light to SFR becomes 4–5 times smaller for $`𝒟6\times 10^3`$ (the Galactic value) than that for $`𝒟6\times 10^5`$ (0.01 times the Galactic value).

$`C_{\mathrm{IR}}`$ is also changed if we adopt a different extinction curve. In the previous discussions, we have adopted the extinction curve of the Galaxy. The current paragraph examines how the relation between $`C_{\mathrm{IR}}`$ and $`𝒟`$ changes when we adopt the extinction curve of the Small Magellanic Cloud (SMC) with the other quantities adopted in this paper unchanged. In the SMC-type extinction law, the ratio of the optical depth in Eq. (4) becomes two times larger, because the UV extinction is enhanced in the SMC extinction. This indicates that the same estimation performed in Sect. 2 is possible if we make the normalization of $`𝒟`$ half. Thus, we can find a rough dependence of $`f`$ on $`𝒟`$ from Fig. 1b by doubling $`𝒟`$ (i.e., the line in Fig. 1b is moved to the left by 0.3 dex). Since $`𝒟`$ of SMC is about $`5\times 10^4`$ (Issa et al. issa90 (1990)), we find $`f0.75`$ for the SMC extinction case, while $`f0.85`$ if we adopt the extinction law of the Galaxy. We can apply the same scaling to $`ϵ`$. For $`𝒟5\times 10^4`$, $`ϵ0.4`$ for the SMC extinction case, while $`ϵ0.2`$ for the Galactic extinction case. As a result, at $`𝒟5\times 10^4`$, $`C_{\mathrm{IR}}3.4\times 10^{10}M_{}\mathrm{yr}^1L_{}^1`$ for the SMC extinction, while $`C_{\mathrm{IR}}4.7\times 10^{10}M_{}\mathrm{yr}^1L_{}^1`$ for the Galactic extinction. Thus, we should be aware of the effect of extinction law on our conversion formula, but the effect is small (a factor of 1.4). If we remember that the uncertainty owing to the IMF is a factor of 2 (e.g., IHK00), we find that the dependence of $`C_{\mathrm{IR}}`$ on the the extinction curve is not very important. We note that we can estimate $`C_{\mathrm{IR}}`$ for the SMC-type extinction by moving the line in Fig. 3 to the left by 0.3 dex (a factor of 2).

### 2.6 Observational check

Here, we examine whether the values of $`f`$ and $`ϵ`$ that we adopted for the Galactic value (i.e., $`f=0.26`$ and $`ϵ=0.9`$) are consistent with the properties of actual H ii regions, by using the ratio between the ionizing-photon luminosity and the IR luminosity. This ratio was extensively studied in 1970s (e.g., Harper & Low harper71 (1971)).

We start with the relation by Petrosian et al. (petrosian72 (1972)) as

$`L_{\mathrm{IR}}=L(\mathrm{Ly}\alpha )+(1f)h\nu _{\mathrm{ion}}N_\mathrm{u}+ϵL_{\mathrm{nonion}},`$ (15)
where $`h\nu _{\mathrm{ion}}`$ is the mean energy of an ionizing photon, and $`L_{\mathrm{IR}}`$, $`L(\mathrm{Ly}\alpha )`$, and $`L_{\mathrm{nonion}}`$ are luminosities of IR dust emission, Ly$`\alpha `$, and nonionizing photons, respectively. We note that the meanings of $`f`$ and $`ϵ`$ in Eq. (15) are the same as those in Eq. (3). IHK00 derived $`L_{\mathrm{nonion}}=1.5h\nu _{\mathrm{ion}}N_\mathrm{u}`$ by adopting the Salpeter’s initial mass function (Salpeter salpeter55 (1955)) with the stellar mass range of 0.1–100 $`M_{}`$ and the mass–luminosity relation of stars (Table 3.13 of Binney & Merrifield binney98 (1998)). Then, according to Mezger (mezger78 (1978)), we define the IR excess (IRE) as

$`\mathrm{IRE}{\displaystyle \frac{L_{\mathrm{IR}}}{h\nu _{\mathrm{Ly}\alpha }fN_\mathrm{u}}}=a+(1f+1.5ϵ){\displaystyle \frac{h\nu _{\mathrm{ion}}}{h\nu _{\mathrm{Ly}\alpha }f}},`$ (16)
where $`\nu _{\mathrm{Ly}\alpha }`$ is the frequency of the Ly$`\alpha `$ radiation, and $`a`$ is the fraction of Ly$`\alpha `$ photons that reach the 2p state and go down to 1s (i.e., $`L(\mathrm{Ly}\alpha )=ah\nu _{\mathrm{Ly}\alpha }fN_\mathrm{u}`$)222The other photons decay to the ground state with a two-photon process.. According to S78, $`a0.67`$. Adopting $`f=0.26`$ and $`ϵ=0.9`$ for the Galactic values (Sects. 2.1 and 2.2), we obtain $`\mathrm{IRE}=11`$, where we have assumed that $`h\nu _{\mathrm{ion}}`$ is equal to $`h\nu `$ at the Lyman limit (912 Å). According to Harper & Low (harper71 (1971)), $`5\mathrm{IRE}10`$ (see also Aannestad aannestad78 (1978); Mezger mezger78 (1978); Maihara et al. maihara81 (1981)) for Galactic H ii regions. Fig. 7.3 in S78 showed that the IR luminosity of H ii regions is larger than predicted from Ly$`\alpha `$ luminosity by an order of magnitude (i.e., $`\mathrm{IRE}10`$). Thus, $`f`$ and $`ϵ`$ adopted in this paper are almost consistent with the properties of the Galactic H ii regions within the scatter of observed IRE.

Another test of our model is possible for the case of H ii regions in the Large Magellanic Cloud (LMC) by using the observational results of DeGioia-Eastwood (degioia92 (1992)). From Table 3 of DeGioia-Eastwood (1992), we obtain systematically lower IRE for the LMC H ii regions; IRE ranges from 1 to 4, and the median is 1.6. This can be naturally explained by our model with a lower value of the dust-to-gas ratio in the LMC. If we adopt $`𝒟=1.2\times 10^3`$ for the LMC (one-fifth of that of the galaxy) according to Issa et al. (1990), we obtain $`f0.7`$ and $`ϵ0.4`$. Then we obtain $`\mathrm{IRE}2.4`$ for the LMC. This value of IRE for the LMC lies in the range of the observed IRE. Thus, we consider Eqs. (11) and (13) to be applicable and useful in spite of simplification in their derivation. Therefore, we discuss the metallicity dependence of $`C_{\mathrm{IR}}`$ and its effect on the cosmic star formation history by using the relation between $`C_{\mathrm{IR}}`$ and $`𝒟`$ as shown in Fig. 3 in the following sections.

## 3 Metallicity dependence

In order to obtain the metallicity dependence of the conversion formula (Eq. (1)), which depends on the dust-to-gas ratio, $`𝒟`$, as shown in the previous section, we must relate $`𝒟`$ and the metallicity. Indeed, $`𝒟`$ and metallicity are positively correlated (Issa et al. issa90 (1990); Schmidt & Boller schmidt93 (1993); Lisenfeld & Ferrara lisenfeld98 (1998)). Here, we adopt the latest relation proposed by Hirashita (1999a, b). He has constructed a new evolution model of dust in the galactic environment, emphasizing the importance that the grains can grow via the accretion of the metal elements in the cool and neutral components of interstellar medium (ISM; see also Dwek dwek98 (1998)). Then, his relation between $`𝒟`$ and metallicity has explained the observational relation for both the giant and dwarf galaxies. Here, we adopt the relation as shown in Fig. 4 (solid line), since this seems to reproduce the relation for dwarf galaxies and spiral galaxies333As for the parameters in Hirashita (1999a, b), we choose the values as $`f_{\mathrm{in},\mathrm{O}}=0.1`$ and $`\beta _{\mathrm{acc}}=2\beta _{\mathrm{SN}}=10`$, where $`f_{\mathrm{in},\mathrm{O}}`$ is the dust mass fraction in the material injected from stars, and $`\beta _{\mathrm{acc}}`$ and $`\beta _{\mathrm{SN}}`$ are defined as gas consumption timescale (gas mass divided by star formation rate) normalized by dust growth timescale in clouds and by dust destruction timescale in supernova shocks, respectively (see Hirashita 1999a for details).. The metallicity is represented by \[O/H\]. We note that $`[\mathrm{O}/\mathrm{H}]=x`$ means that the abundance of oxygen is $`10^x`$ times the solar value (the solar oxygen abundance is assumed to be 0.013 in mass; e.g., Whittet whittet92 (1992)). In Fig. 4, observational data are also shown: The observed relations for nearby spiral galaxies are presented with the square (Issa et al. issa90 (1990)), and the area marked with dwarfs represents a typical locus of dwarf irregular galaxies and blue compact dwarf galaxies (Lisenfeld & Ferrara lisenfeld98 (1998)). We hereafter adopt the solid line in Fig. 4 as the relation between the dust-to-gas ratio and the metallicity.

Once we accept the relation in Fig. 4, we relate $`C_{\mathrm{IR}}`$ and \[O/H\] by using the relation between the dust-to-gas ratio and $`C_{\mathrm{IR}}`$ (Fig. 3). The relation between $`C_{\mathrm{IR}}`$ and \[O/H\] is presented in Fig. 5. We see that if \[O/H\] in the ISM evolves from $`2`$ to 0 via chemical evolution, the coefficient of the conversion from IR light to SFR becomes about 4–5 times smaller, again. Thus, when we would like to determine the SFR precisely within a factor of 4–5, we should not neglect the effect of metallicity.

## 4 Comments on the cosmic star formation history

The cosmological evolution of SFR is derived from the comoving density of galactic light. For example, Madau et al. (madau96 (1996)) applied the conversion formula from UV light to the SFR and showed that the SFR as a function of the redshift $`z`$ seems to have a peak at $`z1`$–2. The cosmic SFR as a function of $`z`$ has been revised and discussed by many authors (e.g., Steidel et al. steidel99 (1999)), commonly referred to as the Madau plot. Takeuchi et al. (t2\_00 (2000)) applied the formula by IHK00 to the determination of the cosmic star formation history from the number-count data of IR galaxies.

The cosmic SFR is also determined from the dust emission at the sub-millimeter (sub-mm) observational wavelength (Hughes et al. hughes98 (1998); Barger et al. barger98 (1998)). We can convert the sub-mm luminosity to the SFR by applying the formula proposed in IHK00. If we take into account the chemical evolution of galaxies on a cosmological timescale (e.g., Pei & Fall pei95 (1995)), however, we must examine the dependence of the conversion formula on metallicity as seen in the previous sections.

According to such a motivation, we apply the results obtained in the previous sections, especially Fig. 5 and its related discussions, to the cosmic star formation history. We consider the evolution of our conversion law along the metal enrichment with the aid of previous research that determined the metallicity as a function of $`z`$. We adopt Pei et al. (pei99 (1999)) as a recent modeling of the cosmic chemical evolution.

Pei et al. (pei99 (1999)) modeled the cosmic star formation history and chemical evolution and calculated the evolution of dust amount. Considering the absorption and reprocessing of light by dust, they determined the cosmic chemical evolution after the calibration with the observed comoving density of light. The resultant metallicity evolution as a function of $`z`$ is shown in the second column of Table 1. Here we adopt the solid line of Fig. 8 in Pei et al. (pei99 (1999)). Though their treatment of the evolution of the dust-to-gas ratio is not just the same as our treatment, for readers’ qualitative understanding, we present a relation between the metallicity and the redshift $`z`$ from their Fig. 8.

We also present $`C_{\mathrm{IR}}`$ in the third column of Table 1 calculated from the relation between metallicity ($`Z`$) and $`C_{\mathrm{IR}}`$ as shown in Fig. 5. Here we assume that $`\mathrm{log}(Z/Z_{})=[\mathrm{O}/\mathrm{H}]`$, where $`Z_{}`$ is the solar metallicity and $`Z`$ is measured in the galactic scale. From Table 1, we see that $`C_{\mathrm{IR}}`$ is 4–5 times larger at $`z3`$ than that at $`z0`$. Thus, we should carefully consider the chemical evolution of the galaxies if we determine the cosmic star formation history from the dust emission within a factor of 4–5. Deriving the cosmic star formation history from IR–sub-mm observation will be made possible by the future observational projects (e.g., Takeuchi et al. t2\_99 (1999)). Our formulation, which takes into account the metallicity dependence, will be useful in determining the cosmic star formation history.

We should mention that Pei et al. (pei99 (1999)) treated averaged quantities for each redshift and focused only on the redshift dependence. Thus, our Table 1 must be applied to the data averaged for each $`z`$. In other words, we should not apply them to each individual galaxy at a certain $`z`$. For each galaxy, Fig. 5 should be used instead after their metallicity is known. Finally, we should keep in mind the possibility that far-IR (FIR) or sub-mm sample of galaxies is biased to metal-enriched system because metal-poor (i.e., dust-poor) galaxies are not easily observed in these wavelengths.

## 5 Conclusion

Based on IHK00’s formulation, we consider the factor $`C_{\mathrm{IR}}`$ in the conversion formula between the IR luminosity and the SFR (Eqs. (2) and (3)). The factor $`C_{\mathrm{IR}}`$ becomes $`1.7\times 10^{10}`$, $`5.6\times 10^{10}`$, and $`7.9\times 10^{10}`$ $`M_{}\mathrm{yr}^1L_{}^1`$ for the metallicity of 1, 0.1, 0.01 times the solar value, respectively. Thus, $`C_{\mathrm{IR}}`$ differs by a factor of 4–5 in the range. Importantly, applicability of our formula is observationally confirmed by IREs of the Galaxy and the LMC. Applying our result to the cosmic star formation history, we have found that $`C_{\mathrm{IR}}`$ may be about 4–5 times larger at $`z3`$ than that at $`z0`$. This factor of 4–5 is larger than the uncertainty caused by the initial mass function of stars (a factor of 2). Thus, we should carefully consider the chemical evolution of galaxies if we determine the cosmic star formation history from the dust emission within a factor of 4–5.

###### Acknowledgements.

We first thank U. Lisenfeld, the referee, for useful comment that much improved the quality of this paper. We acknowledge T. T. Takeuchi and T. Totani for helpful discussions and S. Mineshige for continuous encouragement. One of us (HH) acknowledges the Research Fellowship of the Japan Society for the Promotion of Science for Young Scientists. We made extensive use of the NASA’s Astrophysics Data System Abstract Service (ADS).

## Appendix A dependence of cirrus fraction on metallicity

In the main body of this paper, we have assumed that the cirrus fraction of the IR luminosity, $`\eta `$, is 0.5. Here, we examine how $`C_{\mathrm{IR}}`$ is affected owing to the change of $`\eta `$. There is a large uncertainty about observational estimate of $`\eta `$. But, theoretically, we consider some simple cases where the change in $`\eta `$ cannot be neglected in our conversion formula.

The assumption of $`\eta =0.5`$ breaks if we are interested in starburst galaxies. When we examine a sample of starburst galaxies, indeed, it is reasonable to assume $`\eta 0`$. This is because young stars dominate the radiation field that heats the dust, and the optical depth of dust is so large that almost all of the bolometric luminosity is emitted in the IR (Soifer et al. soifer87 (1987); Kennicutt 1998b ). For the starburst galaxies, thus, the dependence of $`C_{\mathrm{IR}}`$ on the metallicity is obtained by putting $`\eta =0`$ into Eq. (3), while the dependence of $`f`$ and $`ϵ`$ on the metallicity is the same as described in Sect. 3. The result is shown in Fig. 6. We see that $`C_{\mathrm{IR}}`$ in this figure is larger by a factor of 2 than that in Fig. 5.

Next, we consider a possible change in $`\eta `$ of spiral galaxies as a function of metallicity. The cirrus fraction of spiral galaxies is not negligible at the present epoch ($`\eta 0.5`$ for nearby spiral galaxies; Lonsdale Persson & Helou lonsdale87 (1987)). In the beginning of a spiral galaxy’s evolution, when only a little metal is produced, $`\eta =0.0`$ may be reasonable. Thus, the first approximation for the dependence of $`\eta `$ on the dust-to-gas ratio is

$`\eta =0.5\left({\displaystyle \frac{𝒟}{6\times 10^3}}\right).`$ (17)
If we adopt this relation, we obtain Fig. 7 as the relation between $`C_{\mathrm{IR}}`$ and metallicity. We note that our simple treatment as Eq. (17) is only applicable for $`[\mathrm{O}/\mathrm{H}]0.1`$, because $`\eta >1`$ does not make sense. This means that we need a nonlinear modeling of the relation between $`\eta `$ and $`𝒟`$ at the present epoch. We may need to consider a complex mode of evolution of ISM, whose nonlinearity causes intermittent star formation history (Kamaya & Takeuchi kamaya97 (1997); Hirashita & Kamaya hirashitakamaya00 (2000); Takeuchi & Hirashita t2hirashita00 (2000)). Comparing $`C_{\mathrm{IR}}`$ at $`[\mathrm{O}/\mathrm{H}]=2`$ with that at $`[\mathrm{O}/\mathrm{H}]=0`$, we see that $`C_{\mathrm{IR}}`$ changes by an order of magnitude as the chemical enrichment proceeds. Thus, if we want to know the realistic cosmic evolution of galaxies from $`L_{\mathrm{IR}}`$, the metallicity evolution of cirrus must also be examined. We will try this very difficult theme in the near future.


text:
# From Anisotropy to Omega

## 1 Introduction

In the most general and simple inflationary scenario  the overall present energy density of the universe must be equal to the so-called critical energy density ($`\mathrm{\Omega }=\rho /\rho _c=1`$, with $`\rho _c=3H_0^2/8\pi G`$). In fact, the condition for inflation ($`\ddot{a}_{inf}>0`$, and so $`d/dt(H_{inf}^1/a_{inf})<0`$) is precisely that which drives $`\mathrm{\Omega }(t)`$ towards $`1`$ in the Friedmann equation

$$\mathrm{\Omega }(t)1=\frac{k}{a^2(t)H^2(t)}$$
(1)
during the inflationary period. This prediction, however, taken with the standard CDM model of structure formation, is apparently in disagreement with a combined set of observations, such as density-velocity galaxy field comparisons (), pairwise galaxy velocities (), X-ray clusters temperature function evolution (, ) and velocity dispersions 111A more conservative approach would say that the situation is rather unclear and the X-ray clusters data and theoretical modelling can be sufficiently large to prevent an unambigous exclusion of $`\mathrm{\Omega }_m=1`$ (see e.g. , , , )()\], which point in favor of a low density universe 222From here on, with $`\mathrm{\Omega }`$ we will indicate only the present value ($`\mathrm{\Omega }_{matter}<1`$). In order to solve the discrepancy, open ($`\mathrm{\Omega }=\mathrm{\Omega }_{matter}<1`$) inflationary models have been proposed (,  ). This type of model highlights limitations on the predictiveness of the inflationary scenario, which is supposed to have the advantage of removing any dependence on initial conditions from our present observable universe. However, even if this picture leads to a more complicated phenomenology, it generally determines $`\mathrm{\Omega }`$ directly from parameters of the physical theory.

On the other hand, another way to keep low-density models compatible with the simplest model of inflation is to introduce a cosmological constant, $`\mathrm{\Lambda }`$, such that $`\mathrm{\Omega }\mathrm{\Omega }_{matter}+\mathrm{\Omega }_\mathrm{\Lambda }=1`$. The presence of such a cosmological constant, which is compatible with all the above observations since they are practically insensitive to it, is also preferred by measurements of the magnitude-redshift relationship in high-redshift type Ia supernovae (, ). Nonetheless, this ’natural’ solution introduces the cosmological constant problem () that is perhaps even more acute in inflationary cosmology (see  for a review). Thus, an accurate determination of the present overall density parameter $`\mathrm{\Omega }`$, even if not a panacea for the cosmological scenario, is at least extremely important in understanding which theoretical framework could explain the above conflicts.

There is much experimental evidence for the presence of a peak in the CMB angular power spectrum (, ). Furthemore, with the recent release of the BOOMERanG-98  and MAXIMA-1  spectra, the shape and position of the peak has been detected with unprecedented accuracy. This result, apart from being in wonderful agreement with the standard scenario of primordial adiabatic fluctuations, has important consequences on the parameter $`\mathrm{\Omega }`$ (, , , , ). As previously noted and already well explained in the literature, the effect of the curvature is to change the relationship between the physical scales on the Last Scattering Surface (LSS) and the corresponding angular scales (see., e.g. , ). In an open universe, for example, the geodetics focalize in such a way that a particular angular scale will correspond to a greater physical scale on the LSS than the one expected in a flat model. The immediate result is a shift in the radiation angular power spectrum (the so-called $`C_{\mathrm{}}`$’s), and thus a dependence of the position of the first peak $`\mathrm{}_{peak}`$ on the curvature (, , , , ). The CMB power spectrum is therefore a powerful tool for the determination of the curvature and so of the overall energy density.

## 2 $`\mathrm{\Omega }`$ and the shift of the CMB Angular Power Spectrum

In Fig.1, we show the recent data from the BOOMERanG (,) and from the MAXIMA  experiments together with the predictions of the most viable open, flat and closed adiabatic models (from present non–CMB observations). It is quite evident that the open (closed) model predicts a first peak on smaller (larger) scales with respect to the flat model, which is in much better agreement with the data.

In the literature, the dependence of $`\mathrm{}_{peak}`$ on $`\mathrm{\Omega }`$ is often expressed as

$$\mathrm{}_{peak}\frac{220}{\sqrt{\mathrm{\Omega }}}$$
(2)
However, this approximation is not correct in $`\mathrm{\Lambda }`$–dominated universes (see Fig.1, where the peak for the closed model is at $`130`$ instead of $`190`$) and does not take into account the further dependence of $`\mathrm{}_{peak}`$ on other parameters like the Hubble constant and the matter density $`\mathrm{\Omega }_m`$. So, in view of the recent $`\mathrm{}_{peak}=197\pm 6`$ (), further modifications to the above formula are needed (see also , ). Also, the width of the peaks and the inter peak distance vary as functions of $`\mathrm{\Omega }`$ (again see Fig.1) so a more complete expression is needed to describe these effects. This point is rather important because knowing the exact dependence of the CMB spectrum on $`\mathrm{\Omega }`$ will help us in understanding the shape of the probability distribution function for this parameter and, ultimately, how well it can be measured independently.

In order to do this, it is convenient to introduce the ’shift’ parameter $`R`$ defined as

$$\mathrm{}_{peak}=\mathrm{}_{peak}^{flat}/R$$
(3)
where $`flat`$ indicates the flat, pure–CDM, $`\mathrm{\Lambda }=0`$ model. The use of the $`R`$ parameter is more appropriate than the conventional $`\mathrm{}_{peak}`$ because it has a clearer geometrical dependence.

As usual, let us assume the metric of space–time to be of a Friedman–Lemaitre–Robertson–Walker (FLRW) form with curvature $`k`$:

$$ds^2=a(t)^2[d\eta ^2+\gamma _{ij}dx^idx^j]$$
(4)
with

$$\gamma _{ij}dx^idx^j=dr^2+\chi ^2(r)(d\theta ^2+sin^2\theta \varphi ^2)$$
(5)
where the function $`\chi (r)`$ depends on the curvature $`k`$ and is $`r`$, $`\mathrm{sin}(r)`$ or $`\mathrm{sinh}(r)`$ for $`k=0`$ (flat models), $`k=1`$ (closed models) and $`k=1`$ (open models), respectively. The position of the first acoustic peak is determined primarly by the angle subtended by the acoustic horizon $`\lambda _{ac}`$ at decoupling time, $`\eta _{dec}`$. The angle under which a given comoving scale $`\lambda `$ at conformal time $`\eta _{dec}`$ is seen on the sky is given by $`\theta (\lambda )=\lambda /\chi (\eta _0\eta _{dec})`$. As the harmonic number $`\mathrm{}`$ is inversely proportional to the angle $`\theta `$, we have $`R=\theta _{ac}/\theta _{ac}^{flat}`$, with $`\theta _{ac}=c_s\eta _{dec}/\chi (y)`$, where $`c_s=1/\sqrt{3(1+3\mathrm{\Omega }_b/4\mathrm{\Omega }_{rad}(1+z_{dec}))}`$ denotes the adiabatic sound speed of the baryon/photon plasma at decoupling. It is possible to show that :

$$\eta _{dec}=\frac{2\sqrt{|\mathrm{\Omega }_k|}}{\mathrm{\Omega }_m}(\sqrt{\mathrm{\Omega }_{rad}+\mathrm{\Omega }_m/(z_{dec}+1)}\sqrt{\mathrm{\Omega }_{rad}})$$
(6)
where $`z_{dec}1100`$ is the redshift at decoupling and $`\mathrm{\Omega }_k=1\mathrm{\Omega }_m\mathrm{\Omega }_\mathrm{\Lambda }`$. Furthermore, we have:

$$y=\eta _0\eta _{dec}=\sqrt{|\mathrm{\Omega }_k|}_0^{z_{dec}}\frac{dz}{[\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_k(1+z)^2+\mathrm{\Omega }_\mathrm{\Lambda }]^{1/2}}$$
(7)
For a flat, $`\mathrm{\Omega }_m=1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$, universe we have

$$\theta _{ac}^{flat}=c_s\eta _{dec}/(\eta _0\eta _{dec})=c_s(\sqrt{\mathrm{\Omega }_{rad}+1/(z_{dec}+1)}\sqrt{\mathrm{\Omega }_{rad}})$$
(8)
We then find, keeping constant $`\mathrm{\Omega }_bh^2`$ (and so $`c_s`$) and $`\mathrm{\Omega }_mh^2`$ (see next section) :

$$R=\frac{2}{\chi (y)}\sqrt{\frac{|\mathrm{\Omega }_k|}{\mathrm{\Omega }_m}}$$
(9)
which is a quantity just dependent on $`\mathrm{\Omega }_k=1\mathrm{\Omega }_\mathrm{\Lambda }\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_m`$.

In Fig.2 we plot the contours at $`R=const`$ in the $`\mathrm{\Omega }_\mathrm{\Lambda }\mathrm{\Omega }_m`$ plane, together with one dashed line such that $`\mathrm{\Omega }=1`$. As we can see, open models have in general $`R<1`$, while closed models have $`R>1`$. This is the usual result that in an open model the peaks are shifted towards greater $`\mathrm{}`$ values (smaller angular scales) with respect to the flat model case, while for closed models we have the opposite effect. It is worth noting that the CMB angular shift is not linearly related to the curvature: lines at $`\mathrm{\Omega }=const`$ are not parallel with the $`R=const`$ contours, but have multiple intersections, especially in regions far from flatness (again, see Fig.2).

In Section 4, we will build a likelihood distribution function $`L`$ for $`\mathrm{\Omega }`$, using a Bayesian approach to compare the current CMB data with the theoretical predictions and then marginalizing over the remaining cosmological parameters. The above result implies that the probability distribution function for $`\mathrm{\Omega }`$ (if we assume a flat prior distribution on the remaining parameters while marginalizing) will always be ’skewed’ so that $`\mathrm{\Omega }`$ will never be measured at a level better than $`1020\%`$333Another interesting point, is that for $`\mathrm{\Omega }_m0`$ the lines at $`R`$=constant converge towards $`\mathrm{\Omega }_\mathrm{\Lambda }=1`$, but we must warn the reader that in the de Sitter solution the notion of open, flat or closed Universe becomes ambigous ..

## 3 The geometrical degeneracy

With the parameter $`R`$ fixed, the structure and position of the $`C_{\mathrm{}}`$ spectrum is dependent on $`2`$ physical scales: the equality scale and the sound horizon at decoupling scale. These quantities are completely defined once we choose the abundance of cold dark and baryonic matter in our model, by the parameters: $`\omega _{cdm}=\mathrm{\Omega }_{cdm}h^2`$ and $`\omega _b=\mathrm{\Omega }_bh^2`$. The CMB spectrum also depends on the characteristics of the primordial inflationary perturbations. Assuming that we have already selected the primordial power spectrum of our model, both in shape (tilted or ’blue’) and in nature (adiabatic, isocurvature, hybrid), the structure of the CMB angular peaks is completely determined by $`R`$, $`\omega _b`$ and $`\omega _m=\omega _b+\omega _{cdm}`$. cdm This result has an important consequence: if we let $`\mathrm{\Omega }_b`$, $`\mathrm{\Omega }_{cdm}`$ and $`h`$ assume any value but a fixed $`\omega _m`$ and $`\omega _b`$, the lines at $`R=const`$ in the $`\mathrm{\Omega }_m\mathrm{\Omega }_\mathrm{\Lambda }`$ plane correspond to sets of degenerate power spectra with an identical shape on subdegree angular scales ($`\mathrm{}>30`$) .

In Fig.3 we draw a set of degenerate models with $`R`$, $`\omega _m`$ and $`\omega _b`$ set to those of the cosmological concordance model, $`\mathrm{\Omega }_m=0.35`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.65`$, $`h=0.65`$, $`\omega _b=0.021`$, $`n_s=1`$. It is clear from the degeneracy of the models that, once $`\mathrm{\Lambda }`$ is included, the peak position is not directly related to $`\mathrm{\Omega }`$. It is also clear that it seems impossible to obtain any relevant and independent information from the CMB on $`\mathrm{\Omega }_\mathrm{\Lambda }`$: the Integrated Sachs-Wolfe effect on large scales could in principle break the degeneracy but cosmic variance and a possible presence of a gravity wave background make this effect difficult to disentangle.

Our main results are then the following:

* Lines at $`R=const`$ in the $`\mathrm{\Omega }_\mathrm{\Lambda }\mathrm{\Omega }_m`$ plane correspond to sets of degenerate $`C_{\mathrm{}}`$ power spectra.
* Given a flat $`\mathrm{\Omega }_\mathrm{\Lambda }0.65`$ model, a degenerate closed model can be found by decreasing $`\mathrm{\Omega }_\mathrm{\Lambda }`$ and $`h`$ and increasing $`\mathrm{\Omega }_b`$ and $`\mathrm{\Omega }_{cdm}`$.
* $`R`$, $`\omega _m`$ and $`\omega _b`$ are the most meaningful CMB anisotropy observables.
* The CMB spectrum is a useful tool for the determination of $`\mathrm{\Omega }`$ only if we live in a flat universe ($`R=const`$ are not parallel to $`\mathrm{\Omega }=const`$ when $`\mathrm{\Omega }`$ greatly differs from one).
* Assuming that the concordance model describes our real universe we expect that a likelihood analysis for $`\mathrm{\Omega }`$, using only the CMB power spectrum and without including any external information about $`h`$, $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$, will always be skewed towards closed models.

## 4 From CMB to $`\mathrm{\Omega }`$

Let us now describe the standard tools for extracting the cosmological parameters from CMB anisotropy observations. Here we will analyze the recent BOOMERanG97 (), BOOMERanG98 () and MAXIMA () results. The power spectra from these experiments were estimated in $`7`$, $`12`$ and $`10`$ bins respectively, spanning the range $`25\mathrm{}785`$. In each bin, the spectrum is assigned a flat shape, $`\mathrm{}(\mathrm{}+1)C_{\mathrm{}}/2\pi =C_B`$. Following  we use the offset lognormal approximation to the likelihood $`L`$. In particular we define:

$$2\mathrm{l}\mathrm{n}L=(D_B^{th}D_B^{ex})M_{BB^{}}(D_B^{}^{th}D_B^{}^{ex}),$$
(10)
$$D_B^X=\mathrm{ln}(C_B^X+x_B),$$
(11)
$$M_{BB^{}}=(C_B^{ex}+x_B)F_{BB^{}}(C_B^{}^{ex}+x_B^{}),$$
(12)
where $`C_B^{th}`$ ($`C_B^{ex}`$) is the theoretical (experimental) band power, $`x_B`$ is the offset correction and $`F_{BB^{}}`$ is the Gaussian curvature of the likelihood matrix at the peak. Of course, $`C_B^{th}`$ will depend on the various parameters of our cosmological model ($`\omega _m`$,$`n_s`$, …), and so it will be the likelihood function $`L`$. In order to compute the likelihood for a given parameter only $`\alpha `$ we can either marginalize over all the remaining parameters, namely carry out the integral

$$L_{marg}(\alpha )=P_{prior}(\alpha ,\stackrel{}{\mathrm{\Pi }})L(\alpha ,\stackrel{}{\mathrm{\Pi }})𝑑\stackrel{}{\mathrm{\Pi }}$$
(13)
where $`\stackrel{}{\mathrm{\Pi }}`$ is a vector containing all the remaining parameters and $`P_{prior}`$ is the prior probability distribution, or we can maximize i.e. for a fixed $`\alpha `$ find the $`\stackrel{}{\mathrm{\Pi }}_{max}`$ wich maximizes

$$L_{max}(\alpha )=P_{prior}(\alpha ,\stackrel{}{\mathrm{\Pi }}_{max})L(\alpha ,\stackrel{}{\mathrm{\Pi }}_{max}).$$
(14)
The two methods, in general, agree at a level of $`10\%`$. The maximization method is usually based on a search algorithm through the second derivative of the likelihood matrix (). In this approach the $`C_{\mathrm{}}`$ spectrum is computed on the way, without sampling the whole parameter space. A different approach is based on building a database of $`C_{\mathrm{}}`$’s on a discretized grid of the parameter space. $`L(\mathrm{\Omega })`$ is then obtained by maximizing and/or integrating the likelihood computed on the grid (, ). Of course, producing a grid of models can be quite computationally expensive even with the new and fast boltzmann codes like CAMB  or CMBFAST . But this problem can be drastically reduced using morphing  or interpolation  algorithms.

The definition of the database is rather important because it defines the internal $`P_{prior}`$ of our analysis and, in general, it is better to have this prior as flat as possible for each given parameter. This brings us to the choice of the variables in which the database must be sampled. As we saw in the previous section, the CMB anisotropies are mainly sensitive to the physical variables $`\omega _i`$ and $`R`$, so sampling in those variables will avoid degeneracies. Furthemore, the physical baryon density $`\omega _b=\mathrm{\Omega }_bh^2`$ is well determined by independent measurements like primordial nuclide abundances, so this is the optimal choice for extracting information about this parameter (without involving complicated Jacobian transformations) or assume external priors for it. For the same reason, extracting confidence limits on parameters like $`h`$ can be a little more elaborate with this sampling, being the database in $`h^2`$.

Another possibility is to sample the database in cosmological variables, like $`\mathrm{\Omega }_i`$ and $`h`$. Of course, this will introduce degenerate models in our database but this sampling has the advantage of obtaining direct constraints on the commonly used parameters and with flat prior distributions. In most of the recent papers, either a ’hybrid’ variables approach, with sampling in $`\omega _i`$ and $`1\mathrm{\Omega }_m\mathrm{\Omega }_\mathrm{\Lambda }`$, or the cosmological variables approach has been used. Here we will choose the database approach, sampling the parameter space in cosmological  variables as follows: $`\mathrm{\Omega }_m=0.1,\mathrm{},1.1`$; $`\mathrm{\Omega }_\mathrm{\Lambda }=0.0,\mathrm{},1.0`$; $`\mathrm{\Omega }_b=0.01,\mathrm{},0.25`$; $`h=0.25,\mathrm{},0.95`$ and $`n_s=0.50,\mathrm{},1.50`$. We will not consider the possibility of high redshift reionization of the intergalactic medium $`\tau _c>>0`$, a gravity waves contribution or the effect of massive neutrinos.

## 5 Removing the geometrical degeneracy: Results

In Fig.4 we plot our likelihood contours for $`\mathrm{\Omega }`$, $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ using just the intrinsic internal priors of the database plus the quite reasonable age prior $`t_{universe}>10Gyr`$. As we can see, the likelihood for $`\mathrm{\Omega }`$ is skewed towards closed models but is consistent with flatness. It is rather important to note that this skewness is largely due to the $`R`$-degeneracy which makes ’more’ closed models compatible with the observations. The likelihood for $`\mathrm{\Omega }_\mathrm{\Lambda }`$ and $`\mathrm{\Omega }_m`$ are quite flat due to the geometrical degeneracy but they nonetheless feel border effects from the database priors. The likelihood for $`\mathrm{\Omega }`$ starts to be in even more agreement with flatness when a Gaussian prior $`h=0.65\pm 0.2`$ is assumed as in Fig.5. This clearly shows that most of the degeneracies in the $`h<0.5`$ region are well removed by the prior. Including a prior $`\mathrm{\Omega }_m=0.3\pm 0.2`$ (Fig.6) as suggested by the majority of measurements, shrinks the likelihood towards $`\mathrm{\Omega }=1`$ and gives a strong determination for the cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }>0`$ at $`4\sigma `$ level. The complementarity with the supernovae type Ia measurements is even more clear in Fig.7, where a combined CMB+SnIA analysis gives $`\mathrm{\Omega }_\mathrm{\Lambda }>0`$ at more than $`7\sigma `$ and $`\mathrm{\Omega }=1`$ with a few percent uncertainty.

## 6 Conclusions

The BOOMERanG and MAXIMA data support the main prediction of the inflationary paradigm: that the geometry of the universe is flat. The small deviations towards closed models reported in various analyses (,) can be easily explained by the degeneracies in the cosmological parameters, which make more closed models compatible with the data. These conclusions are considerably strengthened by the inclusion of other cosmological data such as measurements of the Hubble constant, the overall matter density $`\mathrm{\Omega }_m`$ and the accelerating expansion rate indicated by observations of distant Supernovae. At the same time, $`\mathrm{\Omega }_\mathrm{\Lambda }`$ and other forms of ’dark energy’ cannot be well determined by the CMB data alone, in spite of their high precision. This does not mean that CMB measurement are not useful in the determination of such parameters: combining the CMB data with constraints from observations of large-scale-structure and from observations of SN-Ia increases the extent to which $`\mathrm{\Lambda }`$ can be quantified, with $`\mathrm{\Omega }_\mathrm{\Lambda }>0`$ at $`7\sigma `$. Future data from the PLANCK and/or SNAP satellites will hopefully enable us to resolve the ’dark energy’ puzzle.

## Aknowledgements

We would like to thank Monique Signore, Paolo de Bernardis, Marian Douspis, Ruth Durrer, Pedro Ferreira, Graca Rocha, Joe Silk, Filippo Vernizzi and Nicola Vittorio. We also acknowledge using CMBFAST, CAMB and Lloyd Knox’s RADPack  packages for the COBE data.


text:
# Counterion correlations and attraction between like-charged macromolecules

## I Introduction

In the last few years a new phenomenon has attracted attention of the community of soft condensed matter physicists — appearance of attraction between like charged macromolecules in solutions containing multivalent ions. The problem is particularly fascinating because it contradicts our well established intuition that like charged entities should repel . The fundamental point, however, is that the electrolyte solutions are intrinsically complex systems for which many body interactions play a fundamental role.

The attraction between like charged macromolecules is important for many biological systems. One particularly striking example is provided by the condensation of DNA by multivalent ions such as $`\text{Mn}^{2+}`$, $`\text{Cd}^{2+}`$ and various polyamines . This condensation provides an answer to the long standing puzzle of how a highly charged macromolecule, such as the DNA, can be confined to a small volume of viral head or nuclear zone in procaryotic cell. Evidently, the multivalent ions serve as a glue which keeps the otherwise repelling like-charged monomers in close proximity . In eukaryotic cells, the cytosol is traversed by a network of microtubules and microfilaments — rigid chains of highly charged protein (F-actin) — which in spite of large negative charge agglomerate to form filaments of cytoskeleton . The actin fibers are also an important part of the muscle tissue, providing a rail track for the motion of molecular motor myosin.

Although the nature of attraction between like charged macromolecules is still not fully understood, it seems clear that the attractive force is mediated by the multivalent counterions . A strong electrostatic attraction between the polyions and the oppositely charged multivalent counterions produces a sheath of counterions around each macromolecule. The condensed counterions can become highly correlated resulting in an overall attraction. It is important to note that the complex formed by a polyion and its associated counterions does not need to be neutral for the attraction to arise. Under some conditions the correlation induced attraction can overcome the monopolar repulsion coming from the net charge of the complexes.

Recently a simple model was presented to account for the attraction between two lines of charges . Each line had $`Z`$ discrete uniformly spaced monomers of charge $`q`$, and $`nZ/\alpha `$ condensed counterions of charge $`\alpha q`$ free to move along the rod. The net charge of such a polyion-counterion complex is $`Q_{\mathrm{eff}}=(Z\alpha n)q0`$. Nevertheless, it was found that if $`n>Z/2\alpha `$ and $`\alpha 2`$, at sufficiently short distances, the two like-charged rods would attract . It was argued that the attraction resulted from the correlations between the condensed counterions and reached maximum at zero temperature. If $`n<Z/2\alpha `$ the force was always found to be repulsive.

Clearly, a one dimensional line of charge is a dramatic oversimplification of the physical reality. If we are interested in studying the correlation induced forces between real macromolecules their finite radius must be taken into account . Thus, a much more realistic model of a polyion is a cylinder with a uniformly charged backbone  or with an intrinsic charge pattern  as, e.g., the helix structure of DNA molecule. Furthermore, the condensed counterions do not move along the line, but on the surface of the cylinder. Unfortunately, these extended models are much harder to study analytically.

In this paper we explore the effects of finite polyion diameter on the electrostatic interactions between the two polyions using Monte Carlo simulations. We find that the finite diameter and the associated angular degrees of freedom of condensed counterions significantly modify the nature of attraction. Thus, although there is still a minimum charge which must be neutralized by the counterions in order for the attraction to appear, this fraction is no longer equal to $`50\%`$ as was the case for the line of charge model. We find that the critical fraction depends on the valence of counterions and is less than $`50\%`$ for $`\alpha 2`$. For monovalent counterions no attraction is found. The crystalline structure of the condensed counterions, as first suggested by simulations of Gronbech-Jensen et al.  and Refs. , is also not very obvious. In particular we find very similar distributions of condensed counterions in the regime of attractive and repulsive interactions.

The structure of this paper is as follows. The model and the method of calculation are described in section II. In section III, we present the results of the simulations. The conclusions are summarized in section IV.

## II Model and method

### A The model

The DNA model considered here is an extension of the one proposed earlier by Arenzon, Stilck and Levin . A similar model has been recently discussed by Solis and Olvera de la Cruz . The polyions are treated as parallel rigid cylinders of radius $`R`$ and $`Z`$ ionized groups, each of charge $`q`$, uniformly spaced — with separation $`b`$ — along the principle axis, Fig. 1. Besides the fixed monomers, each polyion has $`nZ/\alpha `$ condensed counterions with valence $`\alpha `$ and charge $`\alpha q`$, which are constrained to move on the surface of the cylinder. To locate a condensed counterion it is necessary to provide its longitudinal position, $`z`$ ($`0z<Z`$), and the transversal angle, $`\theta `$ ($`0\theta 2\pi `$). To simplify the calculations, the angular and the longitudinal degrees of freedom are discretized, see Fig. 1. The surface of the cylinder is subdivided into $`Z`$ parallel rings with a charged monomer at the center of each ring. Each ring has $`m`$ sites available to the condensed counterions, see Figs. 1 and 2. The hardcore repulsion between the particles requires that a site is occupied by at most one condensed counterion. The two polyions are parallel, with the intermolecular space treated as a uniform medium of dielectric constant $`ϵ`$.

We introduce occupation variables $`n_i^k`$ for the two polyions, so that $`i=1,2,\mathrm{},Z(m+1)`$ and $`k=1,2`$. Thus, $`n_i^k=1`$ if the $`i`$’th site of the $`k`$’th polyion is occupied by a particle of valence $`\alpha _i=\{1,\alpha \}`$ (negative core charge or counterion of valence $`\alpha `$, respectively), otherwise $`n_i^k=0`$. Note that the core charge is always “occupied”, while the counterions are free to move between the $`Zm`$ ring sites of each polyion. The Hamiltonian for the interaction between the two polyions is,

$$\beta H=\xi \underset{k,l}{}\underset{i,j}{}\frac{\alpha _i^kn_i^k\alpha _j^ln_j^l}{|𝒓_i^k𝒓_j^l|},$$
(1)
with the $`ij`$. All lengths are measured in units of $`b`$ (for DNA, $`b=1.7`$Å). The dimensionless quantity, $`\xi =\beta q^2/bϵ`$, is the Manning parameter, which for DNA is $`\xi =4.17`$. The partition function is obtained by tracing over all the possible values of $`n_i^k`$ consistent with the constraint of fixed number of counterions per polyion,

$$𝒬=\underset{\{n_i^k\}}{}\mathrm{exp}(\beta H).$$
(2)
Clearly, this is a very crude model of the interaction between two macromolecules in a polyelectrolyte solution. The molecular nature of the solvent is ignored. Also the number of condensed counterions is fixed instead of being dependent on the separation between the particles. Nevertheless, we believe that this simple model can provide some useful insights for the mechanism of attraction in real polyelectrolyte solutions.

### B The observables

We are interested in statistical averages of observables such as the energy and the force between the two polyions. Furthermore, to understand the nature of the interaction between the two macromolecules it is essential to study the correlations between the condensed counterions on the two polyions.

The force is obtained from the partition function Eq. (2), $`b\beta 𝑭=\mathbf{}(\mathrm{ln}𝒬)`$. From symmetry, only the $`y`$-component is different from zero.

For finite macromolecules the symmetry between the two polyions can not be broken . Hence it is impossible to produce a true crystalline order in a finite system at non-zero temperature. Since within our simplified model the two polyions have exactly the same number of condensed counterions, the average angular counterion distribution $`n_i^k(z,\theta _i^k)`$ must be symmetric with respect to the mid-plane $`y=d/2`$, see Fig. 1. The angle $`\theta _i^k`$ labels the site $`i`$ on polyion $`k`$, see Fig. 2. Thus, $`n_3^2(z,\theta _3^2)`$ denotes the occupation variable for the site 3 on the ring $`z`$ located on polyion 2, with an angle of $`3(2\pi /m)`$. Indeed, Fig. 3 shows that the density profiles are completely symmetric (up to fluctuations). In spite of this symmetry it is possible for the counterions on the two polyions to become highly correlated. Clearly, the strength of these correlations will depend on the product $`\xi \alpha ^2`$ and the separation between the two macromolecules. Considering Fig. 2, it is evident that if the site $`2`$ on the first polyion is occupied, the likelihood of occupation of the site $`8`$ on the second polyion will be reduced.

To explore the nature of electrostatic correlations, we define a counterion-hole correlation function between the adjacent rings on the two polyions,

$`c_{zij}^{12}`$ $`=`$ $`n_i^1(z,\theta _i^1)\left[1n_j^2(z,\theta _j^2)\right]`$ (4)
$`n_i^1(z,\theta _i^1)\left[1n_j^2(z,\theta _j^2)\right].`$
Here $`\mathrm{}`$ denotes the ensemble average. This function should be non-zero when sites on the two polyions are correlated, that is if one is occupied by a condensed counterion there is an increased probability of the second being empty.

### C Simulations

To calculate the force between the two polyions, we have performed a standard Monte-Carlo (MC) simulation with the usual Metropolis algorithm . First, one counterion on polyion 1 is randomly chosen and displaced to a vacant position on the same polyion. This move is accepted or rejected according to the standard detailed balance criterion . We do not permit exchange of particles between the polyions. Next, the same is done for polyion 2. In one Monte-Carlo step (MCS) all $`2n`$ condensed counterions on the two polyions are permitted to attempt a move.

The long-ranged nature of the Coulomb interaction requires evaluation of all the pair interactions in Eq. (1) at every MCS. Due to the limited computational power available to us, we have confined our attention to relatively small systems with $`Z<100`$ and $`m=10`$. We have checked, however, that for $`m=10`$ the force has already reached the continuum limit and did not vary further with increase of $`m`$. Also we note that the “thermodynamic limit” is reached reasonably quickly, so that there is a good collapse of data already for $`Z>50`$, see Fig. 4. 2000 MCS served to equilibrate the system, after which 500 samples were used to calculate the basic observables, namely, the mean force and energy. To obtain the correlation functions, 5000 samples were used with 5000 MCS for equilibration.

## III Results and Discussion

The simulations were performed for $`\xi =4.17`$ and $`R/b=8.2`$, relevant for DNA. For monovalent counterions the simulation results indicate that the force is purely repulsive. This is in complete agreement with the experiments , which do not find any indication of DNA condensation for monovalent counterions.

For divalent counterions the force between the two complexes can becomes negative, indicating appearance of an effective attraction, Fig. 4. The range of attraction is larger than was found for the one dimensional line of charge model, Ref. .

Within the Manning theory  $`88\%`$ of the DNA’s charge is neutralized by the divalent counterions. However, there are indications that even a larger fraction of DNA’s charge can become neutralized by the multivalent ions if the counterion correlations are taken into account . In this case the interaction is purely attractive, with the range of about $`d/b20`$ or $`34\text{Å}`$ ($`7\text{Å}`$ surface-to-surface), Fig. 4.

A minimum number of condensed counterions is necessary for attraction to appear. In Fig. 5 we present the surface-to-surface separation, $`(d_02R)`$, below which the force between the two complexes becomes negative (attractive), as a function of the number of multivalent counterions. For the case of DNA with divalent counterions $`\alpha =2`$, the attraction appears only if $`40\%`$ of the core charge is neutralized. For $`\alpha =3`$ this fraction decrease to $`30\%`$. Furthermore, decrease in the value of the Manning parameter, $`\xi `$, increases the minimum number of condensed counterions necessary for the attraction to appear. This is fully consistent with the fact that the attraction is mediated by the correlations between the condensed counterions. Since raise in temperature tends to disorganize the system, the state of highest correlation between the condensed counterions corresponds to $`T=0`$ or $`\xi =\mathrm{}`$.

The surface-to-surface distance at which the attraction first appears tends to zero as the number of condensed counterions is diminished. We find $`d_02R(\mu \mu _c(\alpha ))^\nu `$, where the average counterion concentration is $`\mu =n/Z`$ and the critical fraction $`\mu _c`$ depends on the valence of condensed counterions $`\alpha `$. From Fig. 5 it is evident that $`\nu =1`$. This should be contrasted with the line of charge model Ref. , for which $`\nu =1/3`$.

In Fig. 6 we show two snapshots of the characteristic equilibrium configurations for (a) $`d/b=32.8`$ and (b) $`d/b=16.8`$. Looking at this figures it is difficult to see something that would distinguishes between them, both appear about the same. There is no obvious crystallization or transversal polarization suggested in previous studies . Yet, the case (a) corresponds to the repulsive, while the case (b) corresponds to the attractive interaction between the polyions. To further explore this point, in Figs. 7 and 8 we present the site-site correlation function, Eq. (4), for macromolecules with $`Z=20`$ and $`n=7`$. For $`d/b=32.8`$ the surface-to-surface distance between the two polyions is sufficiently large for their condensed counterions to be practically uncorrelated, Fig. 7. On the other hand, for $`d/b=16.65`$ strong correlations between the condensed counterions are evident, Fig. 8. The Fig. 8 shows that the sites two and three on the first polyion are strongly correlated with the sites seven and eight on the second polyion, respectively. It is these correlations between the adjacent sites on the two polyions which are responsible to the appearance of attraction between the two macromolecules when they are approximated, Fig. 4.

## IV Summary

We have presented a simple model for polyion-polyion attraction inside a polyelectrolyte solution. It is clear from our calculations that the attraction results from the correlations between the condensed counterions and reaches maximum for $`T=0`$. The thermal fluctuations tend to diminish the correlations, decreasing the amplitude of the attractive force. Consistent with the experimental evidence, the attraction exists only in the presence of multivalent counterions. Our simulations demonstrate that a critical number of condensed counterions is necessary for the appearance of attraction. The fraction of bare charge that must be neutralized for the attraction to arise depends on the valence of counterions. The larger the valence, the smaller the fraction of the bare polyion charger that must be neutralized for the attraction to appear. This result should be contrasted with the line of charge model  for which the critical fraction was found to be equal to $`50\%`$, independent of the counterion charge.

###### Acknowledgements.

We thank J. J. Arenzon for helpful comments on simulations. This work was supported by CNPq — Conselho Nacional de Desenvolvimento Científico e Tecnológico and FINEP — Financiadora de Estudos e Projetos, Brazil.