text:
# I Introduction

## I Introduction

The relativistic mean field (RMF) theory  has been very successful in describing the ground state (zero temperature) properties of nuclei over the entire periodic table in recent years. The binding energies, charge radii and the ground state deformations are reproduced very well; the charge distributions also compare extremely well with the experimental data. This theory has proved to be very fruitful in explaining  various details of exotic nuclei near the drip lines. In contrast to the nonrelativistic models, the RMF theory uses a single set of parameters to explain all these properties. To our knowledge, this approach has not yet been exploited to understand the properties of hot nuclei except some preliminary investigations for closed-shell nuclei . The response of nuclear shapes to thermal excitations, for example, has been experimentally studied from the shapes of the giant dipole resonances (GDR) built on excited states . Theoretically such finite temperature effects have been studied till date in the nonrelativistic framework such as the Hartree-Fock-Bogoliubov (HFB) theory  and the Landau theory of phase transition . Such theories qualitatively explain, for example, the temperature evolution of nuclear shapes. Recent experiments indicate, however, that a quantitative estimate of the persistence of the ground state deformation  with temperature may be missing in some cases.

Against this backdrop, we undertake the study of the thermal response to nuclear properties, particularly the deformation and the level density in the RMF framework. The HFB calculations employ a model hamiltonian in a limited model space with pairing-plus-quadrupole interaction. This may not be very realistic once temperature comes into play. Moreover, to make the calculations numerically tractable, an inert core is assumed. This may be questionable at moderately high temperatures. The RMF approach we use is effectively free from these limitations. The model space employed is sufficiently large and all the nucleons are treated on equal footing.

In the present work, we employ the nonlinear $`\sigma \omega \rho `$ version of the RMF theory . In absence of a simple relativistic prescription for the pairing interaction, it is introduced phenomenologically. We take two representative systems, namely, $`{}_{}{}^{166}Er`$ and $`{}_{}{}^{170}Er`$. We investigate the thermal evolution of the nuclear shapes and the pairing gaps. The temperature dependence of the specific heat as a possible signature of phase transition in pairing and nuclear shapes is explored. The temperature variation of the nuclear level density parameter which has a pre-eminent role in understanding nuclear reactions is further studied. Effects of thermal fluctuations of the nuclear shapes on the deformation and the nature of the phase transition is also discussed.

The theoretical framework used is briefly discussed in Sec.II. The results and discussions are presented in Sec.III and the Sec.IV contains the concluding remarks.

## II Formalism

The Lagrangian density for the nucleon-meson many-body system  is taken as

$``$ $`=`$ $`\overline{\mathrm{\Psi }}_i\left(i\gamma ^\mu _\mu M\right)\mathrm{\Psi }_i+{\displaystyle \frac{1}{2}}^\mu \sigma _\mu \sigma U(\sigma )g_\sigma \overline{\mathrm{\Psi }}_i\sigma \mathrm{\Psi }_i`$ (3)
$`{\displaystyle \frac{1}{4}}\mathrm{\Omega }^{\mu \nu }\mathrm{\Omega }_{\mu \nu }+{\displaystyle \frac{1}{2}}m_\omega ^2\omega ^\mu \omega _\mu g_\omega \overline{\mathrm{\Psi }}_i\gamma ^\mu \omega _\mu \mathrm{\Psi }_i{\displaystyle \frac{1}{4}}\stackrel{}{R}^{\mu \nu }\stackrel{}{R}_{\mu \nu }+{\displaystyle \frac{1}{2}}m_\rho ^2\stackrel{}{\rho }^\mu \stackrel{}{\rho }_\mu `$
$`g_\rho \overline{\mathrm{\Psi }}_i\gamma ^\mu \stackrel{}{\rho }_\mu \stackrel{}{\tau }\mathrm{\Psi }_i{\displaystyle \frac{1}{4}}F^{\mu \nu }F_{\mu \nu }e\overline{\mathrm{\Psi }}_i\gamma ^\mu {\displaystyle \frac{(1\tau _3)}{2}}A_\mu \mathrm{\Psi }_i.`$
The meson fields included are the isoscalar $`\sigma `$ meson, the isoscalar-vector $`\omega `$ meson and the isovector-vector $`\rho `$ meson. The arrows in eq. (3) denote the isovector quantities. The Largrangian contains a nonlinear scalar self-interaction term $`U(\sigma )`$ of the $`\sigma `$ meson:

$$U(\sigma )=\frac{1}{2}m_\sigma ^2\sigma ^2+\frac{1}{3}g_2\sigma ^3+\frac{1}{4}g_3\sigma ^4.$$
(4)
This term is important for appropriate description of the surface properties . The quantities $`M`$, $`m_\sigma `$, $`m_\omega `$ and $`m_\rho `$ are the nucleon, $`\sigma `$, $`\omega `$ and the $`\rho `$meson masses, respectively, while $`g_\sigma `$, $`g_\omega `$, $`g_\rho `$ and $`e^2/4\pi =1/137`$ are the corresponding coupling constants for the mesons and the photon. The field tensors of the vector mesons and of the electromagnetic fields have the following structure:

$`\mathrm{\Omega }^{\mu \nu }`$ $`=`$ $`^\mu \omega ^\nu ^\nu \omega ^\mu ,`$ (5)
$`\stackrel{}{𝐑}^{\mu \nu }`$ $`=`$ $`^\mu \stackrel{}{\rho }^\nu ^\nu \stackrel{}{\rho }^\mu g_\rho (\stackrel{}{\rho }^\mu \times \stackrel{}{\rho }^\nu ),`$ (6)
$`F^{\mu \nu }`$ $`=`$ $`^\mu A^\nu ^\nu A^\mu .`$ (7)
The variational principle gives the equations of motion. The mean field approximation is introduced at this stage by treating the fields as $`c`$numbers or classical fields. This results in a set of coupled equations, namely the Dirac equation with potential terms for the nucleons and the Klein-Gordon type equations with sources for the mesons and the photon. For the static case, along with the time reversal invariance and charge conservation the equations get simplified. The resulting equations, known as RMF equations, have the following form. The Dirac equation for the nucleon is

$$\{i\alpha +V(𝐫)+\beta \left[M+S(𝐫)\right]\}\mathrm{\Psi }_i=ϵ_i\mathrm{\Psi }_i,$$
(8)
where $`V(𝐫)`$ represents the vector  potential

$$V(𝐫)=g_\omega \omega _0(𝐫)+g_\rho \tau _3\rho _0(𝐫)+e\frac{(1\tau _3)}{2}A_0(𝐫),$$
(9)
and $`S(𝐫)`$ is the scalar potential

$$S(𝐫)=g_\sigma \sigma (𝐫),$$
(10)
which contributes to the effective mass as

$$M^{}(𝐫)=M+S(𝐫).$$
(11)
The Klein-Gordon equations for the mesons and the electromagnetic fields with the nucleon densities as sources are

$`\left\{\mathrm{\Delta }+m_\sigma ^2\right\}\sigma (𝐫)`$ $`=`$ $`g_\sigma \rho _s(𝐫)g_2\sigma ^2(𝐫)g_3\sigma ^3(𝐫),`$ (12)
$`\left\{\mathrm{\Delta }+m_\omega ^2\right\}\omega _0(𝐫)`$ $`=`$ $`g_\omega \rho _v(𝐫),`$ (13)
$`\left\{\mathrm{\Delta }+m_\rho ^2\right\}\rho _0(𝐫)`$ $`=`$ $`g_\rho \rho _3(𝐫),`$ (14)
$`\mathrm{\Delta }A_0(𝐫)`$ $`=`$ $`e\rho _c(𝐫).`$ (15)
The corresponding densities are

$`\rho _s`$ $`=`$ $`{\displaystyle \underset{i}{}}n_i\overline{\mathrm{\Psi }}_i\mathrm{\Psi }_i,`$ (16)
$`\rho _v`$ $`=`$ $`{\displaystyle \underset{i}{}}n_i\mathrm{\Psi }_i^{}\mathrm{\Psi }_i,`$ (17)
$`\rho _3`$ $`=`$ $`{\displaystyle \underset{i}{}}n_i\mathrm{\Psi }_i^{}\tau _3\mathrm{\Psi }_i,`$ (18)
$`\rho _c`$ $`=`$ $`{\displaystyle \underset{i}{}}n_i\mathrm{\Psi }_i^{}{\displaystyle \frac{(1\tau _3)}{2}}\mathrm{\Psi }_i.`$ (19)
Here the sums are taken over the particle states only, i.e., the negative-energy states are neglected. The partial occupancies ($`n_i`$) at finite temperature in the constant pairing gap approximation (BCS) is

$$n_i=\frac{1}{2}\left[1\frac{ϵ_i\lambda }{\stackrel{~}{ϵ}_i}\left(12f(\stackrel{~}{ϵ}_i,T)\right)\right],$$
(20)
with $`f(\stackrel{~}{ϵ}_i,T)=1/(1+e^{\stackrel{~}{ϵ}_i/T})`$; $`\stackrel{~}{ϵ}_i=\sqrt{(ϵ_i\lambda )^2+\mathrm{\Delta }^2}`$ is the qusiparticle energy where $`ϵ_i`$ is the single-particle energy for the state $`i`$. The chemical potential $`\lambda `$ for protons (neutrons) is obtained from the requirement

$$\underset{i}{}n_i=Z(N)$$
(21)
The sum is taken over proton (neutron) states. The gap parameter $`\mathrm{\Delta }`$ is obtained by minimising the free energy

$$F=ETS,$$
(22)
where

$$E(T)=\underset{i}{}ϵ_in_i+E_\sigma +E_{\sigma NL}+E_\omega +E_\rho +E_C+E_{pair}+E_{CM}AM,$$
(23)
and

$$S=\underset{i}{}\left[f_ilnf_i+(1f_i)ln(1f_i)\right],$$
(24)
with

$`E_\sigma `$ $`=`$ $`{\displaystyle \frac{1}{2}}g_\sigma {\displaystyle d^3r\rho _s(𝐫)\sigma (𝐫)},`$ (25)
$`E_{\sigma NL}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle d^3r\left\{\frac{1}{3}g_2\sigma ^3(𝐫)+\frac{1}{2}g_3\sigma ^4(𝐫)\right\}},`$ (26)
$`E_\omega `$ $`=`$ $`{\displaystyle \frac{1}{2}}g_\omega {\displaystyle d^3r\rho _v(𝐫)\omega ^0(𝐫)},`$ (27)
$`E_\rho `$ $`=`$ $`{\displaystyle \frac{1}{2}}g_\rho {\displaystyle d^3r\rho _3(𝐫)\rho ^0(𝐫)},`$ (28)
$`E_C`$ $`=`$ $`{\displaystyle \frac{e^2}{8\pi }}{\displaystyle d^3r\rho _C(𝐫)A^0(𝐫)},`$ (29)
$`E_{pair}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }^2}{G}},`$ (30)
$`E_{CM}`$ $`=`$ $`{\displaystyle \frac{3}{4}}\mathrm{}\omega _0={\displaystyle \frac{3}{4}}41A^{1/3}.`$ (31)
Here $`G`$ and $`A`$ are the pairing strength and the mass number respectively. The single-particle energies and the fields appearing in eqs. (23) \- (28) are obtained from the self-consistent solution of eqs. (8) - (14). We generate these self-consistent solutions using the well tested basis expansion method as described in detail in Refs. .

The self-consistent solutions for the RMF equations at a given temperature is obtained by minimizing the free energy $`F`$ which yields the equilibrium or the most probable value of the quadrupole deformation ($`\beta _2^0`$) and the proton (neutron) pairing gaps $`\mathrm{\Delta }_p`$ ($`\mathrm{\Delta }_n`$). To incorporate the effects of the thermal fluctuations, a constrained calculation is performed to generate the free energies for the deformations away from the equilibrium value. Average value of a quantity then can be obtained as

$$\overline{O}=\frac{𝑑\beta _2O(\beta _2)e^{\mathrm{\Delta }F(\beta _2)/T}}{𝑑\beta _2e^{\mathrm{\Delta }F(\beta _2)/T}}$$
(32)
where $`\mathrm{\Delta }F=F(\beta _2)F(\beta _2^0)`$ and $`O(\beta _2)`$ is the expectation value of the operator $`\widehat{O}`$ at a fixed value of $`\beta _2`$ and $`T`$. The quantity $`e^{\mathrm{\Delta }F/T}`$ is a measure of the probablity for the nucleus to have deformation $`\beta _2`$ at the temperature $`T`$.

## III Results and Discussions

We have chosen $`{}_{}{}^{166}Er`$ and $`{}_{}{}^{170}Er`$ as two representative systems for our calculations. The results presented are obtained using the NLSH parameter set  for the values of the coupling constants and masses for the mesons and nucleons. The pairing strength $`G`$ is taken to be $`29/A`$ and $`21/A`$ for protons and neutrons, respectively; they reasonably reproduce the observed pairing gaps at zero temperature obtained from odd-even mass differences. The values of the chemical potential and the pairing gap are determined using all the single particle states up to $`2\mathrm{}\omega _0`$ above the Fermi surface without assuming any core, i.e., the mean field solution is generated by taking into account all the nucleons in the systems considered.

At finite temperatures, the nucleus is not strictly in thermal equilibrium. In order to treat the system as an isolated one in equilibrium, one has to take into account corrections due to excitations of nucleons in the continuum. This is generally done through a subtraction procedure, by treating the liquid plus vapour phase together and then the vapour phase seperately . However, calculations of the nuclear level density parameter at finite temperature show that the results are insensitive  to the continuum corrections for $`T`$ upto $``$ 3 MeV; we have therefore not included the effects due to the continuum in the present calculations as temperatures above it are not relevant.

In subsec.3.1, the results for the most probable values (i.e, the mean field values ) of the different observables are presented. With increasing temperature, thermal fluctions build up which may shift the average values from the most probable ones. In subsec.3.2, we discuss the results with the inclusion of thermal fluctuations. Fluctuations in both the $`\beta _2`$ and $`\gamma `$\- degrees of freedom of the nuclear shape should be incorporated, however, because of simplicity and computational economy, we have included only the $`\beta _2`$\- fluctuations in the present calculations.

### A Mean field results

We have calculated the most probable values of the quadrupole deformation parameter, neutron and proton pairing gaps, heat capacity and the level density parameters for the aforesaid systems as a function of temperature.

In Fig. 1, the variation of the quadrupole deformation parameter $`\beta _2`$ as a function of temperature is shown. The ground state deformations obtained here are almost close to the ones calculated in Refs. using nonrelativistic mean field theory with pairing plus quadrupole (P+Q) interaction. The nature of the temperature dependence of the deformation in the two cases are also not very different. However, whereas in the nonrelativistic case the phase transition from prolate to spherical shape occurs at $`T`$1.8 MeV, in the present case the said transition is found at a higher temperature, $`T`$2.7 MeV. The signature of the phase-transition can be inferred by examining other observables like heat capacity which we consider later.

In Fig. 2, the results for the evolution of the neutron and proton pairing gaps with temperature are displayed. For both the nuclei the pairing gaps monotonically decrease with the increase in temperature. It is seen that for both the nuclei, the neutron pairing gaps vanish almost at the same temperature $`T`$0.4 MeV; similarly the proton pairing gaps vanish at $`T`$0.45 MeV. These results are not much in variance with those obtained in the nonrelativistic P+Q model.

The vanishing of the nuclear deformation and the pairing gaps with temperature indicates that there is a shape transition from prolate to spherical and also a transition from the superfluid phase to the normal phase. To understand the nature of the transitions, the heat capacity at various temperatures are calculated. The specific heat $`C`$ for a given temperature is obtained using

$$C(T)=\frac{E^{}}{T}$$
(33)
where $`E^{}`$ is the excitation energy of the nucleus. In Fig. 3, the temperature variation of the specific heat for $`{}_{}{}^{166}Er`$ and $`{}_{}{}^{170}Er`$ nuclei is plotted. At $`T`$ 0.4 MeV, two closely seperated peaks in the $`C(T)`$ curves are seen. They correspond to the dissolution of the neutron and proton pairing gaps. These are the characteristic signatures of second order phase transitions from superfluid to normal phase. It is further noted that the heat capacities for the two systems have broad bumps at $`T2.7`$ MeV. This signifies a weak second order phase transition corresponding to the transition of the nuclear shape. This result is at variance with those obtained in the P+Q model or in the nonrelativistic calculation  based on the Landau theory of phase transitions where a strong second order phase transition is observed.

We now present the results for the temperature dependent level density parameter $`a`$. The parameter $`a`$ can be obtained using the excitation energy and the entropy as follows,

$`E^{}`$ $`=`$ $`aT^2,`$ (34)
$`S`$ $`=`$ $`2aT.`$ (35)
The parameter $`a`$ obtained using above eqs. (34) and (35) would be equal provided it is independent of temperature . In Fig. 4, the inverse level density parameter $`K=A/a`$ ($`A`$ is the mass number of the nucleus) is plotted as a function of temperature. The subscripts $`E`$ and $`S`$ are used to distinguish the two definitions given by eqs. (34) and (35), respectively. At lower temperature, both $`K_E`$ and $`K_S`$ shoot up due to the collapse of the pairing gaps. At higher temperatures, $`K_E`$ and $`K_S`$ are quite close to each other and there is no appreciable variations in their values in the temperature range considered. This is due to the weak transition of the nuclear shapes from deformed to spherical at $`T2.7`$ MeV.

In order to test the sensitivity of the results presented above to the choice of the parameter set and the model space, we have repeated the calculations for the nucleus $`{}_{}{}^{166}Er`$ with the NL3 parameter set . In comparison to the results for NLSH parameter set it is found that the NL3 parameter set gives slightly larger ($`10\%`$) value of the ground state deformation. For $`T>1.5MeV`$, the values of deformation obtained in both NLSH and NL3 parameter sets are almost identical. The temperature evolution of the pairing gaps and level density parameters are practically the same with the two parameter sets. We extend the model space to include single particle states upto $`3\mathrm{}\omega _0`$ above the Fermi surface instead of $`2\mathrm{}\omega _0`$ as used above. For this extended model space the pairing strength $`G`$ is adjusted to reproduce the ground state pairing gap. We do not find significant changes in the values of deformation parameter and the pairing energy in the temperature range of interest ($`T3`$ MeV). In the extended model space, the value of the inverse level density parameter is also nearly the same. The model space used is thus found to be sufficient for our calculations. To estimate the importance of the continuum corrections to the inverse level density parameter we calculated the occupancy, $`n^{(+)}`$, of the single-particle states with positive energy. For $`T<\mathrm{\hspace{0.25em}1}`$ MeV practically there is no particle in the positive energy states ($`n^{(+)}=0`$) and at the highest temperature of interest studied $`n^{(+)}/A=0.018`$ which is very small. Thus it is expected that the continuum corrections may not play an important role in the temperature range we study. The continuum effects may grow stronger for $`T>3`$ MeV, however, this is beyond the transition temperatures and so we have not taken this into account.

The shape transition temperature and the inverse level density parameter obtained in the present model are higher compared to those calculated in the (P+Q) model . Moreover, in the later model, the deformation falls to zero sharply whereas in the present case, this fall is comparatively a little slower. The origin of these differences can be attributed to the single-particle level spectra. In Fig.5, the energy levels for the different proton orbitals around the Fermi surface for the nucleus $`{}_{}{}^{170}Er`$ at $`T`$=3.0 MeV are displayed. Also the proton single-particle energies in the (P+Q) model are shown at this temperature. In both models the orbits are spherical at this temperature. The striking difference between the orbitals in the two models is the large energy gap across the Fermi surface in the RMF model. Similar is the case for neutrons. These tend to reduce the entropy and the excitation energy and hence the level density parameter. The increase in the shape transition temperature (by $``$ 0.9 MeV) in the RMF model as compared to the P+Q model may be traced back to the specific level structure and the gap across the Fermi surface in the two models and their thermal response.

### B Results with thermal fluctuations

Experimental data on giant dipole resonances built on excited states give estimates of the average values of nuclear deformation at finite temperatures. They can be calculated by taking into account the thermal fluctuation effects around the ’most probable’ value. Nuclei being finite systems, one expects that thermal fluctuations would play an important role in the quantitative estimation of various nuclear observables at finite temperatures. In fact, it has been shown in Ref that the experimental data on GDR built on the excited states can be explained reasonably only if the effects from the thermal fluctuations of all the quadrupole degrees of freedom for nuclear shape are included. In the following, we consider the thermal fluctuations only in the $`\beta _2`$ degree of freedom. Inclusion of fluctuations in the $`\gamma `$-degrees of freedom is very involved and computer-intensive and therefore not considered presently. We perform constrained calculations for free energy at fixed deformations $`\beta _2`$. The free energy surface shows multiple minima for temperatures below the shape transition temperature. The average value of $`\beta _2`$ is calculated using eq.(27). In Fig.6, we have displayed the plots of free energy surface $`\mathrm{\Delta }F`$ for $`{}_{}{}^{170}Er`$ at a few temperatures ranging from 2.0 to 3.0 MeV where $`\mathrm{\Delta }F`$ is measured from the lowest minimum (prolate in the present case) at the corresponding temperature. For $`T`$ = 2.0 MeV, the free energy surface has minima at $`\beta _2=0.16`$ and 0.23. The free energy difference between these two minima is 0.70 MeV (i.e. $`\mathrm{\Delta }F/T=0.35`$). As the temperature approaches the shape transition temperature ($`T`$2.7 MeV), both the minima in the free energy surface merge and flatten the bottom part of the free energy surface. This is evident from the middle panel of Fig.6 corresponding to $`T=`$ 2.65 MeV. With further increase in temperature, the free energy surface with single minimum at zero deformation broadens. The variation of the average deformation ($`\overline{\beta }_2`$)as a function of temperature obtained through eq.(32) for $`{}_{}{}^{170}Er`$ is shown in Fig.7. The average value of $`\beta _2`$ differs significantly from its most probable value (see Fig. 1) for temperatures above 1 MeV. The reason for such differences can be understood as follows. The free energy surface for $`{}_{}{}^{170}Er`$ at low temperatures ($`T<2`$ MeV) shows two well seperated minima. One of these minima lies at prolate deformation ($`F_p`$) and the other one is at oblate deformation ($`F_o`$). Their difference $`\mathrm{\Delta }F^0=F_oF_p`$ as a function of temperature is shown in Fig.8. For $`T>0.5`$ MeV, $`\mathrm{\Delta }F^0`$ decreases monotonically to zero. The initial rise of $`\mathrm{\Delta }F^0`$ can be attributed to the quenching of the pairing correlations; it is seen that without pairing, $`\mathrm{\Delta }F^0`$ decreases uniformly with temperature. The decrease in $`\mathrm{\Delta }F^0`$ enhances the probability of finding the nucleus with oblate deformation and hence the average deformation is smaller than the most probable one for temperatures below the transition temperature. The deviation of the average deformation from the most probable value is also governed by $`\mathrm{\Delta }\beta _2=\beta _2^p\beta _2^o`$. Here, $`\beta _2^p`$ and $`\beta _2^o`$ are the deformations corresponding to the prolate and the oblate minima. With rise in temperature, they come closer and ultimately merge as shown in Fig. 9. The persistence of $`\overline{\beta }_2`$ beyond the transition temperature is attributed to the asymmetry in the free energy surface as seen in the bottom panel of Fig.6. A similar behavior for $`{}_{}{}^{166}Er`$ is also seen.

In the mean field approximation, the thermal evolution of the most probable value of the deformation is seen to be associated with two phase transitions, the first being a superfluid phase to normal phase transition at $`T`$ 0.4 MeV and the other a weak second order phase transition in the shape at $`T`$ 2.7 MeV. It would be interesting to see how thermal fluctuations affect these phase transitions. In Fig. 10, the specific heat for the system $`{}_{}{}^{170}Er`$ is displayed as a function of temperature with inclusion of fluctuations. At low temperatures for $`T<1.0`$ MeV, there is no perceptible change in the specific heat and thus the transition from superfluid to normal phase is not affected by fluctuations. At a little higher temperature, however, the broad bump seen in Fig.3 at $``$ 2.7 MeV is seen to be shifted to $``$ 1.6 MeV and becomes also a little sharper. The phase transition in nuclear shape thus occurs at a lower temperature and is not as weak as seen in the mean field approximation. Normally, one would expect the phase transition to get diluted due to presence of fluctuations. Here one sees the opposite. The evolution of the free energy and the energy profiles with temperature showing two minima (one prolate and the other oblate) in deformation space is responsible for this observation. Here the role of thermal fluctuations becomes significant only in the intermediate temperature domain (temperature in the range 1 - 2.5 MeV, approximately) due to decrease in $`\mathrm{\Delta }F^0`$ and $`\mathrm{\Delta }\beta _2`$ (see Figs. 8 and 9). The delicate interplay of these quantities in the averaging augurs a change in the heat capacity signalling a somewhat sharper shape transition at a lower temperature compared to the most probable one.

The effect of thermal fluctuations of the shape coordinate $`\beta _2`$ on the inverse level density parameter $`K`$ for the system $`{}_{}{}^{170}Er`$ is shown in Fig.11. The full line and the dashed line correspond to results of our calculations for $`K_S`$ and $`K_E`$, respectively. The corresponding lines with circles represent those with inclusion of the fluctuations. It is seen that the fluctuations shift the minima for $`K_S`$ and $`K_E`$ from 2.7 MeV to 1.6 MeV, the corresponding transition temperature. Beyond the transition temperature, the upward slope for $`K`$ is also found to be higher with fluctuations included. Similar is the behavior for $`{}_{}{}^{166}Er`$.

## IV Conclusions

The relativistic mean-field approach together with pairing effects included through the constant gap approximation has been applied to study some properties of axially deformed $`{}_{}{}^{166}Er`$ and $`{}_{}{}^{170}Er`$ nuclei at finite temperature. The temperature dependence of the quadrupole shape for both the nuclei are found to be practically the same; similar is the case for the pairing gaps. The deformation remains close to its ground state value ($`\beta _2`$ 0.3) for $`T`$ less than 1.5 MeV and falls to zero quite sharply at $`T`$2.7 MeV. Pairing gaps vanish at $`T`$ =0.4 - 0.5 MeV leading to transition from a superfluid phase to normal phase. Compared with the results obtained in the nonrelativistic mean field framework, the relativistic calculations yield vanishing pairing gaps at practically the same temperature and the fall is a bit slower. The evolution of the heat capacity with temperature and the transition from superfluid to normal phase do not behave very differently. The nature of the shape transition obtained in the two approaches are also not too different except that in the present case, the approach of the deformation to zero is relatively slower. This slower fall results in a broad bump in the specific heat around the shape transition temperature ($`T`$ 2.7 MeV).

Fluctuations are expected to influence the phase transitions. To explore these aspects, we have taken into consideration thermal fluctuations in the $`\beta _2`$ degrees of freedom. The effects of fluctuations are found to be imperceptible below $`T`$ 1 MeV and therefore the pairing transitions are not affected. However, the influence on the shape transitions is found to be very significant; the transition temperature drops down from $``$ 2.7 MeV to $``$ 1.6 MeV with inclusion of thermal fluctuations. It also becomes a little sharper. The specific nature of the free energy surface in the deformation space and its evolution with temperature are responsible for such a behavior. Fluctuations in the $`\gamma `$ degrees of freedom would also have a role to play in the phase transitions mentioned. We have not included these presently as they are computationally too intensive.

The authors gratefully acknowledge Prof. P. Ring for providing them with the computer code to generate RMF solutions for axially deformed nuclei at zero-temperature.

Figure Captions:

* Variation of the most probable quadrupole deformation $`\beta _2`$ as a function of temperature for $`{}_{}{}^{166}Er`$ and $`{}_{}{}^{170}Er`$.
* Temperature evolution of neutron and proton pairing gaps for $`{}_{}{}^{166}Er`$ and $`{}_{}{}^{170}Er`$.
* Variation of specific heat as a function of temperature for $`{}_{}{}^{166}Er`$ and $`{}_{}{}^{170}Er`$.
* Temperature dependence of inverse level density parameters $`K_E`$ and $`K_S`$ determined from the excitation energy and the entropy, respectively.
* Spectra of the proton single-particle energies around the Fermi surface for $`{}_{}{}^{170}Er`$ at $`T=3.0`$ MeV in the RMF and the (P+Q) models. The single-particle energies $`ϵ`$ are shown relative to the respective chemical potentials.
* Plots for the free energy surface at various temperatures for $`{}_{}{}^{170}Er`$. The free energies at different temperatures are measured from the respective lowest minimum.
* Average value of quadrupole deformation as a function of temperature for $`{}_{}{}^{170}Er`$.
* The free energy difference between the prolate and the oblate minima as a function of temperature for $`{}_{}{}^{170}Er`$.
* The difference $`\mathrm{\Delta }\beta _2`$ between the deformations at the prolate and the oblate minima as a function of temperature for $`{}_{}{}^{170}Er`$.
* Variation of average specific heat as a function of temperature for $`{}_{}{}^{170}Er`$.
* The inverse level density parameter for the system $`{}_{}{}^{170}Er`$. The full and the dashed lines correspond to $`K_S`$ and $`K_E`$, respectively. The lines with filled circles represent the corresponding results with inclusion of thermal fluctuations.


format: text

text:
# Effects of a nano-sized filler on the structure and dynamics of a simulated polymer melt and the relationship to ultra-thin films

## Abstract

We perform molecular dynamics simulations of an idealized polymer melt surrounding a nanoscopic filler particle to probe the effects of a filler on the local melt structure and dynamics. We show that the glass transition temperature $`T_g`$ of the melt can be shifted to either higher or lower temperatures by appropriately tuning the interactions between polymer and filler. A gradual change of the polymer dynamics approaching the filler surface causes the change in the glass transition. We also find that while the bulk structure of the polymers changes little, the polymers close to the surface tend to be elongated and flattened, independent of the type of interaction we study. Consequently, the dynamics appear strongly influenced by the interactions, while the melt structure is only altered by the geometric constraints imposed by the presence of the filler. Our findings show a strong similarity to those obtained for ultra-thin polymer films (thickness $`100`$ nm) suggesting that both ultra-thin films and filled-polymer systems might be understood in the same context.

Significant enhancements in mechanical, rheological, dielectric, optical, and other properties of polymer materials can be obtained by adding fillers such as carbon black, talc, silica, and other inexpensive materials . Applications of filled polymers are diverse, ranging from automobile tires and bumpers to the rapidly expanding area of micro- and nano-electronic devices, which are continually challenged by the progressive miniaturization of devices . In particular, the growing ability to design customized nano-fillers of arbitrary shape and functionality provides an enormous variety of property modifications by introducing specific heterogeneity at the nanoscale . However, detailed knowledge of the effects of fillers on a polymer melt at the molecular level is lacking due to the difficulty of directly probing the polymer structure and dynamics in the vicinity of the polymer-filler interface. In this regard, molecular simulations provide an ideal opportunity for direct insight into filled materials. Additionally, the understanding of ultra-thin polymer films, which also have many important technological applications (e.g. paints, lubricants, adhesives, and electronic packaging), is a topic of continuing discussion ; the present results provide a framework in which to interpret experiments on filled polymers, and also possibly polymer thin films, which report both increases and decreases of the glass transition temperature $`T_g`$ , depending on the details of the system studied .

Our findings are based on extensive molecular dynamics simulations of a single nanoscopic filler particle surrounded by a dense polymer melt (Fig. 1). We simulate 400 chains of 20 monomers each (below the entanglement length), a total of 8000 monomers. The polymers are modeled as chains of monomers, which interact via a Lennard Jones (LJ) potential $`V_{\text{LJ}}(r)=4ϵ((\sigma /r)^{12}(\sigma /r)^6)`$. Additionally, bonded monomers are connected via a FENE anharmonic spring potential $`V_{\text{FENE}}=k(R_0^2/2)\mathrm{ln}(1(r/R_0)^2)`$ . Several authors have studied the pure (unfilled) melt , and the system has been shown to be a good glass-former. This type of “coarse-grained” model is frequently used to study general trends of polymer systems, but does not provide information for a specific polymer. The filler particle is modeled as an icosahedron, not unlike primary particles of graphitized carbon black . We assign ideal force sites at the vertices, at 4 equidistant sites along each edge, and at 6 symmetric sites on the interior of each face of the icosahedron, as shown in the right panel of Fig. 1. We tether a particle to each of these sites by a FENE spring, which maintains a relatively rigid structure but allows for thermalization of the filler. . We consider a filler particle with an excluded volume interaction only, as well as one with excluded volume plus attractive interactions (expected for many fillers), to determine which properties are a result of the steric constraints imposed by the filler, and which properties are affected by polymer-filler attraction. We choose the same parameters for the interaction potential for all filler force sites. Periodic boundary conditions are used in all directions.

Our model filler has several general features typical of a primary carbon black particle (a traditional filler) , as well as some newer nano-fillers : (i) it has a size of order 10 nm and (ii) it is highly faceted, but nearly spherical. The size of the facets is roughly equal to the end-to-end distance $`R_e`$ of the low molecular weight polymers comprising the surrounding melt (for Gaussian chains, $`R_e^2=6R_g^2`$, the radius of gyration). Since we aim to study the changes induced by a single filler on the nanoscopic level, and the implications of those changes for bulk properties, we also consider a pure dense melt for comparison.

We simulate the pure system at density $`\rho =1.0`$ at temperatures ranging from $`T=0.37`$$`1.0`$ . We simulate the filled systems in the range $`T=0.35`$$`1.2`$. Equilibration times range from $`5\times 10^2\tau ^{}`$ at the highest $`T`$ to $`2\times 10^4\tau ^{}`$ (approximately 40 ns in Argon units) at the lowest $`T`$; we use the rRESPA multiple time step algorithm to improve simulation speed . In order to compare the simulations of the filled system with the pure melt, we choose the box size so that the local density far from the filler deviates at most by 0.2% from the density of the pure melt; such a density difference would cause a change in $`T_g`$ in this model less than that reported in Fig. 2 . For attractive monomer-filler interactions, a box size $`L=20.4`$ (all lengths in units of $`\sigma _{\text{mm}}`$) satisfies this constraint at all $`T`$. In the non-attractive case, the characteristic first neighbor distance between the filler sites and monomers is $`T`$ dependent due to the lack of a unique minimum in the polymer-filler interactions. As a result, at each $`T`$ a different $`L`$ is required to achieve the correct $`\rho `$ at large distance from the filler. The box sizes range from $`L=20.49`$ at $`T=1.0`$, to $`L=20.6`$ at $`T=0.4`$.

To quantify the effect of the filler on $`T_g`$ — one of the most important processing parameters — and on dynamic properties, we first calculate the relaxation time $`\tau `$ of the intermediate scattering function

$$F(q,t)\frac{1}{NS(q)}\underset{j,k=1}{\overset{N}{}}e^{i𝐪[𝐫_k(t)𝐫_j(0)]},$$
(1)
where the sum is performed over all $`N`$ monomers and we normalize by the structure factor $`S(q)`$ such that $`F(q,0)=1`$. We define the value of $`\tau `$ by $`F(q,\tau )0.2`$. Relative to the pure system, we find that $`\tau `$ is larger at each $`T`$ for the attractive system (Fig. 2). Conversely, $`\tau `$ appears slightly smaller at low $`T`$ in the non-attractive system, but is nearly indistinguishable from the results for the pure system. The lines are fits to the Vogel-Fulcher-Tammann (VFT) form

$$\tau e^{A/(TT_0)}$$
(2)
where $`T_0`$ is typically quite close to the experimental $`T_g`$ value ; hence changes in $`T_g`$ are reflected in $`T_0`$. Consistent with the changes in $`\tau `$ relative to the pure melt, we find that $`T_0`$ increases in the system with attractive interactions, but clearly decreases in the system with only an excluded volume interaction. Thus, the effect of the steric hindrance introduced by the filler particle decreases $`\tau (T)`$ and $`T_g`$, in spite of the fact that monomers have a reduced number of directions in which to move, and hence degrees of freedom that aid in the loss of correlations. The fact that $`T_g`$ shifts in opposite directions for attractive versus purely excluded volume interactions demonstrates the importance of the surface interactions.

To elucidate how the local dynamics of the monomers are influenced by the filler, we examine the relaxation of the self (incoherent) part $`F_{\text{self}}(q,t)`$ of $`F(q,t)`$ as a function of the monomer distance from the filler. Monomers typically form layers near a surface ; we find well-defined monomer layers surrounding the filler, as seen in the density profile of Fig. 3. Hence, we split $`F_{\text{self}}(q,t)`$ into the contribution from each separate layer. Specifically, we calculate $`F_{\text{self}}^{\text{layer}}(q,t)`$ using the monomers located in each layer at $`t=0`$ such that $`F_{\text{self}}(q,t)=1/N_{\text{layers}}N_{\text{layer}}F_{\text{self}}^{\text{layer}}(q,t)`$, where $`N_{\text{layer}}`$ is the number of monomers in a given layer. We show $`F_{\text{self}}^{\text{layer}}(q_0,t)`$, as well as $`F_{\text{self}}(q_0,t)`$ for one temperature in Fig. 3. In the attractive system, the relaxation of the layers closest to the filler are slowest, consistent with the system dynamics being slowed by the attraction to the filler. Conversely, for the non-attractive system, we find that the relaxation of inner layer monomers is significantly enhanced compared to the bulk, consistent with the observed enhancement of the system dynamics. Preliminary results support the possibility that faster dynamics may also occur with attractive interactions, provided that the polymer-filler attraction is weaker than that of polymer-polymer interactions. The altered dynamics persist for a distance slightly less than $`2R_g`$ from the surface. Our results demonstrate that interactions play a key role in controlling $`T_g`$ and the local dynamics of filled polymers. We expect the role of interactions to be largely the same when many filler particles are present in the melt, but there will be additional effects on dynamic properties due to the more complex geometrical constraints.

We next turn our attention to any structural effect the filler has on the melt. Quantities such as the pair distribution function, average $`R_g`$, average $`R_e`$, and distribution of bond lengths and angles show no significant deviations from the pure system. However, by focusing on the dependence of $`R_g`$ (or $`R_e`$) on the distance $`d`$ from the filler surface, we find a change in the overall polymer structure near the surface. In Fig. 4, we show $`R_g^2`$, as well as the radial component from the filler center $`R_g^2`$ (approximately the component perpendicular to the filler surface) for both attractive and non-attractive polymer-filler interactions at one temperature. $`R_g^2`$ increases by about 50% on approaching the filler surface; at the same time $`R_g^2`$ decreases by slightly more than a factor of 2 for both attractive and non-attractive systems. The combination of these results indicates that the polymers become slightly elongated near the surface, and flatten significantly. Note that not all monomers belonging to a given “surface polymer” are located in the first surface layer, as depicted in Fig. 1. We also point out that the chains retain a Gaussian conformation near the filler surface . We find that the range of the flattening effect roughly spans a distance $`R_g`$ from the surface, and the results depend only weakly on $`T`$. We performed an additional simulation with double the attraction strength between the filler and polymers and did not find any significant further effect on the chain structure. The independence of the chain structure on the interaction suggests that the altered shape of the polymers is primarily due to geometric constraints of packing the chains close ($`dR_g`$) to the surface. For significantly stronger interactions, alteration of the chain structure is expected on theoretical grounds .

We next consider the implications that our results may have for studies of ultra-thin polymer films, where there is long standing debate on the role of interactions versus confinement on $`T_g`$ shifts , local melt dynamics , and melt structure . Our simulations allow us to address the effects of interactions with a surface, without the additional complication of confinement effects present in thin films. It is largely agreed that ultra-thin films with strongly attractive substrates increase $`T_g`$, while weak substrate interactions (or no substrate, as in freely standing films) lead to a downward shift of $`T_g`$; this is consistent with our results. This consistency is reasonable for fillers which have facets that are relatively smooth and large compared to $`R_g`$; for nanoscopic fillers, such as we study, it is surprising that a correspondence occurs even for $`R_g`$ close to the filler size. Not surprisingly, the magnitude of the shifts we observe depends on the relative quantities of polymer and filler; a greater filler concentration would have a more dramatic effect (as observed experimentally in refs. ). Insofar as the magnitude of effects depends only on the ratio of the surface to bulk monomers, the thickness of the film is analogous to the inverse of the concentration of the filler. This is consistent with the experimental observation that $`T_g`$ shifts are more pronounced as film thickness decreases. Recently there have been several experiments on segmental motion in both freely-standing and supported ultra-thin films . The observed segmental dynamics is consistent with a decreased $`T_g`$ find in calorimetric measurements . At this time, it is not clear whether a model with layers of different mobility is applicable to understand $`T_g`$ shifts of thin films ; however, the parallel behavior we observe between the thin films and our simulations of a filled melt support this viewpoint. Finally, the elongation and flattening of polymers we observe near the filler has been observed in thin-film simulations  as well as recent experiments ; the range of the effect found in Ref.  is quantitatively consistent with our results, which show the effect only for a range of roughly $`R_g`$, while the results of Ref.  observed flattening for film thicknesses $`6R_g`$. We also found, as in Ref. , that the chains retain a Gaussian structure near the surface. Thus our findings demonstrate that confinement is not a necessary ingredient for the observed changes in the dynamics and structure of polymers near surfaces. While our results provide strong support for interpreting the results for filled melts and ultra-thin films in the same framework, it is obvious that much care must be used when analyzing specific systems.

We wish to thank E. Amis, J. Douglas, Y. Gebremichael, C. Han, A. Karim, N. Lačević, and A. Nakatani, and W. Wu. We also thank the High Performance System Usage group at NIST for use of the SGI Origin Supercomputer. F.W.S. thanks the National Research Council for support.


text:
# Prompt muon contribution to the flux underwater

## I Introduction

Considerable literature exists on estimating the contribution to cosmic ray muon fluxes that arises from the decay of charmed hadrons . Current data on the high-energy atmospheric muon flux obtained with many surface and underground detectors are too conflicting to provide the means of probing charm production models (see, for example, ).

Both direct and indirect measurements of the atmospheric muon flux at sea level are limited to $`70`$ TeV for the vertical and to $`50`$ TeV for the horizontal. Statistical reliability of these data is still insufficient to evaluate the prompt muon contribution to the high-energy muon flux. Available energies and the accuracy of underground measurements are constrained because of the restricted size of detectors and the uncertainties in the local rock density. Deep-sea installations have substantial advantages just due to large detector volume and homogeneous matter. So it is relevant to discuss the potential of the large underwater neutrino detectors (AMANDA, Baikal), in the context of the prompt muon study, in future high-energy muon experiments.

In this paper, we present calculations on zenith angle dependence of the high energy underwater muon flux taking into consideration the prompt muon fraction obtained in one of the recent perturbative QCD models of Pasquali et al.  (hereafter pQCD) in which small-x behavior of the gluon distributions is probed. This pQCD calculation based on MRSD-  and CTEQ3  parton distribution functions (PDFs) includes the next-to-leading order (NLO) corrections to the charm production cross sections.

Perturbative QCD models differ in the PDF sets being employed in the NLO calculations and in the choice of renormalization and factorization scales. A dependence on these quantities of the vertical sea-level prompt lepton fluxes was studied in Refs. . The predictions of the pQCD model  are comparable to those of the earlier quark-gluon string model  and the recombination quark-parton one (see ). The muon spectra underwater obtained with the pQCD models and other types of the charm production models, the quark-gluon string model (QGSM) and the recombination quark-parton model (RQPM), were partly discussed in Ref. . Here we would like to focus on variations in the expected underwater muon fluxes caused by distinctions between the PDFs used. In addition, we compare the expected underwater muon flux to the zenith angle distributions measured with the Baikal  and the AMANDA  neutrino telescopes.

## II Sea-level muon fluxes

The atmospheric muon energy spectra and zenith angle distributions of the conventional ($`\pi ,K`$) muons, and the RQPM and the QGSM contribution, have been computed using the same nuclear cascade model . Let us glance over its key assumptions.

(i) The all-particle primary spectra and chemical composition are taken according to Ref. . Nuclei of the primary cosmic rays are treated as the composition of free nucleons.

(ii) Feynman scaling is assumed to be valid for hadron produced in collisions of hadrons with nuclei of the atmosphere.

(iii) The inelastic cross sections $`\sigma _{hA}^{\mathrm{𝑖𝑛𝑒𝑙}}`$ for interactions of a hadron $`h`$ ($`=p,n,\pi ^\pm `$) with a nucleus $`A`$ grow logarithmically with the energy:

$`\sigma _{hA}^{\mathrm{𝑖𝑛𝑒𝑙}}(E_h)=\sigma _{hA}^0+\sigma _A\mathrm{ln}(E_h/1\text{TeV}).`$

(iv) The hadron production in kaon-nucleus and in charmed hadron-nucleus collisions is neglected.

(v) Charged pion is considered stable in the kinetic stage of the nuclear cascade (not in the stage of muon production, to be sure).

(vi) Three-particle semileptonic kaon decays $`K_{\mu 3}`$ are taken into account.

The energy spectrum of the conventional muons at the sea level calculated for the vertical can be approximated (see ) for energies $`E_\mu 1`$ GeV by the formula:

$$𝒟_\mu ^{\pi ,K}(E_\mu ,0^{})=\{\begin{array}{cc}A_1E_\mu ^{f(y)},\hfill & E_1E_\mu E_2\hfill \\ A_2E_\mu ^{(1.791+0.304y)},\hfill & E_2E_\mu E_3\hfill \\ A_3E_\mu ^{3.672},\hfill & E_3E_\mu E_4\hfill \\ A_4E_\mu ^4,\hfill & E_\mu >E_4.\hfill \end{array}$$
(1)
Here $`E_\mu `$ in GeV, $`y=\mathrm{log}_{10}(E_\mu /1\text{GeV}),`$

$$f(y)=0.3061+1.2743y0.263y^2+0.0252y^3;$$

$`A_1=2.95\times 10^3,A_2=1.781\times 10^2,A_3=14.35,A_4=10^3`$ (cm-2s-1sr-1GeV-1); $`E_1=1`$ GeV, $`E_2=928`$ GeV, $`E_3=1587.8`$ GeV, $`E_4=4.1625\times 10^5`$ GeV.

The results of the calculations of the muon zenith-angle distributions at sea level are presented in Table I for high energies 1-100 TeV. The differential energy spectra (scaled by $`E_\mu ^3`$) of the conventional muons at sea level are shown (solid) in Fig. 1 for the vertical and near horizontal direction together with the data of the Nottingham spectrograph  (one point at $`E_\mu 1.3`$ TeV), the MUTRON spectrometer , and indirect measurements .

Open circles represent the MACRO best fit for the verical direction . (The detailed comparison between the calculated muon energy spectra for different zenith angles and the sea-level experimental data, in particular for large zenith angles, as well as calculations of other authors, is made in Ref. .)

Dashed and dotted lines in Fig. 1 correspond to the vertical muon flux including the prompt muon contributions calculated  with the CTEQ3 functions (pQCD-2) and the MRSD- set (pQCD-1). Line 1 (dashed) corresponds to the MRSD- set, line 2 (short dotted) corresponds to the CTEQ3 PDF, both with factorization and renormalization scales $`\mu _F=2\mu _R=2m_c`$, and with the charm quark mass $`m_c=`$ 1.3 GeV. As evident from the figure, the pQCD predictions depend strongly on the PDF.

For comparison there are also shown predictions of the charm production model of Volkova et al. (VFGS)  (thin) and the results obtained with the RQPM (dot-dashed) and the QGSM (dotted), both for the vertical direction (lower) and near the horizontal (up). These results enable one to make out the range of prompt muon flux predictions that overspread more than one order of magnitude at $`E_\mu 1`$ PeV. It is interesting to note that old QGSM predictions  coincide practically with those of the pQCD-2 up to $`600`$ TeV, while the RQPM flux appears to be close to the pQCD-1 one.

As is seen from Fig. 1, at $`E_\mu 20`$ TeV none of the above models but the VFGS is consistent with the data of MSU  and Frejus . Conversely, none of the charm production models under discussion contradict the LVD data . The VFGS, differing from the others in the extent of optimism, gives the greatest prompt muon flux that is scarcely compatible with the LVD upper limit .

The “crossing energy” $`E_\mu ^c(\theta )`$ (the energy around which the fluxes of conventional and prompt muons become equal) depends on the choice of the PDF set. The vertical crossing energy $`E_\mu ^c(0^{})`$ is about 200 TeV for the pQCD-1 model, which is close to the RQPM prediction ($`E_\mu ^c150`$ TeV). The vertical prompt muon flux predicted with the pQCD-2 model becomes dominant over the conventional one at the energies $`E_\mu 500`$ TeV. Therefore, in order for the differences between the pQCD models to be found experimentally one needs to measure muon energies above $`100200`$ TeV.

## III Muon fluxes underwater

Muon energy spectra and zenith angle distributions deep underwater are calculated using an analytical method  (see also Ref. ). By this method one can solve the problem of the muon passing through dense matter for an arbitrary ground-level muon spectrum and real energy dependence of differential cross sections of muon-matter interactions. The collision integral on the right-hand side of the muon transport equation describes the “discrete” energy loss of muons due to bremsstrahlung, direct $`e^+e^{}`$ pair production and photonuclear interactions.

In this paper the ionization energy loss and the part of the loss due to $`e^+e^{}`$ pair production with $`v<210^4`$ ($`v`$ is the fraction of the energy lost by the muon) is treated as a continuous one: that is, the corresponding item is subtracted from the collision integral and transferred to the left-hand side as a partial derivative with respect to energy of the mean energy loss rate multiplied by the muon flux.

The calculations of the prompt muon fluxes underwater at different zenith angles were performed with the parameterization of the sea-level muon differential spectra (pQCD-1,2) taken from .

Omitting details, we dwell on a factor that may be useful in correcting the underwater muon flux, provided that it is crudely estimated with the continuous energy loss approximation (see Ref. ). This factor is the ratio $`R_{d/c}`$ of the integral muon flux $`I_\mu ^{\mathrm{disc}}(E_\mu ,h,\theta )`$, computed for discrete stochastic) muon energy losses, to the flux $`I_\mu ^{\mathrm{cont}}(E_\mu ,h,\theta )`$ estimated with the continuous loss approximation.

In Table II the ratio $`R_{d/c}`$ is given as a function of the water depth and zenith angle for muon energies above $`10`$ GeV. As is seen the effect of discrete energy loss for the large depth is far from being small: $`R_{d/c}`$ is about 2 for the depth value of $`10`$ km w. e. The ratio is slightly affected by zenith-angular dependence of the sea-level muon flux. More precisely, the $`R_{d/c}`$ depends both on the “spectral index” of the muon flux and geometric factor of $`\mathrm{sec}\theta `$. The former varies weakly with zenith angle while the latter plays more important role in the $`R_{d/c}`$ defining the thickness of water layer $`X=h\mathrm{sec}\theta `$ that a muon overpasses ($`h`$ indicates the vertical depth in km).

For water the ratio $`R_{d/c}`$ as a function of the slant depth $`X`$ can be approximated with accuracy better than $`10\%`$ as

$`R_{d/c}`$ $`=`$ $`0.99+0.02X+6.7410^4X^3,X=112\mathrm{km};`$ (2)
$`R_{d/c}`$ $`=`$ $`1.43+0.054\mathrm{exp}[(X1.19)/3.64],X=1235\mathrm{km}.`$ (3)
The effect of the discrete loss increases as the muon energy grows. The energy dependence of the ratio $`R_{d/c}`$ is adequately illustrated by the following: for the depth of 12 km w. e. $`R_{d/c}`$2.5 at $`E_\mu =`$10 GeV and $`R_{d/c}`$4.0 at $`E_\mu =1`$ TeV.

In Fig. 2 we present a comparison between the expected muon vertical depth–intensity relation in water and the data obtained in underwater experiments (see for review ,), including recent measurements in the AMANDA-B4 experiment . The computation was performed with water parameters: $`\varrho =1g/cm^3`$, $`<Z>=7.47`$, $`<A>=14.87`$, $`<Z/A>=0.5525`$, $`<Z^2/A>=3.77`$. The muon energy loss per $`1g/cm^2`$ in ice is considered to be equal that in water but $`\varrho _{ice}=0.92g/cm^3`$. The calculations are presented for the muon residual energy (threshold of the detection) $`E_\mu `$1 GeV (solid) and $`E_\mu `$20 GeV (dashed). This difference needs to be considered especially for shallow depth.

Figs. 3, 4 show a comparison of the predicted muon zenith angle distribution (without considering the prompt muon contribution) with the measurements in the neutrino telescopes NT-36  and AMANDA .

The line in Fig. 3 presents the calculation for the muon threshold energy $`E_\mu =`$10 GeV at depth $`h=`$1.15 km. Our calculation is in reasonable agreement with the measurements of the NT-36 at all but the angle range $`8084^{}`$. In Fig. 4, the upper line relates to the flux at the depth $`h=`$1.60 km w. e. calculated for the muon residual energy $`E_\mu `$20 GeV, the lower one relates to $`h=`$1.68 km w. e. for the same energy threshold. The difference illustrates the possible effect of an uncertainty in determining the average “trigger depth”  (relating to the center of gravity of all hit optical modules in the AMANDA-B4 experiment). The computed angle distribution agrees fairly well with the AMANDA-B4 data including zenith angles $`\theta >70^{}`$.

The contributions of the ($`\pi ,K`$) and prompt muons underwater to zenith angle distribution at $`E_\mu >100`$ TeV calculated for four values of depths (of 1.15 to 4 km) are shown in Fig. 5. Here we present the results obtained with the pQCD-1 (dashed) and the pQCD-2 (dotted). It is interesting to note that the dashed line representing the pQCD-1 prompt muon contribution twice intersects the line of the conventional flux at $`h=1.15`$ km: near the vertical and at $`\theta 75^{}`$. This can occur because of the different zenith angle dependence of the conventional muon flux and the prompt muon one. And this means that at a depth of 1.15 km the nearly doubled muon event rate (for $`E_\mu >100`$ TeV) would be observed in the $`075^{}`$ range, instead of the rate expected due to conventional muons alone.

There is no intersection of the pQCD-2 line at $`h=1.15`$ km up to $`\theta 85^{}`$. The intersection point shifts to smaller zenith angles with increasing depth. For a depth of 2 km (nearly the AMANDA depth) it is possible to observe prompt muon fluxes that would be expected with the pQCD-2 at not too large angles ($`70^{}`$). It should be mentioned that the underwater prompt muon flux will be distorted in a large zenith angle region because the angle isotropy approximation considered for the predictions of the pQCD models is valid only at $`\theta 70^{}`$ and $`E_\mu 10^3`$ TeV.

The depth dependence of the muon flux underwater at zenith angle of $`78^{}`$ (Fig. 6) indicates that in the case of the pQCD-1 one can observe the doubling of the muon flux at the Baikal depth of 1.15 km for $`E_\mu 100`$ TeV. At a depth $`2`$ km the same takes place even with the lesser prompt muon flux predicted with the pQCD-2 model.

Fig. 7 shows muon integral energy spectra at a depth $`h=1.15`$ km (Baikal) and 2 km (AMANDA) and for $`\mathrm{cos}\theta =0.2`$ ($`\theta 78.5^{}`$). Also presented are the predictions of the prompt muon flux issued from the pQCD-1 (dashed) and pQCD-2 (dotted). The crossing energies $`E_\mu ^c(\theta )`$ at the AMANDA depth are less than the ones at the Baikal depth by a factor of $`3`$. In particular, the PQCD-1 model gives $`E_\mu ^c(\theta 78.5^{})30`$ TeV at $`h=2`$ km and $`E_\mu ^c(78.5^{})100`$ TeV at a depth of $`1.15`$ km. The same quantity calculated with pQCD-2 is 100 and 250 TeV respectively.

One can see (Fig. 7) that the AMANDA depth ($`2`$ km) gives, in a sense, the definite advantage in comparison with the Baikal one. Indeed, in the former case the assumed threshold energy is less, the muon flux difference between the pQCD-1 model and the pQCD-2 one is larger (up to two orders of magnitude), and the expected event rate remains approximately equal to the rate at the Baikal depth.

It should be pointed out that muon residual energies below $`10`$ TeV and zenith angle $`\theta 75^{}`$ would be available (see Ref.  for a discussion), in the above context, in future high-energy muon experiments with the NESTOR deep-sea detector  which is expected to deploy at a depth of about 4 km.

## IV Summary

Energy spectra and zenith angle distributions of the atmospheric muons at high energies have been calculated for the depths from 1 to 4 km that correspond to the depths of operation of large underwater neutrino telescopes. The estimation of the prompt muon contribution performed with the pQCD-1, 2 shows that the crossing energy $`E_\mu ^c`$ above which the prompt muon flux becomes dominant over the conventional one, is within the range of $`200`$ to $`500`$ TeV at sea-level, depending on the choice of the parton distribution functions. For the flux underwater at a zenith angle $`78^{}`$, the pQCD-1 model leads to the value $`E_\mu ^c30`$ TeV ($`h=2`$ km) and $`E_\mu ^c100`$ TeV ($`h=1.15`$ km). The corresponding crossing energies for the pQCD-2 model are $`E_\mu ^c100`$ and $`E_\mu ^c250`$ TeV.

The absolute value of the muon flux underwater around $`E_\mu ^c`$ depends on the charm production model. This circumstance enables, in principle, bounds to be put on the charm production cross section based on measurements of zenith angle distributions of the muon flux at high energies. In particular, PDF sets under discussion, the MRSD- and the CTEQ3, differ in the small-$`x`$ behavior of the gluon distribution: $`xg(x)x^\lambda `$, where $`\lambda 0.29÷0.35`$ for the CTEQ3 against $`\lambda =0.5`$ for the MRSD- set. (See Ref.  for the $`\lambda `$-dependence of the sea-level prompt muon flux.) These PDFs yield inclusive cross sections of charmed particles produced in nucleon-air collisions and charm production cross sections that diverge rapidly from each other with increasing energy. For muon energies above 100 TeV and for $`\mathrm{cos}\theta =0.2`$ these differences lead to the fact that prompt muon flux predicted with the pQCD-1 exceeds the flux arising from the pQCD-2 model by a factor of about 4 at $`h=1.15`$ km or about 5 at $`h=2`$ km.

In conclusion we outline three probable ways for solving of the prompt muon problem in the underwater experiments. First, one can measure zenith angle dependence of the muon flux in the energy region of $`50100`$ TeV (see Fig. 5): the expected event rate with the Baikal NT-200 is about $`200300`$ per year per steradian, supposing that the effective area of NT-200 is $`10^4`$ m2 for $`E_\mu 100`$ TeV .

Second, the flux with muon energies $`E_\mu 100`$ TeV measured as a function of depth (say, in depth region about $`0.81.2`$ km) at a given zenith angle ($`78^{}`$), could enable the charm production models to be discriminated (see Fig. 6) at the event rate level of about 200 per year per steradian.

And last, one can attempt to extract information on the prompt muon flux underwater from muon integral spectra being measured at a given depth and at a given zenith angle (Fig. 7). In this case the event rate is a factor 7 less than in the previous one (with the NT-200 capabilities). It should be pointed out that the AMANDA depth of $`2`$ km provides some advantage: the threshold energy is less, the muon flux difference between the pQCD-1 model prediction and the pQCD-2 one is greater, and the expected event rate remains approximately equal to that at the Baikal depth.

###### Acknowledgements.

We thank V. A. Naumov for helpful discussions and comments, S. Hundertmark and C. Spiering for kindly providing the table data on the muon zenith-angle distribution and the depth-intensity relation measured in the AMANDA-B4 experiment. This work is supported in part by the Ministry of Education of the Russian Federation under grant No. 015.02.01.04 (the Program “Universities of Russia – Basic Research”).


text:
# Evidence of two-electron tunneling interference in Nb/InAs junctions

## Abstract

The impact of junction transparency in driving phase-coherent charge transfer across diffusive semiconductor-superconductor junctions is demonstrated. We present conductivity data for a set of Nb-InAs junctions differing only in interface transparency. Our experimental findings are analyzed within the quasi-classical Green-function approach and unambiguously show the physical processes giving rise to the observed excess zero-bias conductivity.

PACS numbers: 74.80.Fp, 73.23.Ad, 73.40.Qv

The lack of available states below the superconducting energy gap prevents single-electrons from being injected across normal metal-superconductor (NS) junctions at low bias. A $`two`$-electron tunneling process known as Andreev reflection , however, allows charge transfer across interfaces with non-zero transparency. The two electrons involved in this process upon entering the superconductor form a Cooper-pair and must therefore be linked by specific energy and phase relationships .

The detailed nature of the N region in proximity to the superconductor can have a significant impact on junction conductivity . Indeed the interplay between phase-coherence in the N region and Andreev reflection can lead to enhanced or suppressed conductivity . Finite bias, temperature, and magnetic fields readily destroy phase coherence and consequently these phenomena are observed mostly around zero bias and at very low temperatures.

Experimentally, these coherent phenomena can be studied in the case of diffusive N electrodes. Beenakker  introduced the concept of reflectionless tunneling where the zero bias conductivity is predicted to present an unexpected linear dependence on interface transmission probability (analogous to a one-electron process). In the tunneling regime, Hekking and Nazarov  underlined the role of coherent phenomena in determining the constructive interference of the two electrons entering the superconductor. In the low transparency case, one can identify two regimes depending on the resistance of the disordered (phase coherent) region being larger or smaller than the barrier resistance. The excess conductivity around zero bias as a function of the latter increases or decreases, respectively. Unified theoretical descriptions of these transport phenomena are now available within the quasi-classical Green’s function  and the scattering matrix  approaches.

The first case (interface resistance limiting the conductivity) was invoked for a Cu-Al junction by Pothier et al. , the second (N-region resistivity dominant) for the superconductor-semiconductor junction results of Magnée et al.. Identification of the specific regime was obtained by analyzing the relevant junction parameters. In this article we take advantage of the qualitatively different dependence of conductivity on junction transparency to demonstrate two-electron tunneling interference in a set of InAs/Nb junctions. In our samples, in fact, zero-bias excess conductance develops for increasing the transparency of the junctions: as far as we know this is the first direct experimental demonstration of this characteristic trend. The data are quantitatively analyzed by solving the appropriate quasi-classical equations .

The samples consist of semiconductor-superconductor planar junctions, where the semiconductor part is a n-doped InAs (100) substrate with n = $`1.310^{18}`$ cm-3. The superconductor part is Nb. The relevant semiconductor parameters are the resistivity $`\rho =410^4\mathrm{\Omega }`$ cm and the mobility $`\mu =110^4`$ cm2 V-1 s-1. InAs is one of the best choices for hybrid Sm-S devices, mostly owing to the lack of a Schottky barrier at the contact with a metal. Oxides and other impurities, however, give rise to a residual interface $`barrier`$ region and several technological efforts have been made to improve the transparency of S-Sm junctions and maximize Andreev conversion . In our case, junctions were fabricated after annealing substrates at a moderate temperature (340 C) only marginally affecting the native oxide structure and its nonuniformities. This annealing procedure, however, minimizes the amount of water incorporated by the oxide and reduces the formation of Nb oxides. The superconductor film (100-nm-thick Nb) was e-beam evaporated following substrate annealing and without breaking the vacuum. Substrate temperature during Nb deposition was kept at $``$150 C to promote adhesion of the superconducting film.

Circular junctions (75 $`\mu `$m diameter) were fabricated by standard photolithographic techniques and wet chemical etching in a H2O:H2O2:NaOH solution. Back contacting was provided by metallizing the whole chip back. Samples were measured from 0.3 K to temperatures larger than the critical temperature $`T_c8.0`$ K of our Nb-on-InAs film and for static magnetic fields applied perpendicularly to the plane of the junctions using a 3He closed-cycle cryostat equipped with a superconducting magnet.

Figure 1 shows the differential conductance $`G(V)`$ of several contacts belonging to the same chip at $`T=0.32`$ K with no applied magnetic field. The wide range of conductance values observed is linked to the above mentioned oxide inhomogeneity. Notably high conductivity contacts exhibit zero-bias excess conductance. The latter effect was consistently more pronounced in more conductive junctions. As we shall argue in the following, more conductive junctions correspond to more transparent interfaces. The qualitative behavior of our junctions is therefore that studied by Hekking and Nazarov and excess conductance may be linked to the constructive interference of two-electron tunneling into the superconductor. This conclusion will be substantiated by the following quantitative analysis. High conductance curves are shown in a narrower range owing to heating effects at higher biases, but did exhibit a peak in proximity of $`V=\mathrm{\Delta }`$/$`e`$ (data not shown).

The experimental curves labeled (1) and (2) do not quite follow the expected behavior for a SN junction. The main discrepancies are in the detailed shape and in the position of the peaks at $`\mathrm{\Delta }`$. These deviations are related to the complex structure of the samples in the interface region. In order to describe electronic transport in our InAs/Nb junctions, we adopted the model first proposed in Ref. , namely we take into account the fact that the first layers of Nb are non-superconducting owing to the fabrication process (see Fig. 2). The potential drop still occurs at the metallurgical interface and therefore this interlayer region ($`N_{Nb}`$) can be considered at the same potential as the superconducting electrode. The thickness of $`N_{Nb}`$ ($`L`$) is expected to be of the order of a nm; nevertheless this interlayer strongly modifes the transport properties of the junctions. The various contributions to the current depend on the details of the sample close to the contact region and while the tunneling into the condensate (Andreev current) may not be very sensitive to the interlayer, the normal contribution strongly depends on the single-particle states available close to the interface .

The numerical evaluation of the current-voltage (I–V) characteristics is done in two steps. First we determine self-consistently the BCS gap in the Nb (close to the junction it will be suppressed due to proximity). Since we assume that there is no potential drop in this part of the structure we can use the equilibrium formalism. In the dirty limit, the quasiclassical retarded Green’s function $`\widehat{g}^R(x,E)`$ satisfies the Usadel equation

$$D_x(\widehat{g}^R_x\widehat{g}^R)+(iE+\mathrm{\Gamma }_{in}(x))[\widehat{\tau }_z+\widehat{\mathrm{\Delta }},\widehat{g}^R]=0,$$
(1)
with the constraints $`\widehat{g}^R\widehat{g}^R=1`$ and $`\text{Tr}\widehat{g}^R=0`$. The hat refers to the Nambu notation ($`\widehat{\tau }_z`$ and $`\widehat{\tau }_y`$ are the Pauli pseudospin matrices). The Fermi energy is at $`E=0`$. We assume that the inelastic scattering rate $`\mathrm{\Gamma }_{in}(x)`$ differs from zero only in the thin Nb interlayer close to the interface. The gap matrix $`\widehat{\mathrm{\Delta }}(x)=\mathrm{\Delta }(x)\widehat{\tau }_y`$ is determined self-consistently by means of

$$\mathrm{\Delta }(x)=\frac{1}{2}\lambda (x)_0^{\omega _D}𝑑E\mathrm{Im}[g_{12}(x,E)\mathrm{tanh}(E/2T)],$$
(2)
where $`\lambda (x)=\lambda \mathrm{\Theta }(x)`$ ($`\lambda =\text{arcosh}(\omega _D/\mathrm{\Delta }_{BCS})`$ is the BCS coupling constant, $`\mathrm{\Theta }(x)`$ is the Heaviside step function, $`\mathrm{\Delta }_{BCS}`$ is the BCS gap, $`\omega _D`$ is the Debye cutoff frequency) and $`T`$ is the temperature. The selfconsistent solution close to the interface enables us to determine the appropriate boundary condition for the determination of the I-V curves (along the lines of Volkov et al. )

$$I=\frac{1}{eR_N}𝑑ED(E)[f_0(E+eV/2)f_0(EeV/2)],$$
(3)
where $`R_N`$ the normal state resistance of the sample and $`f_0`$ the Fermi distribution function. The effective transmission coefficient of the structure

$$D(E)=\frac{1+r}{rM^1(E)+\frac{1}{d}_0^dM_T^1(x,E)𝑑x}$$

with $`M(E)=\mathrm{}g_{11}^R(0^+)\mathrm{}g_{11}^R(0^{})+\mathrm{}g_{12}^R(0^+)\mathrm{}g_{12}^R(0^{})`$ and $`M_T(x,E)=[\mathrm{}g_{11}^R(x,E)]^2+[\mathrm{}g_{12}^R(x,E)]^2`$ is determined by the solution of the Usadel equation. The interface between the diffusive region and $`Nb`$ is placed at $`x=0`$ ($`x=0^{+/}`$ indicated the normal/ superconducting side of the interface). The result of this procedure leads to the fits shown in Fig. 3. The theoretical curves in Fig. 3 are obtained assuming the ratio $`r=R_T/R_D=2.5`$ (where $`R_T`$ and $`R_D`$ are the interface and diffusive region contribution to the junction resistance, respectively). The fitting parameters in the $`N_{Nb}`$ region are the length $`L0.2\xi _{Nb}`$ (with $`\xi _{Nb}`$ the experimental coherence length in the dirty S) and the inelastic scattering rate $`\mathrm{\Gamma }_{in}=0.4\mathrm{\Delta }_{BCS}`$. All the curves of the set are fitted with the same parameters, and only the temperature is varied according to its experimental value.

Both $`L`$ and $`\mathrm{\Gamma }_{in}`$ are consistent with the experimental expectations. From the measured value of $`T_c`$ it was possible to extract  $`\xi _{Nb}8.0`$ nm for our Nb-on-InAs film, obtaining in this way a value of the interlayer length $`L1.6`$ nm consistent with the estimated thickness of the oxidized Nb interface region . The large value of $`\mathrm{\Gamma }_{in}`$ is related to the fact that superconductivity is locally suppressed (we stress again that the rest of the Nb is ideal). The presence of the thin layer of ”normal” Nb proved crucial to obtain a good quantitative agreement with the experimental data. In fact, at finite voltages the measured differential conductance is always larger than that predicted by a simple S-I-N system. This increased conductivity stems from the presence of the $`N_{Nb}`$ layer and the availability of states for the single-particle channel of the current. The latter are the subgap states of the proximized normal layer.

From the already given experimental parameters of InAs, we obtained a very good agreement between theory and experiments with $`d=200`$ nm (which corresponds roughly to the estimated coherence length in the normal region) and a mean free path $`l_e=220`$ nm at $`T=0.3`$ K . The fitting parameter $`r`$ allows us to quantitatively show that for the junction corresponding to curve (1) the weight of the tunnel-barrier resistance is predominant with respect to the resistance of the diffusive region. A quantitative fit for curve (2) was also obtained (data not shown). As anticipated above this occurred for a lower value of the parameter $`r`$ ($`r=2`$ compared to $`r=2.5`$ for junction (1), reflecting the larger transparency of the interface in the corresponding contact). Although the various junctions are nominally identical, the native oxide nonuniformities can easily lead to the observed variations. The parameter values extracted from the fits show that our junctions are in an intermediate case as compared to the limit discussed in Ref. . Though not deep in the asymptotic region, our analysis allows to extend the main qualitative features also to the range of parameters typical of the samples considered in this work.

The zero-bias conductance observed in the curves (3), (4), and (5) in Fig. 1 occurs on an increasing normal-state conductance value. The latter variation is influenced also by a concomitant variation of the junction effective area . For the present analysis, however, we must focus on the increased zero-bias contribution which is dominated by the junction-transparency variation. This corresponds to the Hekking-Nazarov regime, in which increasing $`\mathrm{\Gamma }`$ leads to a more pronounced conductance peak at $`V,H=0`$. In terms of our parameters this means that $`r`$ lowers. In fact, this trend is also present in our model calculations: by lowering the barrier resistance the zero-bias excess conductance develops. For this set of curves, however, we were unable to obtain a quantitative agreement with the experimental data.

Further proof on the nature of the enhanced zero-bias conductivity can be gained experimentally by analyzing its temperature and magnetic field dependence. As mentioned above two-electron interference effects can be easily broken by voltage, temperature and magnetic field. Figure 4 shows the temperature and the magnetic-field (inset) dependence of the zero-bias conductance in junction (4). Cut-off values for temperature ($`T1.0`$ K) and magnetic field ($`B6`$ mT) are similar to those reported by other experimental groups for such effects . Additionally we remark that at $`T=0.32`$ K the width of the zero bias peak is of the order of $`k_BT/e=27`$ $`\mu `$eV.

In conclusion, we have analyzed the influence of junction transparency on phase-coherent conductance in superconductor-semiconductor junctions. The dependence of the zero-bias excess conductance on interface transparency was directly observed and junction parameters estimated based on a quasi-classical Green’s function approach.

Very useful discussions with C. W. J. Beenakker and R. Raimondi are gratefully acknowledged. The present work was supported by INFM Sezione E under the PAIS project Eterostrutture Ibride Semiconduttore-Superconduttore.


text:
# 1 Introduction

## 1 Introduction

Why a survey for local QSOs?
Because in the epoch of 2dF and SDSS, of the thousands of QSOs at $`z2`$, it is important to have precise information about the spatial density and the clustering of nearby AGN to provide zero-point and leverage for the study of the QSO evolution. As in other fields, it is paradoxical that we know much better the properties of the high-z QSO population than the local one. The technical problem is however not trivial since, as the PG survey (Schmidt & Green 1983) made clear, it is necessary to cover the whole sky to obtain sufficient statistics where the present data are unsatisfactory: at $`z<0.3`$ and magnitudes around $`1315`$. Such a task has become relatively simpler in recent times, thanks to the availability of large high-quality databases.

## 2 The Photometric Database

What do we need ?
First of all, since our goal is the study of the optical luminosity function (OLF) of QSOs, we need a homogeneous database of optical fluxes covering the whole sky. Of various existing possibilities: APM111http://www.ast.cam.ac.uk/$``$apmcat/, ROE/NRL222http://xweb.nrl.navy.mil/www\_rsearch/RS\_form.html, USNO333http://archive.eso.org/skycat/usno.html, GSC444http://www-gsss.stsci.edu/gsc/gsc.html, DSS555http://www-gsss.stsci.edu/dss/dss.html, none turned out to be entirely satisfactory. After a careful analysis of their sensitivity, completeness, accuracy, we have chosen a combination of the catalogs GSC, USNO and DSS:

1. in the Northern hemisphere objects from the GSC catalog with $`11.0<V_{GSC}14.5`$. The relation between the $`V_{GSC}`$ band and the corresponding Johnson $`V`$ turned out to be: $`V_{GSC}=V0.21`$, with $`\sigma _V=0.27`$ $`mag`$. Seven ($`4\%`$) of the 183 Landolt stars used for this comparison were not found in the GSC catalog.
2. in the Northern hemisphere objects from the USNO catalog with $`13.5<R_{USNO}15.4`$. The relation between the $`R_{USNO}`$ band and the corresponding Johnson-Kron-Cousins $`R`$ turned out to be $`R_{USNO}=R_{JKC}+0.096`$, with $`\sigma _{R_{USNO}}=0.27`$ $`mag`$.
3. in the Southern hemisphere we have derived $`B_J`$ magnitudes from the Digitized Sky Survey (DSS). Small scans (“postage stamps” of $`2^{}\times 2^{}`$) of each object of interest and of 20-50 surrounding objects with known GSC $`B_J`$ magnitudes were extracted from the DSS. The magnitude of the object of interest was then calibrated against the GSC objects. In this way a $`\sigma _{B_J}`$ of $`0.10`$ $`mag`$ was obtained in the interval $`12.0<B_J<15.5`$.

## 3 The QSO Selection

The second step is to select among the millions of objects in the magnitude range of interest the few hundreds of QSOs. We have used the RASS Bright Source Catalogue (RASS-BSC, Voges et al. 1999) to compute the X-optical color, $`\alpha _{ox}`$ which provides a key signature of the AGN phenomenon. If we take from the Véron & Véron catalog (1998) all the non X-ray selected bright QSOs with $`0.04<z<0.3`$, assign to them the optical magnitudes of our databases and plot them vs. the corresponding RASS fluxes, we obtain the diagram of Fig. 1. We can see that, if we select objects with $`\alpha _{ox}<1.7`$, the incompleteness does not depend on the optical flux: we lose a fraction of about $`20\%`$ of the QSOs, which can be accounted for, and we obtain a reasonably short list of candidates: 520 in the North over 8000 sq.deg. and 301 in the South over 5600 sq.deg.

Conveniently, more than $`40\%`$ of the objects have already an identification in SIMBAD666http://simbad.u-strasbg.fr/Simbad or NED777The NASA/IPAC Extragalactic Database (NED) is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. For the remaining ones we are taking spectra with the telescopes at the Asiago, La Silla and Kitt Peak observatories. The success rate is around $`50\%`$. We have observed $`80\%`$ of the candidates in the North, $`50\%`$ in the South, collecting a total of 290 QSOs. At completion we expect to produce a sample of about 450 QSOs.

## 4 Results of the Asiago-ESO/RASS QSO Survey

Preliminary results have been reported in Grazian et al. (2000). The optical counts have been computed in the $`B_J`$ band, correcting our V and R magnitudes according to the average colors of the sample, and are shown in Fig. 2. They confirm the well known incompleteness of the PG survey (about a factor 3) and agree with the extrapolation of other recent surveys to somewhat fainter magnitudes.

The derived OLF confirms and strengthens the claim of La Franca & Cristiani (1997, LC97), i.e. a significant departure of the evolution at low-z from the pattern generally described with a Pure Luminosity Evolution (PLE).

It is debatable whether the best parameterization of the data is provided by a Luminosity Dependent Luminosity Evolution (e.g. LC97) or by a slow-down of the evolutionary rate and/or a mild density evolution at low-z (see Vittorini & Cavaliere, this conference). To obtain a more physical insight it is however advisable to put these results in the context of the models of formation and evolution of galactic structures, for evidence is mounting that the QSO activity, the growth of super-massive black holes and the formation of spheroids are closely linked phenomena.

On the one hand, the increased QSO space density at low-z corresponds better to the predictions of semi-analytical models (Kauffmann & Haehnelt 2000), in which the recurrent activity of short-lived QSOs is related to the merging of halos, the availability of gas for the accretion and its timescale. On the other hand, a weak density evolution and a smoother overall shape of the OLF at low-z are natural predictions of the models of Cavaliere & Vittorini (2000) arguing that at $`z<3`$ efficient black hole fueling is triggered by the encounters of a gas-rich host with its companions in groups.


text:
# Prediction of Astrophysical Reaction Rates: Methods, Data Needs, and Consequences for Nucleosynthesis Studies

## Introduction

Extensive reaction networks have to be employed in the investigation of nuclear energy generation and nucleosynthesis processes in astrophysics. Since stellar and explosive burning involves a considerable number of unstable isotopes which are currently unaccessible by experiments, the prediction of astrophysical reaction rates by means of nuclear reaction theory is unavoidable. There has been progress in the theoretical approaches, especially in the modelling and prediction of specific nuclear properties required for the determination of nuclear reaction rates. However, there is still need for experimental studies testing the predictions and providing more data to further improve theoretical models. After a discussion of statistical model inputs, this will be briefly addressed in the concluding section of the paper.

Considering the still remaining uncertainties in the prediction of nuclear reaction rates, it is of great interest to investigate the sensitivity of abundance yields to variations in the rates. This is also important if one wants to disentangle stellar physics and nuclear effects in the comparison of models which differ in both aspects. A comparison of two sets of reaction rates by employing them in the same stellar model of WW95  is presented in the second part of the paper.

## The Statistical Model

In general, the cross section will be the sum of the cross sections resulting from compound reactions via an average over overlapping resonances (HF) and via single resonances (BW), direct reactions (DI) and interference terms:

$$\sigma (E)=\sigma ^{\mathrm{HF}}(E)+\sigma ^{\mathrm{BW}}(E)+\sigma ^{\mathrm{DC}}(E)+\sigma ^{\mathrm{int}}(E).$$
(1)
Depending on the number of levels per energy interval in the system projectile+target, different reaction mechanisms will dominate RTK97 . Since different regimes of level densities are probed at the various projectile energies, the application of a specific description depends on the energy. In astrophysics, one is interested in energies in the range from a few tens of MeV down to keV or even thermal energies (depending on the charge of the projectile). It has been shown RTK97  that the majority of nuclear reactions in astrophysics can be described in the framework of the statistical model (HF) HF52 . This description assumes that the reaction proceeds via a compound nucleus which finally decays into the reaction products. With a sufficiently high level density, average cross sections

$$\sigma ^{\text{HF}}=\sigma _{\text{form}}b_{\text{dec}}=\sigma _{\text{form}}\frac{\mathrm{\Gamma }_{\text{final}}}{\mathrm{\Gamma }_{\text{tot}}}$$
(2)
can be calculated which can be factorized into a cross section $`\sigma _{\text{form}}`$ for the formation of the compound nucleus and a branching ratio $`b_{\text{dec}}`$, describing the probability of the decay into the channel of interest compared with the total decay probability into all possible exit channels. The partial widths $`\mathrm{\Gamma }`$ as well as $`\sigma _{\text{form}}`$ are related to (averaged) transmission coefficients, which comprise the central quantities in any HF calculation.

Many nuclear properties enter the computation of the transmission coefficients: mass differences (separation energies), optical potentials, GDR widths, level densities. The transmission coefficients can be modified due to pre-equilibrium effects which are included in width fluctuation corrections T74  (see also RTK97  and references therein) and by isospin effects. It is in the description of the nuclear properties where the various HF models differ. A choice of what is thought of being the currently best parametrizations is incorporated in the new HF code NON-SMOKER rt98 , which is based on the well-known code SMOKER thi87 .

## A Reaction Rate Library

Utilizing the NON-SMOKER code, cross sections and reaction rates for reactions with nucleons, $`\alpha `$ particles or $`\gamma `$ rays in entrance and exit channels, respectively, were calculated for all targets between proton and neutron drip line in the range $`9<Z<84`$. Tabulated cross sections and rates can be found at http://quasar.physik.unibas.ch/~tommy/reaclib.html. Analytic fits to these rates, along with further information, can be obtained as an electronic file from the authors or on-line from Atomic Data and Nuclear Data Tables. A selection from these fits is published in RT00 .

In all applications, these rates should be supplemented or replaced with experimental rates as they become available. Such a combination of theoretical and experimental rates is provided, e.g., in the REACLIB compilation. Currently, a new version is being compiled, in which the theoretical rates presented here will be included. Latest information on REACLIB can be found on the WWW at http://ie.lbl.gov/astro.html. Further details on the NON-SMOKER code are presented at http://quasar.physik.unibas.ch/~tommy/reaclib.html.

## Reaction Rate Sensitivity of Nucleosynthesis in Type II Supernovae

When comparing results from different supernova models one faces the difficulty caused by the fact that it is hard to differentiate between influences of differing reaction rate sets and different stellar physics. We tried to segregate the abundance differences between the two models of WW95  (WW95) and TNH96  (TNH) existing because of the dichotomy of stellar models from those reflecting purely the choice of nuclear physics. For that purpose, hybrid calculations were performed, using the same stellar evolution code as in WW95  but with rates from both models. In addition to helping to understand why calculations of the two groups differ, the use of independent rate sets in identical stellar models helps determine the nuclear physics portion of the error bar one should assign to nucleosynthesis studies of this sort. The cause of the differences in the theoretical rates was further investigated, pointing to possibilities for future improvements of rate predictions. In the following, the findings are briefly summarized. A very detailed account of the work can be found in HWWRT99 .

### The Rate Sets

The reaction rates utilized in WW95 were those of WFH78  (WFHZ), TNH used thi87  (TAT). As examples, Figs. 1 and 2 show a comparison of the two sets to each other and to experimental values for 30 keV neutron capture and proton capture at $`T_9=3`$.

Typical differences at astrophysically interesting temperatures are less than a factor of two. There are individual cases, however, where the difference exceeds a factor of 10. Some of the larger differences occur for reactions where scarce experimental information is available and different assumptions were made regarding the photon transmission function, for example, ($`\alpha `$,$`\gamma `$) reactions on $`Z=N`$ nuclei. Different assumptions were also made about the particle transmission functions, nuclear partition functions, and level densities. More modern and complete data used in the TAT rates makes them superior in cases where the partition function is important. WFHZ used an equivalent square well with empirical reflection factors; TAT used a more detailed optical model. Given the quite different values for, e.g., the neutron and proton transmission function, it is perhaps surprising that the rates differ so little. This is because the relevant temperatures for explosive burning are high. For incident particles in the Gamow window, the deviations in the particle transmission functions are typically smaller than a factor of two. In addition, higher partial waves contribute. A comparison of rates at a lower temperature would have revealed larger discrepancies.

Compared to experiment, both sets of theoretical rates give similar agreement, typically to a factor of two. The standard deviations between the two theoretical sets and cross section data are almost identical. In summary, the two rate sets have comparable merit when compared to experiment. All the authors of this paper agree that the new rate set, the “NON-SMOKER” set, will be preferable to both TAT and WFHZ and will be adopted by both groups (WW and TNH) for future work.

### Results and Conclusions

The comparison of the yields obtained with the two reaction rate sets in the WW95 model is shown in Figs. 3 and 4 for a 15 and a 25 $`M_{}`$ supernova, respectively.

When the two current rate sets are included in otherwise identical stellar models we find that the nucleosynthesis, with some interesting exceptions, is not greatly changed. For example, only about a dozen (out of 70) stable isotopes in the mass range 12 to 70 have nucleosynthesis that differ by over 20% in two supernovae of 15 $`M_{}`$ that use the same rate for 12C($`\alpha ,\gamma `$)16O. It can, however, be noticed that most of these isotopes - with one exception 44Ti - are products of hydrostatic burning where individual reaction flows are governed by the cross sections involved. Nevertheless, none differ by more than a factor of 1.7. Given the significantly larger differences that exist in individual reaction rates, one may wonder at the robust nature of the final nucleosynthesis. We see three major causes.

First, as the star burns and becomes hotter, the nuclear flow follows the valley of $`\beta `$ stability making heavier nuclei as it goes. In doing so, it follows the path of least resistance – those reactions having the largest cross section for a nucleon or $`\alpha `$ particle reacting with a given nucleus. These large cross sections are reasonably well replicated by any calculation, normalized to experiment, that treats the Coulomb barrier and photon transmission function approximately correctly. Large differences may exist in rate factors for reactions that are in competition, especially a small channel in the presence of a large one, but these small channels are frequently negligible, at least for the major abundances while they can cause larger differences when one is interested in the abundances of trace isotopes, e.g. in isotopic anomalies of meteoritic inclusions.

If one is, however, interested in accurate abundance predictions resulting from these smaller flows in hydrostatic burning stages, these can in most cases only be obtained by improving the reliability of the cross sections (and reaction rates) that determine these weak flows on light and intermediate mass nuclei. As new experimental information becomes available, a continuous improvement is therefore highly warranted.

Second, beyond oxygen burning (nuclei heavier than calcium), nucleosynthesis increasingly occurs in a state of full or partial nuclear statistical equilibrium. There the abundances are given by binding energies and partition functions. As long as the “freeze-out” is sufficiently rapid, individual rates are not so important.

Third, the reaction rates varied here were only those theoretical values from Hauser-Feshbach calculations for intermediate mass nuclei, i.e., nuclei heavier than magnesium. The really critical reaction rates are, for the most part, those below magnesium. These reactions, like e.g., 12C($`\alpha ,\gamma `$)16O, govern the energy generation, major nucleosynthesis, and neutron exposure in the star. The rest are perturbations on these dominant flows.

This is not to say, however, that the nuclear and stellar details of heavy element synthesis are now well understood. Differences in the stellar model may account not just for 20% variation, but orders of magnitude. That is, uncertainty in stellar physics – especially the treatment of convection and how it is coupled (or not coupled) to the nuclear network – accounts for most of the differences in current nucleosynthesis calculations – provided such calculations use the same nuclear reaction rates below magnesium.

Even in a perfect stellar model though, there will still be interesting nuclear physics issues. Stellar nucleosynthesis is becoming a mature field rich with diverse and highly detailed observational data. The “factor of two” accuracy that was adequate in the past may not do justice to the observations of the future. There are many individual cases where the nuclear physics uncertainty is still unacceptably large. We point out just a few examples for which experimental information would be of interest.

## Nuclear data needs

The suppression of radiative capture reactions into self-conjugate (isospin zero) nuclei is very uncertain. Past Hauser-Feshbach calculations have adopted empirical factors for this suppression. The new NON-SMOKER rates include a significantly improved treatment RGW99 . $`\alpha `$-capture reactions, like 24Mg($`\alpha ,\gamma `$)28Si, 28Si($`\alpha ,\gamma `$)32S, …, 44Ti($`\alpha ,\gamma `$)48Cr, are very important to nucleosynthesis in oxygen and silicon burning. The reaction 40Ca($`\alpha ,\gamma `$)44Ti also directly affects the synthesis of 44Ti. Modern accurate determinations of most of the reaction rates are missing (as well as (p,$`\gamma `$) reactions into the same nuclei). Measurements here would be most welcome.

The Hauser-Feshbach rates are also only as good as the local experimental rates to which the necessary parameters of the calculation are calibrated. In that regard we would point out the near absence of charged particle reaction rate data for $`A>70`$, even for stable nuclei. Charged particle reactions are important, especially on unstable nuclei, at significantly higher atomic weights in the $`r`$ process and in the $`p`$ process.

Extended systematics of other nuclear properties, such as level densities (especially around magic nucleon numbers), low energy behavior of the GDR, and optical potentials, would be highly appreciated for stable as well as unstable nuclei.

### Acknowledgements

This work has been supported by the U.S. NSF (AST 96-17494, AST 96-17161, AST 97-31569, PHY-74-07194) and the Swiss NSF (20-47252.96, 2000-53798.98). T. R. is a PROFIL fellow of the Swiss NSF.