text:
# The Early Universe Odyssey With Gravitational Waves

## 1 Which sources of gravitational waves will be detected first ?

The comprehensive contributions of Coccia  and Giazotto  have described the outstanding experimental effort that is now going on at detecting gravitational waves. As we all know, great expectations are related with the coming in operation laser interferometers: the British-German GEO600 , the two American instruments LIGO , and the French-Italian VIRGO . It is important to know in advance, and to be prepared to monitor, those sources of gravitational waves that are likely to be detected first in the ongoing and forthcoming observations. Before proceeding to the main topic of my contribution - relic gravitational waves - I would like to present a theorist’s view on this issue. I will briefly report on the conclusions of a detailed study undertaken in Ref.. By comparing our theoretical analysis with only 3 out of several experimental programs, we, of course, do not mean to diminish the importance of work of other experimental groups.

It is common to divide possible sources of gravitational waves in three broad categories: quasi-periodic, explosive, and stochastic. We are aware of definitely existing astrophysical sources belonging to each category. To be interesting from the point of view of its detection, the source should be sufficiently powerful, should fall in the frequency band of the detector, and should occur reasonably often during the life-time of the instrument. After having analysed many possible sources, and taking all the factors into account, the authors of Ref. gave their preference to compact binaries consisting of black holes and neutron stars. It is clear that in order to radiate large-intensity gravitational waves at frequencies accesible to ground-based instruments, the objects forming a pair should be massive and should orbit each other at very small separations - a few hundreds kilometers. According to the existing views, such massive objects can only be the end-products of stellar evolution - neutron stars (NS) and black holes (BH). Compact binary systems of our interest are coalescing pairs that are only tens of minutes away from the final merger, i.e., from the formation of a resulting black hole or, possibly, from another spectacular event, such as a gamma-ray burst. The important question is how many such close binaries exist in our Galaxy and at cosmological distances. This determines the event rate, that is, the number of coalescence events that are likely to occur in a given cosmological volume during, say, a 1 year of observations. The event rate is partially constrained by the existing observations of binary pulsars, but its evaluation mostly relies on the numerical modelling of diverse evolutionary tracks of massive binary stars.

The event rates are sensitive to several evolutionary parameters, one of which is the so-called kick velocity parameter $`w_0`$. This parameter enters the distribution function for velocities that may be imparted on a neutron star or a black hole at the moment of its formation. It is believed that $`w_0`$ lies in the interval (200 - 400)$`km/sec`$. In Fig. 1 we plot the NS+NS, NS+BH, and BH+BH merging rates as functions of $`w_0`$. This calculation has been performed for a typical spiral galaxy of the mass $`10^{11}M_{}`$ and under a reasonable choice of other evolutionary parameters .

It follows from this graph that the NS+NS rate is expected to be at the level $`3\times 10^5year^1`$, whereas the BH+BH rate is, at least, one order of magnitude lower, i.e., $`3\times 10^6year^1`$. These rates for a typical galaxy, $`_G`$, define the rates for a given cosmological volume, $`_V`$, which includes many galaxies. When deriving $`_V`$, we have used a conservative estimate for the baryon content of the Universe, and arrived at the relationship

$$_V=0.1_G\left(\frac{r}{1Mpc}\right)^3.$$

Thus, within the volume of radius $`r=100Mpc`$, and during 1 year, one expects that there will be 3 of the NS+NS coalescences and only 0.3 of the NS+BH or BH+BH coalescences. The increase of the radius $`r`$ to $`r=200Mpc`$ increases the volume and the event rates by the factor 8.

The calculated event rates and gravitational wave amplitudes emitted by individual binaries should be compared with the instrumental sensitivities. The important fact is that the mass of a typical neutron star is 1.4 $`M_{}`$, whereas the mass of a typical stellar black hole is (10 - 15) $`M_{}`$. Thus, a pair of black holes is a more powerful source of gravitational waves than a pair of neutron stars. The black hole binaries are less numerous than the neutron star binaries, but they can be seen, by a given instrument, from much larger distances. The observationally required situation is when the detector’s sensitivity allows one to see a few events per year, and at the signal to noise ratio $`S/N`$ level of, at least, 2. Thus, the problem now is to find out the $`S/N`$ for binaries of a given total mass $`M`$.

In Fig. 2 we plot  the $`(S/N)`$ as a function of $`M`$. We have used the sensitivity curves of the initial interferometers in LIGO, VIRGO and GEO. It is assumed that the binaries consist of equal mass components, so the total mass of a NS+NS system is around 2.8$`M_{}`$, while the total mass of a BH+BH system is around (20 - 30)$`M_{}`$. The $`S/N`$ curves are shown for binaries placed at $`r=100Mpc`$. The increase of the distance to $`r=200Mpc`$ decreases the values of the plotted functions by the factor 2.

It is seen from this graph that the $`S/N`$ is less then 1 for the NS+NS binaries. Although the required 3 mergers of neutron stars do occur in the volume with the radius $`r=100Mpc`$, these sources cannot be regarded as detectable by the initial interferometers. On the other hand, for a pair of black holes with the total mass $`M=(2030)M_{}`$, the $`S/N`$ can be at the level of 4, or even larger (in the LIGO detectors). The number of the BH+BH events in the discussed volume is too small, though. However, the $`S/N`$ is still at the level of 2, or better than 2, if the binary black hole is placed at the twice as large distance, i.e. $`r=200Mpc`$. As was shown above, one can expect about 2 - 3 BH+BH events in this larger volume. This is why we believe  that the coalescing black holes will probably be the first sources detected by the initial laser interferometers. It appears that other possible sources of gravitational waves are not likely to be detected at all by these instruments. The situation will become much better, and for all known sources, when the advanced interferometers become operational. And, of course, there exists a nonzero probability that the first detected signals will be totally unexpected to everyone, or something what the authors of Ref.  did not even include in their list of options. For additional information on astrophysical sources of gravitational waves and their detection, see the comprehensive contributions by Ferrari  and Vallisneri .

Relic gravitational waves are not expected to be accessible to the initial laser interferometers. The prospects are much better, though, for the advanced ground-based instruments and for the space-based interferometer LISA . At very large scales, relic gravitational waves should also manifest themselves through the anisotropy and polarisation of the cosmic microwave background radiation (CMBR). It is fundamentally important to discover and study relic gravitational waves. They are the only clue to the processes that have taken place in the extreme conditions of the very early Universe and, indeed, they are the only clue to the origin of the Universe itself. Without observing and studying relic gravitational waves, there exists hardly any other way to test fundamental physical theories, such as superstrings and $`M`$-theory. The understanding that the expansion rate of the very early Universe could be grossly different from the one governed by the radiation-dominated fluid exists for a long time, at least, from the study of Sakharov . In the recent years, however, there exists the unjust tendency to replace the notion and all the complexities of the very early Universe by a single word - “inflation”. The reminder of my contribution is devoted to relic gravitational waves and the associated issues of cosmological perturbations of quantum-mechanical origin.

## 2 The inevitability of relic gravitational waves

The theory of almost all astrophysical sources of gravitational waves is based on classical physics. In contrast, the theory of relic gravitational waves involves the elements of quantum physics. It is the power and beauty of relic gravitational waves that their existence relies only on general relativity and basic principles of quantum field theory, such as the concept of zero-point quantum oscillations. To understand the generation of relic gravitons it is convenient to start from a classical picture. Imagine a classical gravitational wave propagating in a nonstationary universe. The nonlinear character of the Einstein equations provides the coupling of the gravitational wave field to the smooth gravitational field of the evolving universe. The gravitational wave itself can be arbitrarily weak, so that the quadratic terms of the wave field can be neglected. However, the nonlinearity of the Einstein equations leads to the presence of the time-dependent background gravitational field in the coefficients of the wave equation. The background gravitational field plays the role of the pump field capable of increasing the amplitude of the propagating gravitational wave . In fact, the traveling wave is being amplified, with the accompanying appearance of the wave traveling in the opposite direction. Together, they form an (almost) standing wave. This type of coupling of a wave-field to gravity is not something that is true for any wave-field. On the contrary, it is quite unusual. For instance, it does not take place for electromagnetic waves.

The quantum element enters the picture when we realise that if the initial gravitational waves represent only the inevitable zero-point quantum oscillations, the amplification mechanism should be at work. In precise technical terms, the quantum-mechanical Schrödinger evolution transforms the initial vacuum state of gravitational waves into a multi-quantum state known as squeezed vacuum quantum state . It is the variance of phase that gets strongly diminished (squeezed) while the mean number of quanta and its variance get strongly increased. Of course, for this effect to be numerically considerable, the coupling of the wave field to the pump field should be considerable. The coupling is regulated by the ratio of the wave period to the characteristic time-scale of variations of the background gravitational field. If the pump field varies very slowly (adiabatically), this ratio is very small, and practically nothing happens to the wave. In other words, during the adiabatic regime, the vacuum state remains the vacuum state. In the opposite (superadiabatic) regime, the relative increase (as compared with adiabatic behaviour) of the wave amplitude, and the increase of the mean number of particles, is most dramatic. In the cosmological setting, the amplification mechanism is effective during the interval of time when the wavelength is comparable with, and longer than, the cosmological Hubble radius. Thus, the generating process (superadiabatic (parametric) amplification of the zero-point quantum oscillations ) is a local and dynamical process, it does not rely on the existence of globally defined event horizons of the background space-time. Moreover, in cosmological models of practical interest, there is no horizons at all. It is also irrelevant which particular matter source drives the gravitational pump field, be it a scalar field or something else. In some of its basic aspects, the quantum-mechanical generation of cosmological perturbations is similar to the laboratory-tested generation of squeezed light.

The gravitational field of a homogeneous isotropic universe is given by the metric

$$\mathrm{d}s^2=c^2\mathrm{d}t^2+a^2(t)g_{ij}\mathrm{d}x^i\mathrm{d}x^j=a^2(\eta )[\mathrm{d}\eta ^2+g_{ij}\mathrm{d}x^i\mathrm{d}x^j],$$
(1)
where the scale factor $`a(t)`$ (or $`a(\eta )`$) is driven by a matter distribution with some effective (in general, time-dependent) equation of state. The perturbed gravitational field can be written (for simplicity, in a spatially-flat universe) as

$$\mathrm{d}s^2=a^2(\eta )[\mathrm{d}\eta ^2+(\delta _{ij}+h_{ij})\mathrm{d}x^i\mathrm{d}x^j],$$
(2)
$`h_{ij}(\eta ,𝐱)={\displaystyle \frac{𝒞}{(2\pi )^{3/2}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}d^3𝐧{\displaystyle \underset{s=1,2}{}}\stackrel{s}{p}_{ij}(𝐧){\displaystyle \frac{1}{\sqrt{2n}}}[\stackrel{s}{h}_n(\eta )e^{i𝐧𝐱}\stackrel{s}{c}_𝐧+\stackrel{s}{h}_n^{}(\eta )e^{i𝐧𝐱}\stackrel{s}{c}_𝐧^{}].`$ (3)
By expanding the functions $`h_{ij}(\eta ,𝐱)`$ over spatial Fourier harmonics $`e^{\pm i\mathrm{𝐧𝐱}}`$ we reduce the perturbed dynamical problem to the evolution of mode functions $`\stackrel{s}{h}_n(\eta )`$ for each mode $`𝐧`$. Two polarisation tensors $`\stackrel{s}{p}_{ij}(𝐧),s=1,2`$ have different forms depending on whether they represent gravitational waves, rotational perturbations, or density perturbations. In the case of gravitational waves, the polarisation tensors describe two polarisation states which are often called the “plus” and “cross” polarisations.

For the quantized field, the quantities $`\stackrel{s}{c}_𝐧,\stackrel{s}{c}_𝐧^{}`$ are annihilation and creation operators satisfying the conditions

$$[\stackrel{s^{}}{c}_𝐧,\stackrel{s}{c}_𝐦^{}]=\delta _{s^{}s}\delta ^3(𝐧𝐦),\stackrel{s}{c}_𝐧|0=0,$$
(4)
where $`|0`$ (for each $`𝐧`$ and $`s`$) is the fixed initial vacuum state defined at some $`\eta _0`$ in the very distant past, long before the superadiabatic regime for the given mode has started. In that early era, the mode functions $`\stackrel{s}{h}_n(\eta )`$ behaved as $`e^{in\eta }`$, so that each mode $`𝐧`$ represented a strict traveling wave propagating in the direction of $`𝐧`$. In the case of gravitational waves, the normalization constant $`𝒞`$ is $`\sqrt{16\pi }l_{Pl}`$. For cosmological density perturbations, which we will also discuss below, the normalisation constant is $`\sqrt{24\pi }l_{Pl}`$.

The calculation of quantum-mechanical expectation values and correlation functions provides the link between quantum mechanics and macroscopic physics. Using the representation (3) and definitions above, one finds the variance of metric perturbations:

$$0|h_{ij}(\eta ,𝐱)h^{ij}(\eta ,𝐱)|0=\frac{𝒞^2}{2\pi ^2}_0^{\mathrm{}}n^2\underset{s=1,2}{}|\stackrel{s}{h}_n(\eta )|^2\frac{\mathrm{d}n}{n}.$$
(5)
The quantity

$$h^2(n,\eta )=\frac{𝒞^2}{2\pi ^2}n^2\underset{s=1,2}{}|\stackrel{s}{h}_n(\eta )|^2$$
(6)
gives the mean-square value of the metric (gravitational field) perturbations in a logarithmic interval of $`n`$ and is called the (dimensionless) power spectrum. The power spectrum of metric perturbations is a quantity of great observational importance. It defines the temporal structure and amplitudes of the g.w. signal in the frequency bands of direct experimental searches. It is also crucial for calculations of anisotropy and polarisation induced in CMBR by relic gravitational waves and by other metric perturations. The inevitable squeezing makes the spectrum an oscillatory function of time, which can also be thought of as a consequence of the standing-wave pattern of the generated field. At every fixed moment of time, the spectrum, as a function of the wave-number $`n`$, has many maxima and zeros.

To find the power spectrum at any given moment of time (for instance, today or at the moment of decoupling of CMBR from the rest of matter) we need to know the mode functions at those moments of time. The mode functions are governed by the Heisenberg equations of motion whose interaction Hamiltonian participates with the coupling function proportional to $`a^{}/a`$. These equations are, of course, equivalent to the perturbed Einstein equations. The second-order “master equation” describing gravitational waves is

$$\stackrel{s}{\mu }_n^{\prime \prime }+\stackrel{s}{\mu }_n[n^2\frac{a^{\prime \prime }}{a}]=0,$$
(7)
where the functions $`\stackrel{s}{\mu }_n(\eta )`$ are related to the mode functions $`\stackrel{s}{h}_n(\eta )`$ by

$$\stackrel{s}{\mu }_n(\eta )a(\eta )\stackrel{s}{h}_n(\eta ).$$
(8)
One can view Eq. (7) as the equation for a parametrically disturbed oscillator (the term in brackets is the square of its variable frequency), or as the Schrödinger equation for a particle moving in the presence of a potential barrier $`W(\eta )=a^{\prime \prime }/a`$ (while remembering that $`\eta `$ is a time coordinate rather than a spatial coordinate).

As soon as the pump field (represented by $`a(\eta )`$) is known, and since the initial conditions are fully determined, the mode functions are strictly calculable. We know relatively well the behaviour of $`a(\eta )`$ at the matter-diminated $`(m)`$ and radiation-dominated $`(e)`$ stages of expansion. However, we do not know $`a(\eta )`$ at the most important preceeding stage of evolution, which we call the initial $`(i)`$ stage. Ideally, if the full evolution of $`a(\eta )`$ had followed from some fundamental physical theory, the relic g.w. background would be uniquely determined, and we would have to live with what the fundamental theory dictates. In the absence of such a theory, we can only do calculations for various possible models and compare the results with observations. Conversely, the observational restrictions on relic gravitons, or their direct detection, is our clue to the dynamics of the very early Universe.

## 3 Calculating the power spectrum

The scale factor $`a(\eta )`$ at the matter-dominated stage, governed by whatever matter with the effective equation of state $`p=0`$, behaves as $`a(\eta )\eta ^2`$. It is convenient to write $`a(\eta )`$ in the explicit form

$$a(\eta )=2l_H(\eta \eta _m)^2,$$
(9)
where $`l_H`$ is the Hubble radius today and $`\eta _m`$ is a constant explained below. The moment of time “today” (in cosmological sense) is labeled by $`\eta =\eta _R`$ (the subscript $`R`$ denoting “reception”). It is convenient to choose

$$\eta _R\eta _m=1.$$
(10)
With this convention, $`a(\eta _R)=2l_H`$, and the wave, of any physical nature, whose wavelength $`\lambda `$ today is equal to today’s Hubble radius, carries the constant wavenumber $`n_H=4\pi `$. Longer waves have smaller $`n`$’s and shorter waves have larger $`n`$’s, according to the relationship $`n=4\pi l_H/\lambda `$. For example, the ground-based gravitational wave detectors are most sensitive to frequencies around $`30Hz3000Hz`$. The corresponding wavelengths have wavenumbers $`n`$ somewhere in the interval $`10^{20}10^{22}`$.

The matter-dominated era was preceded by the radiation-dominated era with the scale factor $`a(\eta )\eta `$. Without any essential loss of generality, we assume that the transition from $`e`$ era to $`m`$ era was instantaneous and took place at some $`\eta =\eta _2`$. The redshift of the transition is $`z_{eq}`$: $`a(\eta _R)/a(\eta _2)=1+z_{eq}`$. It is believed that $`z_{eq}`$ is somewhere near $`6\times 10^3`$. In its turn, the radiation-dominated era was preceded by the initial era of expansion, whose nature and scale factor are, strictly speaking, unknown. To simplify the analysis, and since the wave equations admit simple exact solutions in the case of power-law scale factors, we assume that the $`i`$ era, similar to the $`e`$ and $`m`$ eras, was also described by a power-law scale factor. Following the early notations , we parametrize the $`i`$ era by $`a(\eta )|\eta |^{1+\beta }`$. The transition from $`i`$ era to $`e`$ era takes place at some $`\eta =\eta _1`$ and at redshift $`z_i`$: $`a(\eta _R)/a(\eta _1)=1+z_i`$. Further analysis shows (see below) that in order to get the right amplitude of the generated perturbations, the numerical value of $`z_i`$ should be somewhere near $`10^{30}`$. Whether the $`i`$ era was governed by a scalar field $`\phi `$ (the central element of inflationary scenaria) or by something else, is irrelevant. What is relevant is the functional form of the gravitational pump field represented by the scale factor $`a(\eta )`$.

We now write the full evolution of the growing scale factor explicitly:

$$a(\eta )=l_0|\eta |^{1+\beta },\eta \eta _1,\eta _1<0,\beta <1,$$
(11)
$$a(\eta )=l_0a_e(\eta \eta _e),\eta _1\eta \eta _2,$$
(12)
$$a(\eta )=2l_H(\eta \eta _m)^2,\eta _2\eta .$$
(13)
The continuous matching of $`a(\eta )`$ and $`a^{}(\eta )`$ at the transition points determines all the participating constants in terms of $`l_H,z_i,z_{eq}`$ and $`\beta `$. Thus, we are essentially left with the two unknown parametrs: $`z_i`$ and $`\beta `$. The $`z_i`$ is primarily responsible for the amplitude of the generated perturbations, while the $`\beta `$ is responsible for the spectral slope of the metric power spectrum. If $`\beta =2`$, the generated spectrum is flat, that is, independent of the wave-number $`n`$. The flat spectrum is also called the Harrison-Zeldovich-Peebles spectrum and is characterised by the spectral index $`\mathrm{n}=1`$ . The relationship between the spectral index $`\mathrm{n}`$ and parameter $`\beta `$ is $`\mathrm{n}=2\beta +5`$. Inflationary models governed by scalar fields $`\phi `$ are uncapable of producing the “blue” power spectra, i.e. spectra with $`\mathrm{n}>1`$ ($`\beta >2`$). Indeed, $`\beta >2`$ requires the effective equation of state at the $`i`$ stage to be $`ϵ+p<0`$, but this cannot be accomodated by any scalar field whatever the scalar field potential $`V(\phi )`$ may be .

For the CMBR calculations one also needs the redshift $`z_{dec}`$ of the last scattering surface $`\eta =\eta _E`$ (with the subscript $`E`$ denoting “emission”), where the CMBR photons have decoupled from rest of the matter: $`a(\eta _R)/a(\eta _E)=1+z_{dec}`$. The numerical value of $`z_{dec}`$ is somewhere near 1000.

Since the scale factor and initial conditions for the mode functions have been strictly defined, the power spectrum (6) is unambiguously calculable. We will present the results of calculation of $`h^2(n,\eta )`$ taken at the time of decoupling $`\eta =\eta _E`$. The exact form of the power spectrum is

$$h^2(n,\eta _E)=\frac{l_{Pl}^2}{4\pi l_H^2}(1+z_{dec})n^4|B|^2g^2(x,b),$$
(14)
where

$$g^2(x,b)=\left[\rho _1(x)j_1(bx)\rho _2(x)j_2(bx)\right]^2,x\frac{n}{2\sqrt{1+z_{eq}}},b\frac{2\sqrt{1+z_{eq}}}{\sqrt{1+z_{dec}}}.$$

$$\rho _1(x)=\frac{1}{x^2}\left[(8x^21)\mathrm{sin}x+4x\mathrm{cos}x+\mathrm{sin}3x\right],$$

$$\rho _2(x)=\frac{1}{x^2}\left[(8x^21)\mathrm{cos}x4x\mathrm{sin}x+\mathrm{cos}3x\right],$$

and $`j_1(bx),j_2(bx)`$ are spherical Bessel functions. The quantity $`n^4|B|^2`$ is, in general, $`n`$-dependent and $`\beta `$-dependent, but for the case of $`\beta =2`$ it simplifies to

$$n^4|B|^2=\frac{4(1+z_i)^4}{(1+z_{eq})^2},\beta =2.$$
(15)
In Fig. 3 we show  the function $`g^2(x,b)`$ for $`b=5`$.

## 4 Detectability of relic gravitational waves

We start from a discussion of manifestations of relic gravitational waves in the CMBR anisotropies. It can be shown (see  and references there) that the angular correlation function has the universal form

$$0\left|\frac{\delta T}{T}(𝐞_1)\frac{\delta T}{T}(𝐞_2)\right|0=\underset{l=2}{\overset{\mathrm{}}{}}\frac{2l+1}{4\pi }C_lP_l(\mathrm{cos}\delta ),$$
(16)
where $`P_l(\mathrm{cos}\delta )`$ are Legendre polynomials for the separation angle $`\delta `$ between the unit vectors $`𝐞_1`$ and $`𝐞_2`$, and the multipole moments $`C_l`$ are explicitely determined by the metric mode functions.

In Fig. 4 we show  by a solid line the graph of the function $`l(l+1)C_l`$ calculated for the g.w. background with the power spectrum (14). The cosmological parameters were taken, for illustration, as $`z_{eq}=10^4,z_{dec}=10^3,\beta =2`$. The parameter $`z_i`$ is adjusted in such a manner ($`z_i=10^{29.5}`$) that the graph goes through the point $`l(l+1)C_l=6.4\times 10^{10}`$ at $`l=10`$, which agrees with observations. For comparison, the dashed line shows the same function, but calculated for the alternative (non-squeezed) background of gravitational waves of the same averaged power density. The remarkable (even if expected) result is that the background of non-squeezed (traveling) gravitational waves does not produce oscillations in the angular power spectrum $`C_l`$, whereas the background of squeezed (standing) gravitational waves does. The peaks and dips of the angular power spectrum are in close relationship with the maxima and zeros of the metric power spectrum of Fig. 3. It is unlikely that the peaks and dips caused by relic gravitational waves can be observationally revealed, but this calculation serves as a guidance for the related explanation (see below) of the actually observed peaks and dips  that lay at a higher level of the signal.

For purposes of direct experimental searches for relic gravitational waves, one needs to know the metric power spectrum in the present era. In Fig. 5 we show the piece-wise envelope of the r.m.s. quantity $`h(n,\eta _R)`$ plotted as a function of frequency $`\nu `$ in $`Hertz`$.

This particular graph is derived under the assumption that the parameter $`\beta `$ is $`\beta =1.9`$, which corresponds to the spectral index $`\mathrm{n}=1.2`$. A few years ago (see, for example, ), this value of $`\mathrm{n}`$ was quoted as the one that follows from the lower-order multipoles of the CMBR data. The most recent analyses, which include data on the large scale structure, seem to prefer a smaller value of $`\mathrm{n}`$. Although the often cited central values of $`\mathrm{n}`$ are still somewhat larger than 1, the 1$`\sigma `$ error-bars do not usually include the value $`\mathrm{n}=1.2`$. One should remember, however, that the current analyses of the data are consistent only if one assumes that about 70 percent of “matter” in the present Universe resides in the cosmological $`\mathrm{\Lambda }`$ term or in a “dark energy/quintessence”, which is not subject to gravitational clustering. It is still not clear whether this conclusion is one of the greatest discoveries about the Universe, or it may be a result of improper evaluation of the role of metric perturbations in the CMBR analysis. We will comment on this issue in the next section. In addition, one should bear in mind that the flat ($`\mathrm{n}=1`$) and “red” ($`\mathrm{n}<1`$) spectra possess a theoretical difficulty: the integral over wavenumbers $`n`$ for the mean square value of the metric (gravitational field) perturbations diverges at the lower limit $`n0`$, i.e. for very long waves. The removal of these divergencies would require some extra assumptions. In any case, we will discuss the prospects of direct detection of relic gravitational waves assuming that the graph in Fig. 5 is correct or close to correct.

According to analytical formulas and Fig. 5, the theoretically expected r.m.s. value at $`\nu =100Hz`$ is $`h_{th}=10^{25}`$. The target value of the LIGO-II  at this frequency is $`h_{ex}=10^{23}`$. The remaining gap of 2 orders of magnitude can be covered by the cross-correlation of sufficiently long strips of data from two or more detectors. The $`S/N`$ ratio will be better than 1, if the common integration time exceeds $`\tau =10^6sec`$. This does not look like a hopeless task. The situation is somewhat better in the LISA frequency interval. According to evaluations of Ref. , a 1-year observation time can provide a $`S/N`$ ratio around 3 in the frequency interval $`(2\times 10^310^2)Hz`$. These estimates do not include a possible improvement in the $`S/N`$, which can be achieved if one succeeds in finding a sophisticated method of exploiting, at these frequencies, the ever-present squeezing. This improvement would certainly be possible in a very narrow-band detector, but it is unclear whether it can ever be realised in a realistic broad-band detector (for a discussion, see  and ). At the present time, it is difficult to say at which scales the relic gravitational waves will be first discovered. This can happen in direct experiments at relatively small scales, or, indirectly, at cosmological scales, through the anisotropy and polarisation measurements of the CMBR.

## 5 Relic gravitational waves and primordial density perturbations

The superadiabatic (parametric) mechanism is always operative for gravitational waves, but its applicabilty to density perturbations requires special assumptions about properties of matter. Of course, the notion of density perturbations includes the associated metric perturbations $`h_{ij}`$ of Eq. 3; the matter and gravitational field (metric) variables are always linked by the perturbed Einstein equations. If the $`i`$ stage was governed by a scalar field $`\phi `$ coupled to gravity in a special manner (a hypothesis usually assumed in inflationary scenaria), the quantum-mechanical generation of primordial density perturbations, in addition to the inevitable generation of relic gravitational waves, becomes possible. The underlying physical theory for these two respective fields is exactly the same, and therefore the general properties of relic gravitational waves and primordial density perturbations (squeezing, standing-wave pattern, spectral slopes, numerical values of the amplitudes, etc.) are almost identical. In particular, relic gravitational waves and primordial density perturbations should provide approximately equal contributions to the lower-order multipoles (starting with the quadrupole moment $`C_2`$) of the CMBR anisotropies . Some differences arise, however, at the late stage of their evolution. In the matter-dominated era, and for wavelengths comfortably shorter than the Hubble radius, numerical values of the metric components associated with density perturbations grow in proportion to $`a(\eta )`$ (gravitational instability), whereas the amplitudes of gravitational waves decrease in proportion to $`1/a(\eta )`$ (adiabatic behaviour). This explains why the metric amplitudes of these two fields, being of the same order of magnitude at scales comparable with the prersent-day Hubble radius, are significantly different at much smaller scales. This difference makes gravitational waves subdominant in their contribution to the CMBR multipoles near the peak at $`l200`$, and makes them so small and difficult to detect at laboratory scales.

As a benchmark for the spectral index one can take $`\beta =2`$, i.e. flat spectrum. The CMBR anisotropies caused by gravitational waves are shown in Fig. 4. The lowest excited multipole is $`l=2`$ (quadrupole). In contrast, density perturbations produce all multipoles, including $`l=0`$ (monopole) and $`l=1`$ (dipole). In general, CMBR anisotropies caused by density perturbations have contributions from 3 different ingredients (for a recent detailed discussion, see ). Specifically, the ingredients are: the accompanying metric perturbations (mostly driven by the gravitationally dominant matter, which is, according to the existing views, a presureless cold dark matter (CDM)), the intrinsic temperature variations at the last scattering surface, and the velocities (with respect to the coordianate system comoving with the CDM) of the last scattering electrons.

The angular power spectrum $`l(l+1)C_l`$ is usually presented as an outcome of numerical calculations, where all contributions are being mixed up. However, the evolutionary equations, on which the numerical codes are based, allow analytical calculations for the lower multipoles, i.e. for $`l`$’s up to $`l(3040)`$. It is often claimed that the analytical limit of the employed evolutionary equations produces the lower-order multipoles $`C_l`$ according to the formula

$$C_l=\frac{const.}{l(l+1)}.$$
(17)
This formula predicts an infinitely large monopole moment $`C_0`$, and the dipole moment $`C_1`$ only a factor 3 greater than the quadrupole moment $`C_2`$; both predictions are in a severe conflict with observations. Moreover, this formula misses the quite steep growth, beginning from $`l=2`$, of the function $`l(l+1)C_l`$, which must take place for purely gravitational reasons, that is, when the metric perturbations are treated correctly .

It was demonstrated  that formula 17 is incorrect and cannot be the outcome of correct equations. Specifically, formula 17 arises only if one unjustifiably neglects 3 out of 4 terms in the exact formula (43) of the Sachs-Wolfe paper . The correct handling of metric perturbations makes the monopole moment $`C_0`$ finite and small, whereas the dipole moment $`C_1`$ becomes 5 orders of magnitude greater than the number following from formula 17, and, thus, the theoretical (statistical) $`C_1`$ comes in agreement with the observed value of $`C_1`$. It was also shown that the growth of the function $`l(l+1)C_l`$ would be unlimited, if one were allowed to extrapolate the flat spectrum to arbitrarily large wave-numbers $`n`$. This is illustrated in Fig 6.

In reality, however, the modulating (transfer) function of the metric perturbations bends the angular spectrum down and introduces a series of peaks and dips, analogous to the peaks and dips caused by gravitational waves in Fig. 4. We will discuss peaks and dips shortly, but first we shall say a few words about the importance of theoretical calculations of the expected (statistical) value of the dipole moment $`C_1`$. Because of statistical properties of the generated perturbations, all statistical multipoles $`C_l`$ are independent of the observer’s position.

The derivation of formula 17 takes into account the perturbations of all wavelengths. The integral over $`n`$, which gives rise to this formula, is formally extended from $`n=0`$ to arbitrarily large $`n`$’s. If this formula were correct, it would be applicable all over the linear regime of cosmological perturbations and could possibly break down only due to nonlinearities at relatively short scales. The gap in 5 orders of magnitude between the statistical prediction of this formula and the actually observed $`C_1`$ would need to be attributed either to extremely improbable realisation of the perturbed gravitational field, or to the quite local conditions where nonlinearities set it and where the linear approximation is not valid. In the latter case, numerical values of $`C_1`$’s observed by inhabitants of the not-too-distant galaxies would not need to have anything in common with each other and with the $`C_1`$ observed by us. In contrast, calculations of Ref.  show that numerical value of the statistical $`C_1`$ is most sensitive to scales where the metric power spectrum bends down. In other words, the corresponding integral is primarily accumulated at scales of the order of $`(100200)Mpc`$, i.e. well within the linear regime. This result comes about from the “gradient of gravitational potential” terms (see also ), which are totally neglected in the derivation of formula 17. This means that the galaxies within, at least, the radius of $`(100200)Mpc`$, and, at least, the field galaxies, are expected to have approximately equal observed $`C_1`$’s. This number should also be of approximately the same value as our own $`C_1`$ (as soon as our measured $`C_1`$ turns out to be close to the statistical $`C_1`$). In principle, this picture is observationally testable. Of course, we cannot place ourselves at different galaxies and measure the $`C_1`$ there. But, it is not excluded that the CMBR dipoles seen at other galaxies can be evaluated indirectly, possibly with the help of the Sunyaev-Zeldovich effect. It would be very interesting if these observational evaluations could be done. In any case, the calculation of the theoretical (statistical) $`C_1`$ is as important as the calculation of other multipoles. It would be embarassing if some of currently popular cosmological models would predict the statistical $`C_1`$ way out of what is actually observed by us and by observers in other galaxies (assuming that we can evaluate their dipoles).

It is argued in Ref.  that the (squeezed) metric perturbations associated with density perturbations are mostly responsible for the recently observed peaks and dips in the angular power spectrum $`l(l+1)C_l`$. This interpretation is different from the concept of “acoustic peaks” based on the existence of the plasma sound waves at the last scattering surface. The rigorous evolution of the quantum-mechanically generated density perturbations through the radiation-dominated and matter-dominated eras, allows one to find the metric power spectrum, defined by Eq. 6, at the last scattering surface $`\eta =\eta _E`$. The analytical formula for $`h^2(n,\eta _E)`$ is given by

$$h^2(n,\eta _E)=\frac{𝒞^2}{2\pi ^2}\frac{n^4|B|^2(1+z_{eq})}{48l_H^2}\left(\frac{\mathrm{sin}y_2}{y_2}\right)^2\frac{(30020p^2y_2^2+p^4y_2^4)}{200},$$
(18)
where

$$y_2=\frac{n}{2\sqrt{3}\sqrt{1+z_{eq}}},\mathrm{and}p\frac{2\sqrt{3}\sqrt{1+z_{eq}}}{\sqrt{1+z_{dec}}}.$$

The graph of the oscillating function

$$f^2(x,p)\left(\frac{\mathrm{sin}x}{x}\right)^2\left[30020p^2x^2+p^4x^4\right],$$

where $`xy_2`$ and, for illustration, we take $`p=8`$, is shown in Fig. 7.

The oscillations in the metric power spectrum are a direct consequence of squeezing (standing-wave pattern) of the primordial density perturbations. These oscillations find their reflection in the peaks and dips of the angular power spectrum. The main difference with the concept of “acoustic peaks” lays in the two facts: the amplitudes of metric perturbations are larger than the amplitudes of matter density variations in the photon-electron-baryon fluid, and the wave-number periodicity of metric perturbations is governed by a sine function, instead of a cosine function for acoustic perurbations. It is argued  that the actually observed features are better described by the derived analytical formulas than by the concept of “acoustic peaks”.

The study of density perturbations is directly relevant to the problem of relic gravitational waves. If cosmological perturbations are indeed generated quantum-mechanically (what the present author believes is likely to be true), many properties of primordial density perturbations and relic gravitational waves should be common. By studying density perturbations (for instance, by proving the presence of squeezing) we will be more confident about the properties of relic gravitational waves.

One should be aware that the story of quantum-mechanically generated cosmological perturbations is dramatically different in inflationary scenario. Inflationists love to speak about “stretching the quantum fluctuations of the inflaton to cosmological scales”. As a result of this “stretching”, they manage to produce the arbitrarily large amplitudes of today’s density perturbations. Indeed, the “standard result” of inflationary scenario predicts the infinitely large density perturbations, in the limit of the flat spectrum ($`\mathrm{n}=1`$), through the set of evaluations: $`\delta \rho /\rho h_S\zeta H^2/\dot{\phi }V^{3/2}(\phi )/V^{}(\phi )1/\sqrt{1\mathrm{n}}`$. \[The “quantization”, as understood by inflationists, of the necessarily coupled system of scalar field and metric perturbations is such, that the so-called Bardeen’s gauge-invariant quantity $`\zeta `$ is arbitrarily large already at the beginning of the amplifying (superadiabatic) regime, and then, being “conserved”, this arbitrarily large number translates to the end of the amplifying regime, and to the present time.\] Of course, this “standard result” is in full disagreement not only with the theoretical quantum mechanics, but with available observations too, as long as the error-boxes of the observationally derived spectral index $`\mathrm{n}`$ are centered at $`\mathrm{n}1`$ and include the value $`\mathrm{n}=1`$. However, by composing the ratio of the gravitational wave amplitude $`h_T`$ to the predicted divergent amplitude of the scalar metric perturbations $`h_S`$ (the so called “consistency relation”: $`h_T/h_S\sqrt{1\mathrm{n}}`$), inflationary theorists substitute their prediction of arbitrarily large density perturbations for the claim that it is the amount of gravitational waves that should be zero, or almost zero, at cosmological scales and, hence, down to laboratory scales. This claim has led to many years of mistreatment of a possible g.w. contribution to the CMBR data. The “standard inflationary result” is maintained by inflationists and their followers until now. For instance, the recent review paper  claims that the initial spectrum of gravitational waves is “constrained to be small compared with the initial density spectrum”. It is only in a few recent papers (for example, ) that the inflationary “consistency relation” is not being used when analyzing the CMBR and large scale structure observations, with some interesting conclusions. The inflationary claim of the “negligible” amount of gravitational waves influences also the LISA community. It is sometimes opined that the LISA parameters should be better optimised for detection of g.w. noise from the multitude of binary white dwarfs, rather than for a search for relic gravitational waves. There seems to be some logic behind this opinion. If one hears so often about “excellent agreement” of inflationary predictions with observations, and if the central of these predictions is “negligible” amount of relic gravitational waves, there is no need to bother about them. Many of the contemporary writers manage to do both: to praise inflationary predictions and to speculate about “testing inflation” with the help of gravitational waves, including the CMBR polarisation measurements and post-LISA missions. Apparently, they plan to do this at the signal to noise ratio level of $`10^{20}`$, or so, for “inflation predicted” gravitational waves. For some reason, those authors do not suggest to test inflationary predictions right from the predicted divergent amplitudes of density perturbations. Sometimes, inflationary calculations based on the “standard result” and “consistency relation” arrive at gravitational-wave amplitudes which are only 1-2 orders of magnitude smaller than, or even comparable (for the spectral index $`\mathrm{n}`$ sufficiently smaller than 1) with the metric amplitudes of density perturbations. Of course, this does not make the fundamentally wrong “standard result” a correct one. One can find more reading about the “standard inflationary result” and its critique, including the earlier claim of inflationists that their divergent amplitudes of density perturbations are justified by the (incorrect and now abandoned) concept of “big amplification during reheating”, in the end of Sec. 6 of Ref.  and in references there.

## 6 Conclusions

It is important to understand that the existence of relic gavitational waves relies only on general relativity and basic principles of quantum field theory. The detection of relic gravitational waves will be a direct probe of the expansion rate of the very early Universe, and not a test of some extra assumptions, such as the dominance of a scalar field or the form of its potential $`V(\phi )`$. If the observed large-scale CMBR anisotropy $`\delta T/T`$ is indeed caused by cosmological perturbations of quantum-mechanical origin, the contribution of relic gravitational waves to $`\delta T/T`$ must be considerable. The extrapolation of the g.w. metric power spectrum to shorter wavelengths shows that, depending on the still existent uncertainty in the spectral index, relic gravitational waves may be measurable by the advanced ground-based or space-based laser interferometers. The ever-present squeezing must manifest itself in the periodic structure of the metric power spectra and, at cosmological scales, in the oscillatory behaviour of the CMBR multipoles $`C_l`$ as a function of $`l`$. When evaluating the scientific importance and prospects of discovery of relic gravitational waves, it is necessary to respect the conclusions of quantum mechanics and general relativity, and not the statements of popular fallacies.

## Acknowledgments

I am grateful to my coauthors whose work was reviewed in this contribution.


format: text

text:
# The 12C/13C ratio in the Planetary Nebula NGC 3242 from Hubble Space Telescope STIS observationsBased on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy (AURA), Inc., under NASA contract NAS 5–26555.

## 1 Introduction

In principle, the abundance of 3He in our Galaxy can be used to test the predictions of Big-Bang nucleosynthesis and to provide constraints on the baryon density of the Universe (e.g. Bania et al., 2002). However, the usefulness of 3He as a cosmic baryometer is limited by a poor understanding of its evolution in our Galaxy, one of the major unsolved mysteries in the field of cosmic abundances. The problem was first recognized by Rood et al. (1976) who noted the strong discrepancy between the theoretical yields of low-mass stars, $`X(^3\mathrm{He})10^3`$ for a 1 $`M_{}`$ star (Iben, 1965), and the measured abundances in Hii regions, $`X(^3\mathrm{He})10^5`$ (Balser et al., 1999a). Despite the increase in the number and quality of observational constraints and in the accuracy of stellar evolutionary calculations over the past 25 years, the 3He problem is still with us today (e.g. Tosi, 2001).

An interesting solution to the puzzle of 3He based on stellar physics has been suggested by various groups (Hogan, 1995; Charbonnel, 1995). Accordingly, a non-standard mixing mechanism acting during the RGB and/or AGB phases of stars of mass up to $`2M_{}`$ can effectively suppress the production of 3He. In this way, low-mass stars are not 3He producers and may even become net 3He destroyers, lessening the tension between the results of chemical evolution models and observations. If this mechanism is indeed at work in low-mass stars, an unavoidable consequence is that the ratio of 12C/13C in the ejecta of planetary nebulae (PNe) should be much lower than in the standard case. For a $`1M_{}`$ star, the predicted ratio is $`5`$ against the standard value of 25–30 (Boothroyd & Sackmann, 1999). Thus, it is very important to have an accurate measure of the carbon isotopic ratio in those PNe where a high 3He abundance has been determined. Should these nebulae show a 12C/13C ratio close to 25–30, then no modifications to the standard scenario would be required. Otherwise, one has to invoke another selective process (circulation, diffusion etc.) that operates on 13C but not on 3He.

NGC 3242 is the only PN with a firm detection of the 3He+ hyperfine line at 8.67 GHz. The initial discovery indicated a value of 3He/H $`10^3`$ by number (Rood et al., 1992). Subsequent observations confirmed the detection and refined the estimate of the abundance to lower values, 3He/H $`=(25)\times 10^4`$ (Balser et al., 1999b). In both cases, the resulting values of 3He agree with the predictions of standard stellar nucleosynthesis models (Rood et al., 1976).

In an attempt to derive the 12C/13C ratio in PNe, we observed at millimeter wavelengths a sample of 28 objects, determing the CO isotopic ratio in 14 of them and obtaining robust upper limits in 6 other nebulae (Palla et al., 2000). Unfortunately, no CO emission was detected in NGC 3242, so that we could not draw any conclusion about the occurrence of mixing in this PN. Following Clegg et al. (1997), we decided to adopt the complementary approach of observing the C III\] multiplet near 1908 Å, exploiting the unique spectroscopic capabilities of HST/STIS at ultraviolet wavelengths. Clegg et al. (1997) measured for the first time the carbon isotopic ratio in the ionized gas of NGC 3918 and SMC N2, using the GHRS aboard the HST. The big advantage of this method over observations in the millimeter band is that the simultaneous presence of the 12C and 13C lines in this spectral region allows a direct estimate of the carbon isotopic ratio without having to worry about problems such as isotopic fractionation and optical depth effects that make the derivation of abundances from CO lines more uncertain (Sahai et al., 1994; Schöier & Olofsson, 2000).

## 2 Observations and Results

NGC 3242 is a multiple-shell, attached halo PN, located at a distance of about 1 kpc (Stanghellini & Pasquali, 1995). The inner and outer shells have diameters of 15′′ and 29′′, respectively. The inner shell (or rim) is characterized by the presence of two prominent spots in O I, O III, and N II. The post-AGB mass of the central star is estimated to be $`M_{\mathrm{CS}}=(0.56\pm 0.01)`$ $`M_{}`$, corresponding to a main-sequence mass of $`M_{\mathrm{MS}}=(1.2\pm 0.2)`$ $`M_{}`$ (Galli et al., 1997). The metallicity of NGC 3242 is $`ZZ_{}/2`$ (Barker, 1985).

We present high quality STIS observations of 12C and 13C in the ionized gas of NGC 3242 by using the C III\] multiplet near 1908 Å. The multiplet has two lines, 3P$`{}_{1}{}^{}{}_{}{}^{0}`$1S0 at 1908.7 Å and 3P$`{}_{2}{}^{}{}_{}{}^{0}`$1S0 at 1906.7 Å which are intercombination and magnetic quadrupole transitions, and a third line, $`{}_{}{}^{3}{}_{1/2}{}^{}`$P$`{}_{0}{}^{}{}_{}{}^{0}`$$`{}_{}{}^{1}{}_{1/2}{}^{}`$S0 at 1909.6 Å, which is generally completely forbidden except in the case of the 13C atom as a result of the non zero nuclear spin. The first two lines have different transition probabilities and their ratio has been used to determine the electron density in PNe and other low-density objects. The third line is weaker than the others (typically, less than one percent), and originates only in 13C.

We observed NGC 3242 with STIS on the HST over 5 orbits in March 2001 for a total of 13000 seconds of integration time. All of the spectra were obtained through a 52$`{}_{}{}^{\prime \prime }\times 0`$$`^{\prime \prime }.`$1 slit to maximize throughput without a significant loss in spectral resolution. The slit was positioned to cut through the two conspicuous spots of the rim, avoiding the central star. We used the MAMA detectors and G230M gratings whose range well encompasses the spectral region of interest, and whose pixel scale (0.09 Å per pixel) is adequate for our observations. The prime tilt was centered at 1933 Å, achieving a spectral resolution of 0.36 Å ($`25`$ km s-1) which is enough to resolve the C III\] line at 1909.6 Å from the stronger C III\] line at 1908.7 Å.

The spectra were recalibrated using the STIS pipeline and the best reference files available in the archive as of November 2001. The individual spectra from each orbit were independently calibrated in wavelength and flux, and resampled to a linear wavelength scale (retaining the same average dispersion). Careful examination of the spectra from each orbit reveals no discernible change in flux over time or as a function of position on the detectors. We therefore averaged the spectra (weighted by exposure time) to obtain the final versions of the G230M spectra. The spectrum was then obtained by coadding $`250`$ rows selected from the highest surface brightness regions coinciding with the two bright spots of the rim.

The resultant spectrum and its residuals are shown in Figure 1. The 12C lines at 1906.7 Å and 1908.7 Å have fluxes of $`(5.36\pm 0.04)\times 10^{13}`$ and $`(3.91\pm 0.04)\times 10^{13}`$ erg cm-2 s-1 arcsec-2, respectively, yielding a ratio of 1.37$`\pm 0.02`$. For the 13C line at 1909.6 Å, we have obtained a $`1\sigma `$ upper limit of $`1.5\times 10^{15}`$ erg cm-2 s-1 arcsec-2. Thus, C III\] 1909.6/1906.7$`<2.8\times 10^3`$.

The 12C lines cannot be described with gaussians profiles because they present faint blue and red wings. As shown in Table 1, we fit the 12C lines with two gaussian components with the same central wavelength but different fluxes and line widths. The bright narrow component has a FWHM of 52 km s-1, significantly larger than the instrumental resolution of the G230M grating ($`30`$ km s-1); the broadening is probably due to the presence of two unresolved components coming from the receding and approaching edges of the nebula. The faint wings can be described with a gaussian profile with 290 km s-1 FWHM, indicating the presence of fast moving gas.

Following Clegg et al. (1997), the 1906.7/1908.7 intensity ratio yields the electron density, while the 1909.6/1908.7 line ratio is used to determine the 12C/13C abundance ratio (see their Fig. 2). For a value of $`1906.7/1908.7=1.4`$, we obtain an electron density $`n_e=8.9\times 10^3`$ cm-3. This value is larger than that reported by Stanghellini & Kaler (1989) from a variety of forbidden lines of lower ionization ($`n_e=2.7\times 10^3`$ cm-3), but still below the threshold for collisional quenching of the 13C line ($`n_e2\times 10^4`$ cm-3). Then, the expected 1909.6/1908.7 ratio computed assuming equal abundances of 12C and 13C is 0.15. Using the observed limit on $`1909.6/1906.7<2.8\times 10^3`$, we infer a lower limit to the carbon isotopic ratio of 12C/13$`>0.15/(1.4\times 2.8\times 10^3)38`$.

In Figure 2, we compare the observed 3He/H (Balser et al., 1999b) and 12C/13C ratio (this work) of NGC 3242 with the predictions of standard and non standard models for the advanced phases of the evolution of stars in the mass range 1 to 3 $`M_{}`$. Note that the largest discrepancy between the models occurs below $`2M_{}`$, i.e for stars that do experience the helium flash. It is evident that standard models (left panels) reproduce very well the results in both cases, whereas the inclusion of extra-mixing processes fails by factors between 3 and 5. As noted above, the metallicity of NGC 3242 is lower than solar by about a factor of 2. For standard models, the 3He yields are rather insensitive to $`Z`$, as shown by the calculations of Dearborn et al. (1996) (see their Fig. 1). Similarly, a 50% decrease of $`Z`$ produces only a $``$20% increase of the 12C/13C ratio (e.g. van den Hoeck & Groenewegen, 1997).

Considering the sensitivity of the non-standard models to metallicity, the effect is shown in Fig. 2 where we plot two curves for $`Z=0.02`$ and $`Z=0.007`$. At lower metallicities, the depth of the convective zone is larger, implying higher temperatures and more efficient nuclear burning at the base of the convection region. Therefore, the predicted 3He abundance and 12C/13C ratio are lower than in the solar case, and the discrepancy with the observed values becomes even larger.

## 3 Discussion

The HST observations reported here provide an independent support to the 3He observations at 9.8 GHz and allow to draw firm conclusions about the nucleosynthesis history of NGC 3242. The lack of the 13C line and the presence of the 3He line in the spectrum indicate that the progenitor star did not undergo a phase of deep mixing during the late stages of its evolution. At first glance, this result may suggest that there is no need to revise standard nucleosynthesis models, since Fig. 2 shows that there is excellent agreement between their predictions and the values derived from the observations of both 3He and 12C/13C . However, the question is whether NGC 3242 is a special case among PNe.

The existence of isotopic anomalies in low-mass stars (below $`3M_{}`$) is well documented for red giant stars (e.g. Charbonnel et al., 1998; Gratton et al., 2000). From the analysis of a large sample of field and cluster RGB stars, Charbonnel & do Nascimento (1998) conclude that more than 90% of them present 12C/13C ratios inconsistent with the results of standard models. The case for AGB stars is more controversial. Evidence for mixing comes from data on oxygen and carbon isotopic ratios (Wasserburg et al., 1995) and on 7Li abundances (Charbonnel & Balachandran, 2000), but its statistical significance is hard to assess. For example, J-type carbon stars clearly show 12C/13C ratios smaller than 10, well below the theoretical predictions (e.g. Abia & Isern, 1997; Ohnaka & Tsuji, 1999). Their progenitors appear to be low-mass stars with $`M23M_{}`$, suggesting the operation of mixing processes (Abia & Isern, 2000). However, J-type represent only $``$15% of the whole C-star population. The occurrence of mixing in this class of stars finds support from the analysis of $`s`$-process elements in a sample of C(N)-type stars (Abia et al., 2001). These stars account for the majority of the C-stars and have values of 12C/13C greater than 10 but less than the solar value (89), consistent with model results only for $`M1.5M_{}`$ (cf. Fig. 2). For $`M1.5M_{}`$, standard models predict a sharp rise of the 12C/13C ratio (up to $``$100), which is quite sensitive to metallicity effects. The reduction of this ratio induced by the occurrence of deep mixing in the AGB phase cannot be quantified at present, and detailed modelling is obviously needed to clarify the situation.

Owing to the lack of non standard AGB yields, we show in the bottom right panel of Fig. 2 the results for RGB deep mixing. How representative are these numbers for the more advanced phases? From our previous estimate of the 12C/13C ratio in 20 PNe from CO line emission, we have found that the majority of them have 12C/13C ratios between 10 and 25 for masses 1.5 to 3 $`M_{}`$, These values are lower than those expected from standard AGB models, even in the case where the 12C/13C value was very low ($``$10) at the base of the AGB, as a result of deep mixing during the RGB. The strong 12C production by third dredge-up along the AGB raises the isotopic ratio to much higher values (Abia et al., 2001). Therefore, the overall evidence seems to indicate that active mixing occurs also in the early AGB phase. Returning to the case of NGC 3242, we may conclude that it must belong to the small group of “standard” low-mass stars, consisting of $``$10% of stars with $`M2`$ $`M_{}`$.

This point is of direct relevance to the 3He probem in the context of Galactic chemical evolution. A very small percentage of low-mass stars producing large quantities of 3He is required to explain the 3He abundances measured in the pre-solar nebula, in the local ISM, and in Hii regions (Palla et al., 2000). All these objects show 3He/H ratios around $`10^5`$, lower than that measured in NGC 3242 by at least an order of magnitude.

The upper panel of Figure 3 shows the predictions of chemical evolution models assuming standard and non standard stellar nucleosynthesis and the most recent abundances measured in the solar system and local ISM. The model computed with the extra-mixing prescriptions of Sackmann & Boothroyd (1999) for 90% of low-mass stars fits very well the observations, whereas the standard model is systematically higher.

The failure of standard models is more dramatic when considering the Galactic radial distribution of 3He , as shown in the bottom panel of Fig. 3. The mixing hypothesis improves significantly the agreement, and provides an upper envelope to the distribution of simple Hii regions, i.e. those which are well described by homogeneous spheres (Bania et al., 2002). Bania et al. have argued that only these objects provide sufficiently accurate 3He/H estimates for comparison with Galactic evolution models. Since simple Hii regions have 3He abundances systematically lower than the average, the implication is that almost 100% of low-mass stars should undergo non-standard nucleosynthesis. In that event, NGC 3242 would stand as a unique object in the Galaxy.

The research of D.G., F.P. and M.T. has been supported by grants COFIN 1998-MURST and COFIN 2000-MURST at the Osservatorio Astrofisico di Arcetri and Osservatorio Astronomico di Bologna.


text:
# Responses of Ensemble Neurons to Spike-Train Signals with Independent Noises: Stochastic Resonance and Spike Variability

## I INTRODUCTION

It has been controversial how neurons communicate information by spikes -. Much of debates on the nature of the neural code has been mainly focused on the two issues. The first issue is whether information is encoded in the average firing rate of neurons (rate code) or in the precise firing times (temporal code). Andrian  first noted the relationship between neural firing rate and stimulus intensity, which forms the basis of the rate code. Actually firing activities of motor and sensory neurons are reported to vary in response to applied stimuli. In recent years, however, the alternative temporal code has been proposed in which detailed spike timings are assume to play an important role in information transmission: information is encoded in interspike intervals (ISIs) or in relative timings between firing times of spikes -. Indeed, experimental evidences have accumulated in the last several years, indicating a use of the temporal coding in neural systems -. Human visual systems, for example, have shown to classify patterns within 250 ms despite the fact that at least ten synaptic stages are involved from retina to the temporal brain . The transmission times between two successive stages of synaptic transmission are suggested to be no more than 10 ms on the average. This period is too short to allow rates to be determined accurately.

The second issue is whether information is encoded in the activity of single (or very few) neurons or that of a large number of neurons (population or ensemble code). The population rate code model assumes that information is coded in the relative firing rates of ensemble neurons, and has been adopted in the most of the theoretical analysis . On the contrary, in the population temporal code model, it is assumed that relative timings between spikes in ensemble neurons may be used as an encoding mechanism for perceptional processing -. A number of experimental data supporting this code have been reported in recent years -. For example, data have demonstrated that temporally coordinated spikes can systematically signal sensory object feature, even in the absence of changes in firing rate of the spikes .

The strong criticism against the temporal code is that spikes are vulnerable to noise while the rate code performs robustly in the presence of noise but with limited information capacity. It is well known that although firings of single neocortical neurons in vitro are precise and reliable, those in vivo are quite unreliable . This is due to noisy environment in viro, which makes the reliability of neurons firings worse. In recent years, however, much studies have been made for the stochastic resonance (SR)  in which information transmission of signals is enhanced by background noises, against our conventional wisdom. SR in a neural system has been theoretically investigated -. The transmission fidelity for weak external signals, which is evaluated by the signal-to-noise ratio (SNR) or the peak height of the interspike-interval (ISI) distribution, is enhanced by added noises. SR is supported in some physiological experiments for biological systems such as crayfish , cricket  and rat .

Although SR itself is a genetic phenomenon, its detailed character is determined by the three important factors: (a) kinds of systems (neurons), (b) input signal and (c) noises. As for the first factor of (a) neurons, SR in single neurons has been studied by using various theoretical models such as the integrate-and-fire (IF) model -, the FitzHough-Nagumo (FN) model - and the Hodgkin-Huxley (HH) model . SR in coupled or ensemble neurons has been also investigated by using the IF model -, FN model - and HH model -. The transmission fidelity has maximum when the noise magnitude or the coupling strength is changed. It has been pointed out that the transmission fidelity of ensemble neurons for sinusoidal inputs with independent white noises is improved as the size of an ensemble is increased .

As for the second factor of (b) input signals, most of theoretical studies have been made for analog inputs with periodic (mostly sinusoidal) or aperiodic amplitude modulation . This is because these studies have been motivated by a fact that peripheral sensory neurons play a role of transducers receiving analog stimulus and emitting spikes. In central neural systems, however, cortical neurons are reported to play a role of data-processors receiving and transmitting spike trains . There are only a few theoretical studies on SR for spike-train inputs -. The response of single IF neurons to coherent spike-train (ST) inputs is shown to be enhanced by an addition of weak ST noises characterized by the Poisson   or gamma distribution . Quite recently, the present author has studied SR of ensemble HH neurons for transient spike-train inputs with independent Gaussian noise by using the wavelet analysis .

As for the third factor of (c) noises, it has been reported that noises are ubiquitous in neural systems. The origin of these noises is not clear at the moment. We may suppose, however, several conceivable origins of noises: (i) cells in sensory neurons are exposed to noises arising from the outer world, (ii) ion channels of the membrane of neurons are known to be stochastic , (iii) the synaptic transmission yields noises originating from random fluctuations of the synaptic vesicle release rate , and (iv) synaptic inputs include leaked currents from neighboring neurons . Most of existing studies on SR have simulated noises of the items (i)-(iii) by the Gaussian white noise - - or Orstein-Uhrenbeck (OU) noise  . ST noise is employed in Refs.- taking account the item (iv). In our study, we will include the ST and white noises which may be regarded, in a crude sense, as analog and digital noises, respectively, with the rather different character.

One of the controversial issues concerning cortical neurons is how the neurons may communicate information by spikes with a large variability. It has been reported that the variability of $`c_v=0.51.0`$ is observed in spike trains of non-bursting cortical neurons in visual V1 and MT of monkey , which is in strong contrast with a small $`c_v(=0.050.1)`$ in motor neurons . There have been much discussions how to understand the observed large variability: a balance between excitatory and inhibitory inputs , the high physiological gain in the plot of input current vs. output frequency , correlation fluctuations in recurrent networks , the active dendrite conductance , input ISIs with the distribution of a slow-decreasing tail , and input ISIs with large $`c_v`$ . We expect that although there might be several origins responsible to the observed, large variability in ISI, noises may be one of conceivable mechanisms.

The purpose of the present study is to investigate responses of ensemble neurons to spike-train inputs subject to ST and white noises, in order to get some insight to the following issues.

(1) Is the effect of ST noise on responses to spike-train signals, particalrly on SR, different from or same as that of white noise?

(2) Is a population or ensemble of neurons important for the fidelity of signal transmission?

(3) Is a large variability of spikes compatible with information transmission with a fairly good SNR?

Responses of single HH neurons to various types of spike-train inputs with deterministic, chaotic and stochastic ISIs without noises, have been investigated . SR for coherent spike-train inputs has been theoretically studied with single IF neurons -. We should note, however, that the response to applied, external stimulus of the realistic HH model  is rather different even qualitatively from that of the IF neuron . For an excitatory dc input $`I`$, the IF neuron which is classified as the type I, evokes the self-excited oscillation showing the continuous $`f_oI`$ relation with a wide range of frequency $`f_o`$. On the other hand, the HH neuron, which is classified as the type II, has the discontinuous $`f_oI`$ relation at the critical current, above which it shows the oscillation with a fairy narrow range of $`f_o`$. For an inhibitory dc input current, the HH neuron can fire with the so-called rebound process while the IF neuron cannot. Since the threshold-crossing behavior of the neuron is important in determining the behavior of its SR, it is necessary to re-examine SR for the spike-train inputs with the use of the realistic HH neuron model. Furthermore, since SR of single neurons is generally different from that of ensemble neurons, it is also necessary to investigate SR not only of single HH neurons but also of ensemble HH neurons.

The present paper is organized as follows. In Sec. II, an adopted model for an ensemble of $`N`$-unit ($`N`$ = 1, 10, 100) HH neurons is described. Simulations for responses of ensemble neurons to ST signals with added ST noise, white noise, and ST plus white noises are reported in Sec. IIIA, IIIB and IIIC, respectively, where SR and the variability of ISIs are discussed. The final Sec. IV is devoted to conclusion and discussion.

## II Ensemble Neuron Model

We assume a network consisting of $`N`$-unit HH neurons which receive the same ST signals but independent, ST and Gaussian noises through excitatory synapses. Spikes emitted by the ensemble neurons are collected by a summing neuron. A similar model was previously adopted by several authors studying SR  . Dynamics of the membrane potential $`V_i`$ of the HH neuron i is described by the non-linear differential equations given by

$$\overline{C}dV_i(t)/dt=I_i^{\mathrm{ion}}+I_i^{\mathrm{ps}}+I_i^\mathrm{n},\text{(for }1iN\text{)}\text{}$$
(1)
where $`\overline{C}=1\mu \mathrm{F}/\mathrm{cm}^2`$ is the capacity of the membrane. The first term $`I_i^{\mathrm{ion}}`$ of Eq.(1) denotes the ion current given by

$$I_i^{\mathrm{ion}}=g_{\mathrm{Na}}m_i^3h_i(V_iV_{\mathrm{Na}})+g_\mathrm{K}n_i^4(V_iV_\mathrm{K})+g_\mathrm{L}(V_iV_\mathrm{L}),$$
(2)
where the maximum values of conductivities of Na and K channels and leakage are $`g_{\mathrm{Na}}=120\mathrm{mS}/\mathrm{cm}^2`$, $`g_\mathrm{K}=36\mathrm{mS}/\mathrm{cm}^2`$ and $`g_\mathrm{L}=0.3\mathrm{mS}/\mathrm{cm}^2`$, respectively; the respective reversal potentials are $`V_{\mathrm{Na}}=50`$ mV, $`V_\mathrm{K}=77`$ mV and $`V_\mathrm{L}=54.5`$ mV. Dynamics of the gating variables of Na and K channels, $`m_i,h_i`$ and $`n_i`$, are described by the ordinary differential equations, whose details have been given elsewhere .

The second term $`I_i^{\mathrm{ps}}`$ in Eq.(1) denotes the post-synaptic current given by

$$I_i^{\mathrm{ps}}=\underset{m}{}g_s(V_aV_s)\alpha (tt_{\mathrm{i}m}),$$
(3)
which is induced by an input spike with the magnitude $`V_a`$ given by

$$U_i(t)=V_a\underset{m}{}\delta (tt_{\mathrm{i}m}),$$
(4)
with the alpha function $`\alpha (t)`$:

$$\alpha (t)=(t/\tau _\mathrm{s})e^{t/\tau _\mathrm{s}}\mathrm{\Theta }(t).$$
(5)
In Eqs.(3)-(5) $`t_{\mathrm{i}m}=(m1)T_\mathrm{s}`$ is the $`m`$-th firing time with the input-signal ISI of $`T_\mathrm{s}`$, the Heaviside function is defined by $`\mathrm{\Theta }(t)=1`$ for $`x0`$ and 0 for $`x<0`$, and $`g_s`$, $`V_\mathrm{s}`$ and $`\tau _\mathrm{s}`$ stand for the conductance, reversal potential and time constant, respectively, of the synapse.

The third term $`I_i^\mathrm{n}`$ in Eq.(1) denotes added, independent noises which consist of two terms:

$$I_i^n(t)=\underset{m}{}\sqrt{C}\alpha (tt_{im}^n)+\sqrt{2D}\xi _i(t).$$
(6)
The first term of Eq.(6) expresses Poisson ST noise, whose magnitude, $`C`$, is hereafter expressed by $`\sqrt{C}g_n(V_aV_s)`$ in terms of $`g_n`$ as in Eq.(3) for a later purpose, and $`t_{im}^n`$ is the $`m`$-th firing time of the ST noise of the neuron $`i`$ with the average ISI of $`\mu _n`$. The second term of Eq.(6) denotes Gaussian white noises with the magnitude of $`D`$ given by

$$<\overline{\xi _j(t)}>=0,$$
(7)
$$<\overline{\xi _j(t)\xi _k(t^{})}>=\delta _{jk}\delta (tt^{}),$$
(8)
where the overline $`\overline{X}`$ and the bracket $`<X>`$ denote the temporal and spatial averages, respectively.

We should remark that our ensemble neuron model given by Eqs.(1)-(6) does not include synaptic couplings among constituent HH neurons, in contrast with the coupled ensemble models ; related discussions being given in Sec. IV.

We assume that information is carried by firing times of spikes. Dividing the time scale by the width of time bin of $`T_\mathrm{b}`$ as $`t=t_{\mathrm{}}=(\mathrm{}1)T_\mathrm{b}`$ ($`\mathrm{}`$: integer), we define input and output signals summed over ensemble neurons within the each time bin by

$$W_\mathrm{i}(t)=\underset{m}{}\mathrm{\Theta }(T_\mathrm{b}/2tt_{\mathrm{i}m}),$$
(9)
$$W_\mathrm{o}(t)=(1/N)\underset{i=1}{\overset{N}{}}\underset{n}{}\mathrm{\Theta }(T_\mathrm{b}/2tt_{\mathrm{o}in}).$$
(10)
In Eqs. (9) and (10) $`\mathrm{\Theta }(t)`$ stands for the Heaviside function, $`t_{\mathrm{i}m}`$ the $`m`$-th firing time of inputs, and $`t_{\mathrm{o}in}`$ the $`n`$-th firing time of outputs of the neuron $`i`$ when $`V_i(t)`$ crosses $`V_z=0`$ mV from below. The time bin is chosen as $`T_\mathrm{b}`$ =2.5 ms in our simulations. The fast Fourier transformation (FT) is performed for $`W_\mathrm{o}(t)`$ in order to get the SNR defined by

$$SNR=10\mathrm{log}_{10}(A_s/A_n)\text{(dB)},$$
(11)
where $`A_s`$ is the signal power spectrum at a given frequency of $`1/T_\mathrm{s}`$ and $`A_n`$ the background noise level.

Differential equations given by Eqs.(1)-(6) are solved by the forth-order Runge-Kutta method by the integration time step of 0.01 ms with double precision. The initial conditions for the variables are given by $`V_i(t)`$= -65 mV, $`m_i(t)`$=0.0526, $`h_i(t)`$=0.600, $`n_i(t)`$=0.313 at $`t=0`$, which are the rest-state solution of a single HH neuron. Hereafter time, voltage, conductance, current, and $`D`$ are expressed in units of ms, mV, $`\mathrm{mS}/\mathrm{cm}^2`$, $`\mu \mathrm{A}/\mathrm{cm}^2`$, and $`\mu \mathrm{A}^2/cm^4`$, respectively. We have adopted parameters of $`V_a=30`$, $`V_c=50`$, and $`\tau _\mathrm{s}`$ =2. Adopted values of $`g_s`$, $`g_n`$, $`D`$, $`\mu _\mathrm{n}`$ and $`N`$ will be described shortly. The simulation time for each run for a given set of parameters is $`T_{sim}`$ 1500 ms (150000 $`\times N`$ time steps) and initial 3000 $`N`$ time steps are discarded to get asymptotic solutions. The size of sample data for FT analysis becomes $`N_{FT}`$=512 when the input ISI is chosen to be $`T_\mathrm{s}`$ = 25 ms. A single simulation with $`N=100`$ requires the CPU time of about 150 minutes by DOS/V PC with 900 MHz processor. The calculated SNR is expected to be improved if the simulation time $`T_{sim}`$ is increased. Unfortunately, we have not been able to adopt larger $`T_{sim}`$ because of a limitation in our computer facility (related discussion will be presented in Sec. IV). We hope that our simulation reveals some of the interesting features of SR of ensemble neurons for spike-train inputs with added, multiple noises.

## III CALCULATED RESULTS

### A Spike-Train Noises

#### 1 SNR and SR

Firstly we discuss the case in which ensemble HH neurons receive input signals with independent ST noise only. The input ISI is assume to be $`T_\mathrm{s}`$ = 25 ms because spikes with this value of ISI are reported to be ubiquitous in cortical brains . We study sub-threshold inputs with $`g_s<g_{th}`$ for which neurons cannot fire without added noises. The threshold value of $`g_{th}`$ generally depends on the input ISI, $`T_s`$: for example, $`g_{th}`$=0.088, 0.085 and 0.095 for $`T_s`$=20, 25 and 30 ms, respectively .

Raster in Fig. 1 shows firings of $`N=100`$ ensemble neurons when spike-train signals of $`T_s=25`$ ms are applied to ensemble neurons with added ST noises of $`g_n=0.10`$, $`D=0.0`$ and $`\mu _n=25`$ ms. Neurons fire when a spike-train input plus noise or noises exceed the threshold level. At a glance, firings in Fig. 1 seem random. We may realize it is not true when firings are summed over ensemble neurons as shown in Fig. 2(b) where the output signal $`W_o(t)`$ summed over a $`N=100`$ ensemble is plotted. We note that $`W_o(t)`$ includes a periodic component with a period of $`T_s`$=25 ms as an input signal shown in Fig. 2(a). This regularity is more clearly seen in its FT power spectrum shown in Fig. 3(a), which shows clear peaks at a fundamental frequency of $`1/T_s`$ = 40 Hz and its harmonics. Information may be transmitted with aid of SR in a $`N=100`$ neuron ensemble.

When the size of ensemble neurons is small, however, the information transmission is much degraded. Figures 2(c) and 2(d) show $`W_o(t)`$ for $`N=10`$ and $`N=1`$, respectively, other parameters except $`N`$ being the same as in Fig. 2(b). Figure 2(d) for $`N=1`$ , for example, shows that a single neuron intermittently fires, yielding 18 firings for 40 input spike inputs. The FT power spectra for $`N=10`$ and $`N=1`$ are expressed in Fig. 3(b) and 3(c), respectively. Figures 3(a)-3(c) show that the magnitude of a peak at the fundamental frequency divided by that of background noises, i.e. SNR, increases as the size of ensemble neurons is increased. These mean that despite of ostensibly irregular firings of single neurons, output signals summed over an ensemble, $`W_o(t)`$, become more regular when $`N`$ becomes sufficiently large.

Figure 4 shows SNR as a function of $`g_n`$ calculated for $`g_s=0.06`$, $`D=0`$ and $`\mu _n=25`$ ms with $`N`$ = 1, 10 and 100. Results for $`N`$= 1and 10 are averages of ten runs. When the noise intensity $`g_n`$ is increased from zero, neurons begin to occasionally fire by a cooperative action between input signals and added noises. As $`g_n`$ becomes greater than the threshold level of $`g_{th}=0.085`$, ST noise alone is sufficient to trigger firings without signal inputs. When $`g_n`$ is further increased, SNR of outputs is gradually degraded. Then SNR has a maximum at $`g_ng_{th}`$, which is the characteristics of SR. For the case of $`N=10`$ and 100, we realize SR with the maximum in SNR for a weak noise of $`g_n0.07`$. On the contrary, for the case of $`N=1`$, the maximum in SNR is not evident although SNR is enhanced at $`g_n>0.04`$.

In many non-linear systems, SNR for sinusoidal input signals is reported to obey the noise-intensity dependence given by

$$SNR=10\mathrm{log}_{10}[(A/X^\alpha )\mathrm{exp}(B/X)],$$
(12)
with a maximum at $`X_{max}=B/\alpha `$ where $`X=D`$ (the intensity of the white noise), $`\alpha =2`$, and A and B are constants depending on the model parameters. Substituting $`X=g_n^2`$ \[$`C`$ in Eq.(6)\], we have tried to analyze the $`g_n`$ dependence of SNR for $`N=100`$ obtained in our simulation. Dashed curves in Fig.4 show the results adopting two sets of parameters in Eq.(12): ($`\alpha ,A,B`$)=(2, $`1.77\times 10^2`$, $`8.45\times 10^3`$) and (3, $`2.00\times 10^4`$, $`1.27\times 10^2`$), which are chosen such as to locate the maximum of SNR at $`g_n=0.065`$. The latter choice of parameters with $`\alpha =3`$ yields the much better agreement with the data obtained in our simulations than the former with $`\alpha =2`$, although the plateau around the maximum is not well reproduced even in the latter.

We have investigated the effect of input ISI on SR by changing $`T_s`$. Solid curves in Fig. 5 denote SNRs as a function of $`g_n`$ for $`T_s`$=20, 25 and 30 ms obtained in our simulation and dashed curves those analyzed by Eq.(12) with $`X=g_n^2`$ and sets of parameters of ($`\alpha ,A,B`$) =(3, $`5.01\times 10^4`$, $`1.27\times 10^2`$) for $`T_s=20`$, (3, $`2.00\times 10^4`$, $`1.27\times 10^2`$) for $`T_s=25`$, and (3, $`7.94\times 10^5`$, $`1.27\times 10^2`$) for $`T_s=30`$. The index $`\alpha =3`$ is realized for all the $`T_s`$ values investigated. The maximum SNR value is decreased as $`T_s`$ is increased, although the maximum position at $`g_n0.065`$ in SNR is not changed.

So far our discussion is confined to the sub-threshold case. We have performed simulations also for the supra-threshold inputs, adopting $`g_s`$ larger than the threshold value of $`g_{th}`$. Figure 6 shows the $`g_s`$ dependence of SNR calculated for ST noise only ($`g_n=0.10`$ and $`D=0`$) with $`\mu _n=25`$ ms, and $`N=1`$, 10 and 100. When $`g_s`$ is increased across the threshold value of $`g_{th}`$, SNR is discontinuously increased. SNR for the supra-threshold inputs is better than that for sub-threshold inputs, as expected. We note, however, that ensemble neurons with large $`N`$ is fairly robust against weak noises relevant to SR.

#### 2 ISI

Next we discuss the distribution of output ISIs of individual neurons given by

$$T_{oin}=t_{oin+1}t_{oin},$$
(13)
where $`t_{oin}`$ is the $`n`$-th firing time of outputs of the neuron $`i`$. Figures 7(a)-(d) show histograms of output ISIs for $`g_s`$=0.02, 0.06, 0.10 and 0.14, respectively, when ST noise of $`g_n=0.10`$ is added to ST signal (see Fig. 6). In the cases of $`g_s=0.02`$, ISI histograms (ISIH) nearly obey the exponential distribution as shown by dashed curve, but it vanishes at $`T_{oin}<15`$ ms because a HH neuron cannot properly respond to small-ISI inputs due to its refractory period . On the contrary, in the cases of $`g_s=0.10`$ and 0.14, ISIH has larger magnitudes at $`T_{oin}T_s`$, which may be approximately expressed by the gamma-type distribution. In the case of $`g_s=0.06`$, ISIH includes not only the truncated exponential distribution but also finite contributions at multiples of $`T_s`$. This change in the distribution is more clearly seen in the $`g_s`$ dependence of the average ($`\mu _o`$) and RMS values ($`\sigma _o`$) of output ISIs, which are plotted in Fig. 8. Both $`\mu _o`$ and $`\sigma _o`$ are almost constant at $`g_s<0.08`$, and at $`g_s\stackrel{>}{}0.10`$ they are suddenly decreased. When ISIH of output spikes for $`g_s<g_{th}`$ is expressed by the truncated exponential distribution given by

$$P(T)\mathrm{\Theta }(TT_L)\mathrm{exp}(sT),$$
(14)
$`\mathrm{\Theta }()`$ being the Heaviside function and $`T_L`$ the lower bound, the average and RMS values are given by $`\mu _o=(1+sT_L)/s`$ and $`\sigma _o=1/s`$, which yield $`\mu _o=40`$ and $`\sigma =25`$ ms for $`1/s=25`$ and $`T_L=15`$ ms. These figures are a little different from $`\mu _o50`$ and $`\sigma 35`$ ms for $`g_s<0.8`$ shown in Fig. 8. The difference may be attributed to the extra contribution at $`2T_s`$ obtained in our simulation, as shown in Figs. 7(a) and (b). The variability defined by $`c_{vo}=\sigma _o/\mu _o`$ is 0.655, 0.703, 0.307 and 0.212 for $`g_s`$=0.02, 0.06, 0.10 and 0.14, respectively. When comparing Fig. 6 and Fig. 8, we note that even if the variability of ISIs of individual neurons is considerable, output signals summed over a large-scale ensemble \[$`W_o(t)`$ in Eq.(10)\] may carry information with a fairly good SNR. For example, in the case of $`g_s=0.06`$, we get SNR=17.4 dB for $`c_{vo}`$ = 0.703. This is an advantage of a neuron ensemble.

### B White Noises

#### 1 SNR and SR

Next we discuss SR for Gaussian white noise, which is applied to our ensemble neurons instead of ST noise. Neurons occasionally fire when signal plus noise exceed the threshold level. Firings of neurons for white noises are similar to those to the ST noises shown in Fig. 1. SNR calculated as a function of $`D`$, the intensity of the white noise, for $`N=1`$, 10 and 100 is plotted in Fig. 9. When the white noise intensity is increased from zero, SNR is rapidly enhanced with a maximum at $`D2`$ followed by a gradual decrease, which is a typical SR. Although the calculated SNR shows the SR behavior irrespective of N, it is more evident for larger $`N`$.

The dashed curve in Fig. 9 expresses SNR calculated by Eq.(12) with a set of parameters of ($`\alpha ,A,B`$)=(2, 5890, 5), which are chosen such as to agree with the maximum position at $`D=2.5`$ in SNR for $`N=100`$ obtained by our simulation. The agreement between the result of $`\alpha =2`$ and our data seems satisfactory. This value of $`\alpha =2`$ agrees with the results of SR for sinusoidal input signals in HH neurons  as well as those realized in many non-linear systems .

Fig. 10 shows the $`g_s`$ dependence of SNR for white noises only ($`D=2`$ and $`g_n=0`$) with $`N=1`$, 10 and 100. SNR is gradually increased as increasing $`g_s`$. In contrast with the case for ST noise shown in Fig. 6, there is no significant changes in SNR at the threshold level shown by the vertical, dashed line. A comparison of Fig. 6 and Fig. 10 shows that $`g_s`$ dependence of SNR for ST noise is different from that for white noise.

#### 2 ISI

Figures 11(a)-(d) show histograms of output ISIs for $`g_s`$=0.02, 0.06, 0.10 and 0.14, respectively, when white noise of $`D=2`$ is added to ST signal (see Fig. 10). ISIHs for $`g_s=0.02`$ and $`g_s=0.06`$ show a typical behavior with peaks at multiples of $`T_s`$ whose magnitudes decrease exponentially . It is noted that the distribution extends up to 200 ms. As $`g_s`$ is increased across $`g_{th}`$, magnitudes of the main peak at $`T_s`$ are much increased as expected. Figure 12 shows the $`g_s`$ dependence of the average ($`\mu _o`$) and RMS values ($`\sigma _o`$) of ISI and the variability ($`c_{vo}`$). As increasing $`g_s`$, both $`\mu _o`$ and $`\sigma _o`$ are gradually decreased with no sudden changes at $`g_sg_{th}`$. The variability is 0.813 and 0.710 for $`g_s`$=0.02 and 0.06, respectively. As $`g_s`$ becomes larger than $`g_{th}`$, ISIH has a larger magnitude at $`T_{oin}T_s`$ and then the variability is decreased: $`c_{vo}`$ becomes 0.528 and 0.343 for $`g_s`$=0.10 and 0.14, respectively. Comparing Fig. 12 to Fig. 8, we note that $`g_s`$ dependence of $`\mu _o`$, $`\sigma _o`$ and $`c_{vo}`$ for white noises is rather different from that for ST noise.

### C Spike-Train plus White Noises

Since neurons in real neural systems are in the environment with various kinds of noises as discussed in Sec. I, it is necessary to examine various effects of multiple types of noises. Taking into account independent Gaussian and OU noises, Liu, Hu and Wang have investigated the effect of spatial correlation on SR in coupled HH neurons . Recently Lindner and Schimansky-Geier have included, in the ensemble IF model, the additive and signal-coded noises, which are expressed by $`\sqrt{2D_1}\xi _1(t)`$ and $`\sqrt{D_2s(t)}\xi _2(t)`$, respectively, in terms of the external sinusoidal signal $`s(t)`$ and Gaussian noises $`\xi _n`$ with the magnitudes of $`D_n`$ ($`n=1,2`$).

In previous Sec. IIIA and Sec. IIIB, we have separately discussed ST and white noises. Now we simultaneously add both the noises to our ensemble neurons. Figure 13 shows the three-dimensional plot of SNR for $`N=100`$ as functions of $`g_n`$ and $`D`$. In the case of $`D=0`$ (ST noise only), SNR has a maximum at $`g_n0.07`$, as shown in Fig.4. In the case of $`g_n=0`$ (white noises only), on the other hand, SNR has a maximum at $`D2`$ as shown in Fig. 9. The contour plot depicted in the base of Fig. 13 shows that SNR is rapidly increased from zero as ST or white noise is increased. SNR in the presence of weak white noise of $`0<D<1`$ is enhanced by a further addition of ST noise and it depends considerably on $`g_n`$. With stronger white noise of $`D>1`$, SNR is slightly enhanced by an addition of weak ST noise of $`g_n0.01`$ although it only weakly depends on $`g_n`$ at $`g_n>0.02`$. SNR in the presence of weak ST noise with $`g_n<0.05`$ is much enhanced by a further addition of white noise. We note, however, that white noise enhances SNR even in the presence of stronger ST noise of $`g_n>0.1`$, where SNR is decreased by excess ST noise for $`D=0`$. These clearly show that white noise is more effective for SR than ST noise.

## IV CONCLUSION AND DISCUSSION

In previous Sec. III, the simulation time of each run has been limited to be 1500 ms (150000 $`\times N`$ time steps) because of a limitation in our computer facility. It is, however, possible, to extend the simulation time when the simulation is made only for single ($`N=1`$) neurons. Figure 14 shows SNR as a function of $`g_n`$ for ST noise added to single HH neurons, which are calculated with the simulations times of $`T_{sim}=`$ 1500, 3000, 6000 and 12 000 ms, yielding FT-data sizes ($`N_{FT}`$) of 512, 1024, 2048 and 4096, respectively. We note that as $`T_{sim}`$ becomes larger, SNR is improved, in particular, its maximum becomes more evident. The dashed curve expresses SNR calculated by using Eq.(12) with $`X=g_n^2`$ and a set of parameters of $`(\alpha ,A,B)`$ =(3, $`6.31\times 10^5`$, $`9.75\times 10^3`$), which are chosen to reproduce SNR calculated for $`T_{sim}=12000`$ ms with the maximum at $`g_n=0.057`$. Again the index of $`\alpha =3`$ is realized for ST noise added to single neurons. We note from Figs. 4 and 14 that SR becomes more significant if the duration of applied, coherent spike trains is longer and/or the size of ensemble neurons is larger. Even if the duration of applied signal is not long, SNR may be improved if the size of ensemble neurons is sufficiently large. This is more evident when the input signal is transient spike train, as recently demonstrated in Ref..

In a summary, we have numerically investigated SR responses of an ensemble of HH neurons to spike-train signals with added ST and/or white noises. Our conclusions against the three issues raised in the Introduction are summarized as follows:

(1) Comparisons of Figs. 4, 6, 7 and 8 for ST noise with Figs. 9, 10, 11 and 12 for white noise, respectively, clearly show both the difference and the similarity between the effects of ST and white noises. Although SR is a genetic phenomenon, its detailed behavior depends on kinds of the input signal and added noises. When analyzing SNR obtained in our simulations with the use of Eq.(12), we get $`\alpha =3`$ for ST noise, which is different from $`\alpha =2`$ for white noise. ST noise is less effective for SR than the white noise (Fig. 13).

(2) SNR is more improved as the size of ensemble is larger. In a large neuron ensemble, the transmission fidelity for supra-threshold inputs is not significantly degraded by weak noises responsible to SR for sub-threshold inputs (Figs. 6 and 10).

(3) The variability of ISIs of individual neurons for sub-threshold inputs is rather large ($`c_{vo}\stackrel{<}{}0.8`$). Nevertheless the output $`W_o`$ summed over an ensemble may carry information with a fairly good SNR.

The item (2) is consistent with the results of SR for transient spike-train signal  and for analog signals , showing that a population of neurons plays a very important role for the transmission of spike-train inputs both with sub- and supra-threshold levels. It is worth to note that this enhancement in SNR is due to the pooling effect  because our ensemble neuron model have no couplings among HH neurons. The item (3) shows that noise may be one of conceivable mechanism yielding a large variability observed experimentally . Even when the variability of firings of individual neurons is considerable, firings summed over a ensemble may carry information with a fairly good SNR enhanced by SR and pooling effects. Thus the large variability and high SNR are not incompatible in a large-scale neuron ensemble.

The present study entirely relies on simulations. It would be interesting to theoretically elucidate the dependence of the index $`\alpha `$ on a kind of added noises mentioned above. However, conventional approaches having been employed for a study of SR such as the rate-equation and linear-response theories -, do not work on our case. We leave its analytical study as our future problem.

## Acknowledgements

This work is partly supported by a Grant-in-Aid for Scientific Research from the Japanese Ministry of Education, Culture, Sports, Science and Technology.