text:
# Charged-current neutrino-²⁰⁸𝑃𝑏 reactions

\[

## Abstract

We present theoretical results on the non flux-averaged $`{}_{}{}^{208}Pb(\nu _e,e^{})^{208}Bi`$ and $`{}_{}{}^{208}Pb(\nu _\mu ,\mu ^{})^{208}Bi`$ reaction cross sections, obtained within the charge-exchange Random-Phase-Approximation. A detailed knowledge of these cross sections is important in different contexts. In particular, it is necessary to assess the possibility of using lead as a detector in future experiments on supernova neutrinos, such as OMNIS and LAND, and eventually detect neutrino oscillation signals by exploiting the spectroscopic properties of $`{}_{}{}^{208}Bi`$. We discuss the present status on the theoretical predictions of the reaction cross sections.

\]

The study of reactions induced by neutrinos on nuclei is at present an active field of research. A detailed knowledge of the reaction cross sections is interesting for different domains, going from high energy physics to astrophysics . For example, they are necessary in the interpretation of current experiments on neutrinos as well as in the evaluation of possible new detectors for future experiments. The importance of neutrino-nuclei reactions in astrophysical processes, such as the r-process nucleosynthesis, is also being attentively studied . In particular, $`\nu Pb`$ reactions have attracted much interest recently. Lead has been used as a shielding material in the recent experiments on neutrino oscillations performed by the LSND collaboration  so that estimates of the $`\nu Pb`$ reaction cross sections are necessary for the evaluation of backgrounds in these experiments; also projects on lead-based detectors , such as OMNIS  and LAND , are being studied for the purpose of detecting supernova neutrinos. These detectors might provide information on neutrino properties, such as oscillations in matter  or the mass  by measuring the time delay and/or spreading in the neutrino signal  as well as help in testing supernova models. From the practical point of view, lead-based detectors seem to present several of the characteristics required to be supernova observatories, namely high sensitivity to neutrinos of all flavors, simplicity, reliability with inexpensive materials . Large cross sections for neutrinos in the supernova energy range are also an important condition since they determine the possible rates and therefore the maximum observable distance. Actually, $`\nu `$-nucleus reaction cross sections increase strongly with the charge of the nucleus. For example, if the neutrinos come from the Decay-At-Rest (DAR) of $`\mu ^+`$, the cross sections of the flux-averaged charged-current (CC) reaction $`\nu _e+_ZX_N_{Z+1}X_{N1}^{}+e^{}`$ goes from about $`14.10^{42}`$ cm2 for $`{}_{}{}^{12}C`$ , to $`2.5610^{40}`$ cm2 in $`{}_{}{}^{56}Fe`$  and is estimated to be $`3.6210^{39}`$ in $`{}_{}{}^{208}Pb`$ . Besides these practical features which are essential in the choice of the nucleus to use to detect neutrinos, another important feature is the spectroscopic properties which may suggest attractive signals of supernova neutrino oscillations. In , for example, it has been shown that the measurement of events where two neutrons are emitted by $`{}_{}{}^{208}Bi`$ excited in the reaction $`\nu _e+^{208}Pb^{208}Bi+e^{}`$ is both flavor-specific and very sensitive to the mean energy of the $`\nu _e`$. In case when $`\nu _\mu ,\nu _\tau \nu _e`$ oscillations take place, the hotter $`\nu _e`$ would increase the number of two neutron events by a factor of forty . Another possible signal has been proposed in , that is that the energy distribution of the neutrons emitted in the same CC reaction should have a peak at low energy more or less pronounced according to whether the oscillations occur or not. This peak would come from the excitation of a peak at around 8 MeV in the Gamow-Teller strength distribution. (One should however note that this peak has never been observed experimentally). Both the estimate of the CC $`\nu Pb`$ reaction cross section in  and the microscopic calculations of  show that a possible oscillation signal relies strongly on the knowledge of the spectral properties of $`{}_{}{}^{208}Bi`$. In fact, the CC reaction cross section induced by $`\nu _e`$ scales almost as the square of the electron energy and is particularly sensitive to the detailed structure of the excitation spectrum as was already pointed out for the case of $`{}_{}{}^{12}C`$ . It is then important either to get the cross sections directly from the experiment or/and to obtain different theoretical estimates in order to know the theoretical uncertainties and how they affect the reaction cross sections. This is crucial when the impinging neutrino energy increases because not only the allowed Gamow-Teller (GT) and Isobaric Analogue State (IAS) contribute significantly to these cross sections but also forbidden transitions, of first-, second-, third-order (which are not very well known experimentally).

In this paper, we present new theoretical results for the CC $`\nu _e+^{208}Pb^{208}Bi+e^{}`$ reaction cross section. Our calculations, as opposed to , are performed in a self-consistent charge-exchange Random-Phase-Approximation (RPA) with effective Skyrme forces. Contrary to all the previously published calculations, we present non flux-averaged cross sections, obtained for both low-energy $`\nu _e`$ and high-energy $`\nu _\mu `$. These reaction cross sections given as a function of neutrino energy span a large energy range. They can be used to convolute with different neutrino fluxes in various contexts, for example for future experiments with astronomical neutrinos which are at present under study, for the very recent terrestrial experiments such as the LSND ones  to estimate the background, or in the r-process nucleosynthesis.

We will emphasize the importance of the contribution of forbidden transitions and how it evolves as a function of neutrino energy. This is often not taken into account in many present r-process nucleosynthesis calculations and so the neutrino-nuclei cross sections are underestimated (in  only the importance of first forbidden transitions in neutron-rich nuclei was emphasized).

We will compare our results with presently available calculations . With this aim, we will present two different flux-averaged cross sections, where the neutrino fluxes are given by either the DAR of $`\mu ^+`$ and Decay-In-Flight (DIF) of $`\pi ^+`$; or by a Fermi-Dirac spectrum for a supernovae explosion. Finally, we will discuss our results in relation to the suggested possible oscillation signals that would use the spectroscopic properties of $`{}_{}{}^{208}Bi`$.

The general expression for the differential cross section as a function of the incident neutrino energy $`E_\nu `$ for the reaction $`\nu _l+^{208}Pbl+^{208}Bi`$ ($`l=e,\mu `$) is

$$\sigma (E_\nu )=\frac{G^2}{2\pi }cos^2\theta _C\underset{f}{}p_lE_l_1^1d(cos\theta )M_\beta ,$$
(1)
where $`Gcos\theta _C`$ is the weak coupling constant, $`\theta `$ is the angle between the directions of the incident neutrino and the outgoing lepton, $`E_l=E_\nu E_{fi}`$ ($`p_l`$) is the outgoing lepton energy (momentum), $`E_{fi}`$ being the energy transferred to the nucleus, $`M_\beta `$ are the nuclear Gamow-Teller and Fermi type transition probabilities .

In a nucleus as heavy as $`Pb`$ the distortion of the outgoing lepton wavefunction due to the Coulomb field of the daughter nucleus becomes large and affects the integrated cross section considerably. In our treatment of this effect we follow the findings of ref.. In ref. it is found that the “Effective Momentum Approximation” (EMA) works well for high energy neutrinos. This approximation consists in using an effective momentum $`p_l^{eff}=\sqrt{E_{eff}^2m^2}`$ where $`E_{eff}=EV_C(0)`$ ($`V_C(0)`$ is the Coulomb potential at the origin) in calculating the angle integrated cross section and multiplying eq.(1) by $`(p_l^{eff}/p_l)^2`$. It is also shown that the Modified EMA (MEMA) works better than EMA for $`\nu _\mu `$ of low and high energies. In this approximation eq.(1) is multiplied by $`p_l^{eff}E_{eff}/p_lE_l`$. We use therefore this method in all our calculations of the $`(\nu _\mu ,\mu ^{})`$ cross sections. In the case of the $`(\nu _e,e^{})`$ process, the situation is somewhat more complicated. The Fermi function works only for very low energies, namely $`E_e10MeV`$ (where $`p_eR1`$, $`R`$ is the nuclear radius), whereas the EMA seems to be a good approximation for most energies of the outgoing electrons . As in ref., for $`\nu _e`$ , we treat Coulomb corrections by interpolating between the Fermi function at low electron energies and the EMA approximation at high lepton energies.

To get flux-averaged cross sections it is necessary to convolute (1) by the neutrino flux $`f(E_\nu )`$, that is

$$\sigma _f=_{E_0}^{\mathrm{}}𝑑E_\nu \sigma (E_\nu )f(E_\nu ),$$
(2)
$`E_0`$ being the threshold energy. The choice of $`f(E_\nu )`$ depends on the neutrino source and can be taken for example equal to the supernova neutrino energy spectrum given by transport codes or the neutrino fluxes produced by a beam dump.

The nuclear structure model used to evaluate the transition probabilities $`M_\beta `$ in (1) is the charge-exchange Random-Phase-Approximation (RPA). The details of the approach can be found in . The calculations we present have been obtained in a self-consistent approach: the HF single-particles energies and wavefunctions as well as the residual particle-hole interaction are derived from the same effective forces, namely the SIII  and SGII  Skyrme forces. We have found that the model configuration space used is large enough for the Ikeda and Fermi sum rule to be satisfied as well as the non energy-weigthed and energy-weighted sum rules for the forbidden transitions . The GT strength distribution we have obtained is peaked at $`19.2MeV`$, in agreement with the experimental value. This main peak exhausts about $`60\%`$ of the Ikeda sum rule. The IAS results at $`18.4MeV`$ and this value compares again well with the experimental finding ($`18.8MeV`$). Apart from these two resonances and the spin-dipole, the experimental knowledge about states of higher multipolarity is rather poor. The recent experiment of Ref.  shows that isovector monopole strength exists in $`{}_{}{}^{208}Bi`$ between $`30`$ and $`45MeV`$ and in the present calculation we find some strength in the same energy region.

In figs.1 and 2 we show the non flux-averaged $`{}_{}{}^{208}Pb(\nu _e,e^{})^{208}Bi`$ and $`{}_{}{}^{208}Pb(\nu _\mu ,\mu ^{})^{208}Bi`$ inclusive cross sections as a function of the neutrino energy, for a mesh of energies, namely $`\mathrm{\Delta }E=2.5MeV`$ for $`E_{\nu _e}`$ and $`\mathrm{\Delta }E=5.0MeV`$ for $`E_{\nu _\mu }`$. The dashed line in fig.1 shows the cross section obtained when only the Fermi function is used to include the Coulomb corrections. The results shown have been obtained with the SIII force, but we have found that with the SGII force we get quite similar results. All the multipolarities with $`J6`$ are included. We have checked that the contribution coming from $`J=7`$ is small. (Note that, for higher multipolarities, a mean field description, neglecting the particle-hole residual interaction, can be used to evaluate the transition probabilities (1)). In the calculations we present the axial vector coupling constant has been taken equal to 1.26. Note that the use of an effective $`g_a`$ to take into account the problem of the “missing” GT strength will reduce the reaction cross section by $`1015\%`$ as it was already discussed in .

Figure 3 shows the contribution of the different multipolarities to the total cross section (fig.1), for the impinging neutrino energies $`E_{\nu _e}=15,30,50MeV`$, which are characteristic average energies for supernova neutrinos. When $`E_{\nu _e}=15MeV`$ (fig.3, up), $`\sigma _{\nu _e}`$ is dominated by the allowed Gamow-Teller ($`J^\pi =1^+`$) transition. As the neutrino energy increases (fig.3, middle), the allowed IAS and other forbidden transitions start to contribute significantly. Finally, when $`E_{\nu _e}=50MeV`$ (fig.3, bottom), the GT and IAS transitions are not dominating at all, the cross section is being spread over many multipolarities. These results suggest that r-process nucleosynthesis calculations such as , which include neutrino-nuclei reactions, should take into account forbidden transitions. This may be even more important if $`\nu _\tau ,\nu _\mu \nu _e`$ oscillations occur, because in this case electron neutrino may have a higher average energy than it is usually expected from current supernovae models.

Let us now come to the comparison with other available calculations. Table 1 shows our flux-averaged cross sections, in comparison with those of refs.. The low-energy neutrino flux is given by a Fermi-Dirac spectrum

$$f(E_\nu )=\frac{1}{c(\alpha )T^3}\frac{E_\nu ^2}{exp\left[(E_\nu /T)\alpha \right]+1}$$
(3)
where $`T,\alpha `$ are fitted to numerical spectra and $`c(\alpha )`$ normalizes the spectrum to unit flux. The values of the parameters $`T`$ and $`\alpha `$ have been chosen to be able to compare our results with those of . As we can see from table 1, our predictions are in close agreement (the difference is at most $`2030\%`$) with . The results of  have been obtained in a CRPA approach. A variation of $`2030\%`$ is actually to be expected for calculations based on the same approach but using different parametrization (for example for single particle wavefunctions and effective particle-hole interaction), because of the sensitivity of the flux-averaged cross sections to the detailed strength distributions , as we will discuss further. On the contrary, our results and those of ref. present significant differences with those of , obtained using the allowed approximation and including the IAS, the GT and the first-forbidden contributions treated on the basis of the Goldhaber-Teller model.

We have checked that the differences do not come from the higher order forbidden transitions which are not included in the calculations of . The three calculations satisfy the same constraints, namely they reproduce the centroid of the resonances and satisfy the sum-rules.

We believe that the significant differences (by a factor of 2) with  may have two origins. The first possible origin might be the way the Coulomb corrections are treated. In , the Coulomb distortion of the outgoing electron wave function was taken into account by multiplying the cross section (1) by a Fermi function. In order to see the effect of using only the Fermi function instead of making an interpolation between the Fermi function and the EMA approximation, we have calculated the reaction cross sections using these two possible corrections. As figure 4 shows, the two cross sections have a quite different behaviour as a function of the neutrino energy so that this difference on the flux-averaged cross section may vary according to the particular neutrino flux considered. To get a quantitative idea of the variation, we have calculated the flux-averaged cross sections by convoluting the two curves of fig.4 with (3). If we use the Fermi function only, the reaction cross sections increase, on average, by $`50\%`$.

The second possible origin of the discrepancies between our work,  and  might be the sensitivity of the flux-averaged cross sections to the detailed strength distributions in $`{}_{}{}^{208}Bi`$. In fact, it has already been discussed in , that for low-energy neutrinos, the flux-averaged cross sections are very sensitive to the energy of the excited states in the final nucleus. The reason is twofold. First, due to the small electron mass, the non flux-averaged cross section (1) scales as the square of energy of the states. Second, the energy dependence of the neutrino flux may emphasize differences in the non flux-averaged cross sections due to variations in the energy of the states. As it was discussed in , these two effects may modify the flux-averaged reaction cross sections by $`2030\%`$.

To complete our comparison with the calculations of , we have calculated two more flux-averaged cross sections, using the neutrino fluxes of both $`\nu _\mu `$ coming from the DIF of $`\pi ^+`$ and $`\nu _e`$ coming from the DAR of $`\mu ^+`$. The neutrino fluxes $`f(E_\nu )`$ were taken from . These neutrino fluxes have been used in the recent experiments $`\nu _\mu \nu _e`$ , $`\overline{\nu }_\mu \overline{\nu }_e`$  or $`\nu _\mu \nu _x`$  performed by the LSND and KARMEN collaborations. The $`DAR(\nu _e,e^{})`$ cross section calculated is $`\sigma _{DAR}=44.3910^{40}cm^2`$ which is very close to $`36.210^{40}cm^2`$ obtained in . On the contrary, our $`DIF(\nu _\mu ,\mu ^{})`$ is $`\sigma _{DIF}=399.210^{40}cm^2`$; whereas the one of  is $`11510^{40}cm^2`$. We believe that some of the disagreement may come from differences in the strength distributions of the high order (higher than 2) forbidden transitions. In fact, contrary to the reactions of neutrinos on light nuclei such as carbon, where these states contribute only by $`20\%`$ to the total DIF cross section, their contribution represents $`65\%`$ of the total cross section when the nucleus is as heavy as lead.

Let us finally discuss the two possible neutrino oscillation $`\nu _\mu ,\nu _\tau \nu _e`$ signals based on the spectroscopic properties of $`{}_{}{}^{208}Bi`$ excited in the CC reaction that have been proposed recently. In , it was shown that the 2-neutron events associated with the deexcitation of $`{}_{}{}^{208}Bi`$ are very sensitive to the mean electron neutrino energy. This signal relies on the fact that most of the IAS, GT and first-forbidden strength distributions are above the $`2n`$ emission threshold ($`14.98MeV`$) in $`{}_{}{}^{208}Bi`$. Our results show that not only the allowed and spin-dipole strengths are above this threshold, but also a fraction of the strength distributions associated with other forbidden transitions (fig.3) will contribute to the $`2n`$ decay. All the arguments given in  are based on the statistical calculations of $`1n`$ and $`2n`$ decays. The direct $`1n`$ emission represents about $`50\%`$ of the total width in the case of the IAS, and $`510\%`$ in the case of the GT .

In , it was pointed out that the energy distribution of the neutrons in the $`1n`$ events should form a peak at low energy, more or less pronounced according to the occurence or absence of oscillations. This peak comes from the GT strength distribution at around $`7.6MeV`$ which is located above the $`1n`$ threshold emission at $`6.9MeV`$. Our GT distribution also shows a peak at around $`7.5MeV`$. We have checked that its location is not sensitive to the choice of the effective forces used. Still one should be careful about conclusions, because predictions of different models about the energy location and strength of that peak are at variance.

In summary, we have presented the non flux-averaged $`{}_{}{}^{208}Pb(\nu _e,e^{})^{208}Bi`$ and $`{}_{}{}^{208}Pb(\nu _\mu ,\mu ^{})^{208}Bi`$ reaction cross sections, calculated in a self-consistent charge-exchange Random-Phase-Approximation with Skyrme effective forces. These predictions can be employed for very different purposes, such as for the interpretation of the recent experiments on neutrino oscillations performed by the LSND collaboration (where reactions induced by neutrinos on lead contribute significantly to the background) and to evaluate the feasibility of future projects in which lead should be used as detector for supernova neutrinos. We have emphasized that forbidden transitions contribute significantly to the neutrino-nuclei reaction cross sections even at the “astrophysical neutrino energies” and they should be included in present r-process nucleosynthesis calculations. We have discussed the present status on the theoretical predictions on the reaction cross sections for the $`\nu _e`$ having typical energies from present models on supernovae. If on one hand our calculations agree with those of ref., which are also based on RPA; on the other hand, they both significantly disagree with those of ref.. We point out that the origin of the discrepancy might be mainly the different treatment of Coulomb corrections, but also the sensitivity of the reaction cross sections to the detailed energy spectrum of the final nucleus. We have also compared our flux-averaged reaction cross sections with $`\nu _\mu `$ coming from the DIF of $`\pi ^+`$ and with $`\nu _e`$ coming from the DAR of $`\mu ^+`$, with the ones of . As expected, the DAR cross sections are very close. On the contrary our DIF cross section differs significantly from the one of . We have pointed out that the two predictions may differ because of differences in the strength distributions of forbidden transitions of high multipolarity which represent the main contribution in reactions of neutrinos on nuclei as heavy as lead. Finally, we have discussed our results in relation with recently proposed signals to measure supernova neutrino oscillations.

This work was supported by the US-Israel Binational Science Foundation.


format: text

text:
# Fractional and unquantized dc voltage generation in THz-driven semiconductor superlattices

## Abstract

We consider the spontaneous creation of a dc voltage across a strongly coupled semiconductor superlattice subjected to THz radiation. We show that the dc voltage may be approximately proportional either to an integer or to a half-integer multiple of the frequency of the applied ac field, depending on the ratio of the characteristic scattering rates of conducting electrons. For the case of an ac field frequency less than the characteristic scattering rates, we demonstrate the generation of an unquantized dc voltage.

The theoretical analysis of nonlinear transport properties of strongly-coupled semiconductor superlattices (SSLs) irradiated by a high-frequency electric field began already in the mid 1970’s . Recently, experimental progress in creating powerful sources of THz radiation, the development of a coupling technique , and improvement in the fabrication technology of microstructures leading to very high carrier mobility  have stimulated many new theoretical investigations of this long-studied problem. Among these works have been studies of strongly nonlinear effects including multistability , short pulse generation , chaos , and spontaneous generation of a dc voltage in a purely ac-driven SSL .

In this paper we investigate the last effect in further detail and discuss the appearance of new dc voltage states, which are a generalization of the integer dc voltage states in SSLs previously described . These states are related to the complex dynamics of miniband electrons in an SSL and the formation of the Wannier-Stark ladder in a purely ac driven SSL. Two distinct mechanisms are known for the spontaneous generation of a dc bias in purely ac-driven SSLs . Both these nonlinear mechanisms work in SSLs with a high mobility and a relatively high level of doping, when the effects of a self-consistent electric field generated by an electron’s motion become significant. The first mechanism arises if the ac field frequency, $`\omega `$, is much greater than the plasma frequency, $`\omega _{pl}`$, and is related  to an instability caused by absolute negative conductivity in the ac-driven SSL . The effect has been attributed  to the phenomena of dynamical localization of electrons  and miniband collapse in a collisionless SSL . It was reported that the generated dc bias is such that the “induced Bloch frequency” $`\omega _B=eaE_{dc}/\mathrm{}`$ ($`E_{dc}`$ is the spontaneously generated dc electric field and $`a`$ is the SSL period) is approximately equal to the ac field frequency $`\omega `$ .

The second mechanism is responsible for dc bias generation when the ac field frequency is near the plasma resonance, $`\omega \omega _{pl}`$; it arises for a smaller ac field strength than is required in the previous situation . The instability responsible for the dc bias generation in SSLs with strong enough electron scattering also results in chaotic motion in the case of small scattering rates or for a collisionless SSL . The creation of a dc bias may be qualitatively explained  and classified  using the semiclassical theory of wave-mixing in SSLs .

In this paper we re-examine the problem of spontaneous dc voltage generation in an SSL subjected to a THz electric field. We show that, depending on the relative values of the scattering rates and the ac field frequency, a variety of different dc voltage states can exist, including both integer and half-integer quantized states, for which the induced Bloch frequency is approximately an integer or half-integer multiple of the ac field frequency, and completely unquantized states. In particular, if the electron velocity relaxation rate, $`\gamma _v`$, is sufficiently different from the electron energy relaxation rate, $`\gamma _\epsilon `$, and $`\omega \gamma _v`$, we find integer states with $`\omega _Bn\omega `$ ($`n=\pm 1,\pm 2,\mathrm{}`$); while for $`\gamma _v=\gamma _\epsilon `$, the states are close to the half-integer states $`\omega _Bn\omega /2`$. In contrast, in the case of low-frequency driving or high damping, $`\omega <(\gamma _v,\gamma _\epsilon )`$, the dc voltage states are unquantized.

We study electron transport through a single miniband, spatially homogeneous SSL with period $`a`$ and miniband width $`\mathrm{\Delta }`$, which is subjected to an ac electric field $`E(t)=E_0\mathrm{cos}\omega t`$ along the SSL axis. For the tight-binding energy-quasimomentum dispersion relation $`\epsilon (k)=(\mathrm{\Delta }/2)\left[1\mathrm{cos}(ka)\right]`$ ($`k`$ is the electron wave vector along the axis of SSL), the dynamics of electrons is described by the superlattice balance equations

$`\dot{v}`$ $`=`$ $`uw\gamma _vv,`$ (1)
$`\dot{w}`$ $`=`$ $`uv\gamma _\epsilon (ww_{eq}),`$ (2)
$`\dot{u}`$ $`=`$ $`\omega _{pl}^2v\alpha u+I_{ext}(t),`$ (3)
where $`v=m_0\overline{V}a/\mathrm{}`$, $`w=(\overline{\epsilon }\mathrm{\Delta }/2)(\mathrm{\Delta }/2)^1`$ and $`w_{eq}`$ are a scaled electron velocity, a scaled electron energy, and an equilibrium value of scaled electron energy, respectively, and $`m_0=(2\mathrm{}^2)/(\mathrm{\Delta }a^2)`$ is the effective mass at the bottom of miniband. The scaled variables $`v(t)`$ and $`w(t)`$ are proportional to the variables $`\overline{V}(t)`$ and $`\overline{\epsilon }(t)`$, which are the electron velocity and energy averaged over the time-dependent distribution function satisfying the Boltzmann equation. The lower (upper) edge of the miniband corresponds to $`w=1`$ ($`w=+1`$), and the value of $`w_{eq}`$ is a function of the lattice temperature (for a thermal equilibrium $`w_{eq}<0`$). The variable $`u(t)`$ is related to the electric field inside the SSL $`E(t)`$ as $`u=eaE/\mathrm{}`$. In deriving Eqs. (2), we assumed that the electrical properties of an SSL of total length $`l`$ and cross-section $`S`$ can be modeled by an equivalent high-quality circuit  which consists of a capacitor $`C=(ϵ_0S)/(4\pi l)`$ ($`ϵ_0`$ is the average dielectric constant for the SSL) driven by an ac current of the form $`I_{ext}=\omega _s\omega \mathrm{sin}\omega t`$, where $`\omega _s=eE_0a/\mathrm{}`$, in our scaled units.

The first Eq. of set (2) describes an acceleration of electrons under the action of electric field $`E(t)`$ and their slowing down caused by an effective friction due to scattering. In original dimensional variables, the term $`uw`$ in r.h.s. of the first Eq. looks like $`eE/m(\overline{\epsilon })`$, where $`m(\overline{\epsilon })=m_0/(12\overline{\epsilon }/\mathrm{\Delta })`$ is energy-dependent effective mass of the electrons in SSL’s miniband. The second Eq. describes a balance of electron’s energy gain under the action of electric field and energy loss due to scattering. Finally, the third equation describes a balance of diffusive, external and displacement currents in the SSL. The degree of nonlinearity in Eqs. (2) is controlled by the value of miniband plasma frequency, $`\omega _{pl}=\left(4\pi e^2N/m_0ϵ_0\right)^{1/2}`$, which is a function of the electron doping density $`N`$, while the parameter $`\alpha `$ determines the quality of the effective circuit ($`\alpha \omega _{pl},\omega `$).

The relaxation processes for miniband electrons are characterized by an average energy scattering rate, $`\gamma _\epsilon `$, as well as by an average velocity scattering rate, $`\gamma _v=\gamma _\epsilon +\gamma _{el}`$, where $`\gamma _{el}`$ is an average rate of elastic collisions . The scattering rates, $`\gamma _v`$ and $`\gamma _\epsilon `$, can have different values depending on the material, the doping density, the temperature, etc. In particular, for microstructures with modulation doping , the ionized impurities are spatially separated from electrons, which greatly reduces the elastic scattering rate $`\gamma _{el}`$ so that $`\gamma _v\gamma _\epsilon `$. In contrast, for many vertical SSLs operating at room temperature, the scattering rate for electron velocity is about of order of magnitude greater than the characteristic scattering rate of electron energy, $`\gamma _v/\gamma _\epsilon 10`$ .

We solve the nonlinear balance Eqs. (2) numerically for the initial conditions $`v(0)=0`$, $`w(0)=w_{eq}=1`$, with the circuit damping rate $`\alpha /\omega _{pl}=0.01`$, and for two typical sets of relaxation constants: (i) $`\gamma _v/\gamma _\epsilon =10`$, $`\gamma _\epsilon /\omega _{pl}=0.01`$, and (ii) $`\gamma _v/\gamma _\epsilon =1`$, $`\gamma _\epsilon /\omega _{pl}=0.1`$. After removing the transients, we calculate the time-average value $`u`$, which gives the value of the Bloch frequency determined by the spontaneously generated dc bias $`E_{dc}`$: $`\omega _BeaE_{dc}/\mathrm{}=u`$. Figs. 1 and 2 present the results of computations of $`u`$ for 201 value of $`\omega _s`$ equally distributed in the range $`0\omega _s2\omega _{pl}`$ for each driving frequency.

For the first set of relaxation rates and for $`\omega >2\gamma _v`$, plateaus of integer-quantized states are clearly observable: $`\omega _Bn\omega `$, with $`n=\pm 1,2,3,4`$ in Fig.1. However, at low frequencies, $`\omega <2\gamma _v`$, instead of steps there is a region of unquantized dc states. The dependence of the induced Bloch frequency $`u`$ on the ac field frequency $`\omega `$ for the second set of damping parameters is presented in Fig.2. Here, qualitative differences from Fig. 1 appear, including 1) the existence of half-integer states, specifically states with $`n\pm 1/2`$, and 2) the nonzero width of the $`n=0`$ plateau. One can expect the width of a plateau to be equal to the scattering rate, $`0.1\omega _{pl}`$ in this case. As a result, we have plotted the lines $`n\pm 0.1\omega _{pl}/\omega `$ for $`n=0,\pm 1,\pm 2,\pm 3,4,5`$ in Fig.2. We found that, indeed, most points for a given plateau lie in the region demarcated by the two lines with same $`n`$. In contrast to Fig. 1, there are states with $`u/\omega <0.1`$, but $`u0`$. Also, some points fall very near the line $`u=\pm 0.5\omega `$. As an example we refer to the solution of Eqs. (2) for $`\omega /\omega _{pl}=0.6`$ and $`\omega _s/\omega _{pl}=0.6`$ (other parameters are same as in Fig.2), which in the phase space corresponds to a symmetry-broken limit cycle*** This is a limit cycle whose projection on the $`vu`$ plane is not symmetric about the origin, in contrast to a symmetric limit cycle (see ). with $`u=0.289\omega _{pl}`$.

If $`\omega /\gamma `$ is less or order of unity but $`I_{ext}/\omega _{pl}^2`$ is large enough and $`\gamma _v=\gamma _\epsilon \gamma `$, then the existence of nonquantized and half-integer dc voltage states can be demonstrated analytically. In these limits, the motion on an attractor of dynamical system (2) is governed by the following pendulum equation

$$\ddot{\theta }+\gamma \dot{\theta }+(w_{eq})\omega _{pl}^2\mathrm{sin}\theta =2I_{ext}(t),$$
(4)
where we made the substitutions: $`v=(w_{eq}/2)\mathrm{sin}\theta `$, $`w=(w_{eq}/2)(1+\mathrm{cos}\theta )`$ and

$$u=\dot{\theta }/2+\gamma \mathrm{tan}(\theta /2).$$
(5)
The dc voltage $`u`$ can be obtained as a result of averaging over time of r.h.s. of Eq. (5). The equation (4) looks same as a motion equation of the well-known Stewart-McCumber model from the theory of ac-driven Josephson juctions . However, the principal difference from this model also exists. In the Stewart-McCumber model the voltage across the Josephson junction is proportional to the velocity of pendulum , while in our case the voltage is a function of both velocity and co-ordinate (see Eq. (5)). This difference plays a principal role in the explanation of existence in our case of both unquantized and only approximately quantized dc voltage states, which are absent in the Josephson junction model.

The damped and driven pendulum (Eq. (4)) have two distinct types of attractors with regular dynamics: rotating and oscillating . Majority of rotating states are phase-locked, i.e., $`\dot{\theta }=(n/l)\omega `$ ($`n`$ and $`l`$ are integer numbers); they are observable mainly at $`\omega \gamma `$ . As it is evident from Eq. (5), such phase-locked rotational pendulum states are responsible for the generation of half-integer dc voltage states in SSL: $`u(n/2l)\omega `$. Note that these corrections to these quantized states arise from the contribution of $`\theta `$-dependent term in Eq. (5). The corrections can give the dependence of $`u`$ on ac current amplitude, $`\omega _s`$, its frequency, $`\omega `$, and scattering constant, $`\gamma `$. However, it can be shown that a relative contribution into these corrections is controlled by the parameter $`\gamma /\omega _{pl}`$ . Therefore, the dependence of approximately quantized voltage on ac field amplitude and frequency is weak until $`\gamma /\omega _{pl}1`$.

For oscillating attractors, the equality $`\dot{\theta }=0`$ is always valid . If additionally stationary pendulum oscillations are symmetric, i.e. $`\theta =0`$, then the dc voltage can not be generated, $`u=0`$. However, the pendulum (Eq. (4)) can demonstrate symmtry-broken oscillations, for which $`\theta 0`$ . The symmetry-broken oscillations mainly exist at the low frequencies of external force; they can survive at $`\omega <\gamma `$ even for a strong damping, when $`\gamma /\omega _{pl}1`$ . As it is evident from Eq. (5), just the symmetry-broken oscillations corresponding to $`\dot{\theta }=0`$ and $`\theta 0`$, are responsible for the generation of the unquantized dc voltage in SSL.

Unquntized dc voltage states exist in the low frequency region for different ratious of the scattering constants. In order to understand the transition from quantized to unquantized dc bias states as the scattering rates increase while maintaining $`\gamma _v\gamma _\epsilon `$, we present in Fig.3 the dependence of $`u`$ on $`\omega `$ for strong damping. As it is evident from this figure, the strong damping destroys all quantized states; dc bias generation persists only for some unquantized states.

It is instructive to consider the mechanism of unquantized dc bias generation in the terms of energy levels structure. We offer the following qualitative explanation of the underlying physics: when the electron scattering rates are sufficiently small and the amplitude of the ac field is large enough, the SSL spontaneously creates a Wannier-Stark ladder with the spacing, $`\omega _B`$, that makes multiphoton absorption of ac field most effective, i.e. $`\omega _Bn\omega `$. In the pendulum analogy this means the appearance of the pendulum rotations, which frequency of rotation is proportional to the generated dc bias. However, for larger damping, the rotations are ceased and therefore the symmetry broken oscillations remain, that means that only a small bias can be generated in the SSL. Hence the spacing of the Wannier-Stark ladder states is less than the ac field frequency except for very low frequencies $`\omega \gamma `$. In this case, the broadening of the self-organized ladder levels is quite comparable with their spacing; it is therefore practically impossible to achieve quantized values of the voltage (or a rotation of the pendulum). In a real space picture such a situation should correspond to an appearence of some kind of semi-localized Wannier-Stark wave functions; this structure of the ladder is reminiscent of the wave function of biased short superlattice , where the spacing between energy levels depends non-linearly on the bias voltage . In our case, a dc bias is created by ac field, thus we have a weak dependence of ladder spacing on ac field strength via the self-induced bias.

We have performed a systematic numerical study of the positions and widths of different plateaus and of the unquantized states for many values of the driving amplitude and frequency and several different initial conditions at different damping levels ; the results are in a qualitative agreement with situation described above. We now make several remarks on the results of this search. First, we found no indication of the noninteger (fractional) dc states for $`\gamma _v\gamma _\epsilon `$. However, for $`\gamma _v=\gamma _\epsilon `$, we find that the $`1/2`$-dc-states are quite common. Moreover, for weak enough damping, we additionally saw a few dc states that are very close to fractional states of the form $`n/k`$ with $`n`$ being an integer and $`k`$ always being an even integer. Such dc states are formed by the symmetry-broken limit cycles with large even periods. As examples, we refer to $`3/2`$-states formed by period-12 and period-24 limit cycles, as well as to the $`7/6`$-state formed by a period-12 limit cycle; both occur for $`\gamma _v=\gamma _\epsilon =0.05\omega _{pl}`$ and $`\alpha /\omega _{pl}=0.01`$. We should note, however, that for weak damping, chaotic behaviour is quite typical  and that both the stable limit cycles and the asymmetric chaotic attractors, which are responsible for the generation of stable quantized dc voltage states in the SSL, occupy only a small amount of parameter space of the system .

Importantly, it appears possible to achieve the parameter values used in our simulations in the conditions of experiments. To begin with we refer to the experiments with a heavily doped SSL ($`N=8\times 10^{16}`$ $`\mathrm{cm}^3`$) having at room temperatures scattering constants $`\gamma _v10^{13}`$ $`\mathrm{s}^1`$ and $`\gamma _\epsilon 0.1\gamma _v`$ . Using the values $`a=4.8`$ nm, $`\mathrm{\Delta }50`$ meV and $`ϵ_013`$ , we get $`\omega _{pl}=1.2\times 10^{13}`$ rad/s and $`\gamma /\omega _{pl}0.84`$ (compare, e.g., with data of Fig.3).

Longer scattering time provide SSLs based on the cleaved edge overgrown technique: $`\gamma ^13`$ ps at low temperatures . In conditions of experiment , for $`a=10`$ nm, $`\mathrm{\Delta }20`$ meV, $`N=N_sb3\times 10^{15}`$ $`\mathrm{cm}^3`$ ($`N_s=3\times 10^{11}`$ $`\mathrm{cm}^2`$ is an electron areal density and $`b10^4`$ cm is sample’s thickness), the miniband plasma frequency is $`\omega _{pl}=3.1\times 10^{12}`$ rad/s providing $`\gamma /\omega _{pl}0.1`$ (cf. with data of Fig. 2). In our numerical simulations we used the ac field strengths satisfying $`\omega _s2\omega _{pl}`$, what in the physical units correspond to the realistic values $`E_05`$ kV/cm. Necessary for the observation of our effects ac field frequencies belong to the THz range. Finally, we should also note that for all listed cases, the standard condition of validity of the single-miniband approximation, $`\mathrm{\Delta }\mathrm{}\omega _s`$ , is well satisfied because $`\omega _s\mathrm{}/\mathrm{\Delta }0.1`$ .

In summary, we have shown that a semiconductor superlattice irradiated by a high-frequency electric field can spontaneously generate a dc bias, which can be quantized in approximately integer or particular fractional ratios of the driving frequency, or completely unquantized. In this respect, the effects in semiconductor superlattices are no less rich than their counterparts in Josephson junctions subjected to a microwave field, where the exactly integer and the exactly fractional dc voltage states (“phase-locked states”) are known .

We thank Anatoly Ignatov, Pekka Pietiläinen and Karl Renk for discussions. This research was partially supported by the Academy of Finland (grant 163358) and NorFA.


text:
# Disk-Assisted Spindown of Young Radio Pulsars

## 1 Introduction

The magnetic dipole model for radio pulsars has proven extremely successful in explaining the observational properties of these objects (Pacini 1967; 1968; Gunn & Ostriker 1969; Gold 1968; Goldreich & Julian 1969; see, e.g., Shapiro & Teukolsky  for an overview; Manchester & Taylor  for a detailed review). A longstanding problem of the theory, however, has been that the predicted braking index (see §2.1 for definition) is $`n=3`$, while measurements of braking indices in five young radio pulsars are all less than 3 (see §3 for the most recent measured values). Various effects have been proposed to explain this discrepancy, such as magnetic axis wandering (Macy 1974), non-standard vacuum dipole theory (Melatos 1997), and non-dipolar field structure or field strength secular evolution (Blandford, Applegate & Hernquist 1983; Blandford & Romani 1988).

In this Letter, we explore the consequences of a recent proposal by Chatterjee, Hernquist & Narayan (2000) according to which young radio pulsars could be surrounded by remnant disks produced by supernova fallback (see also Michel & Dessler 1981, 1983; Michel 1988; Yusifov et al. 1995; Alpar 1999, 2000). The disk may indeed be able to contribute to the neutron star spindown even in the radio pulsar phase, provided that matter does not fall through the magnetosphere inside the light cylinder radius. We note that Marsden, Lingenfelter & Rothschild (2001a,b) have recently considered some of the consequences of the presence of such a disk around a radio pulsar, with a special emphasis on explaining discrepancies between a pulsar’s true age and its timing age (see also Gvaramadze  for the case of an interaction with dense circumstellar material). The outline of this Letter is as follows. In §2, we derive a theory for pulsar spindown in the presence of a fallback disk and the associated propeller torques. We apply the theory to the Crab pulsar and other young radio pulsars in §3.

## 2 Pulsar Spindown Model

We certainly lack a proper understanding of the magnetospheric interaction between a radio pulsar and a surrounding disk. However, simple arguments can be used to derive a reasonable rate at which the pulsar is spun down by this interaction. Implicit to our model is the assumption that the disk can survive the high-energy environment created by the intense relativistic pulsar wind. For the pulsar mechanism to operate, we require the disk to be located beyond the light cylinder radius, at $`R_{lc}=c/\mathrm{\Omega }`$ (where $`\mathrm{\Omega }`$ is the neutron star angular speed and $`c`$ is the speed of light). On the other hand, the magnetic field lines are open beyond $`R_{lc}`$ and the field structure is no longer dipolar. This is the main motivation for our assumption that the disk extends down to $`R_{lc}`$. At $`R_{lc}`$, the inflowing material encounters a dipolar field structure rotating at nearly the speed of light: an efficient propeller effect (Illarionov & Sunyaev 1975) may be expected to operate there, spinning down the neutron star and preventing mass inside the pulsar magnetosphere at the same time.111We note that, if the disk were truncated beyond $`R_{lc}`$ (e.g. because of a large energy density right outside the light cylinder), the torque on the neutron star would likely be largely reduced or entirely suppressed.

Following Menou et al. (1999; see also Daumerie 1996), we adopt an efficient propeller torque

$$\dot{J}I\dot{\mathrm{\Omega }}=2\dot{M}R_m^2\mathrm{\Omega },$$
(1)
where $`I`$ is the moment of inertia of the neutron star, a dot denotes a time derivative and $`\dot{M}`$ is the disk mass inflow rate at the magnetospheric radius, $`R_m`$. This expression for the torque assumes that the material flung away by the propeller effect has been accelerated to an angular speed corresponding to that of the star. In the pulsar case, we effectively take $`R_m=R_{lc}`$ and the material flung away only approaches the speed of the light. Under these simple assumptions, the torque for this mildly relativistic propeller is given by

$$\dot{J}=2\dot{M}c^2\mathrm{\Omega }^1,$$
(2)
which is essentially what one expects from dimensional arguments assuming the torque acts at $`R_{lc}`$. A more efficient, relativistic propeller is energetically possible; it would correspond to a reduced value of the disk mass inflow rate for a given value of the torque.

The pulsar loses rotational energy at a rate:

$$\dot{E}I\mathrm{\Omega }\dot{\mathrm{\Omega }}=\beta \mathrm{\Omega }^4\gamma ,$$
(3)
where the first term on the RHS corresponds to magnetic dipole radiation losses and the second term to propeller losses. We have

$$\gamma =2\dot{M}c^22\times 10^{37}\mathrm{erg}\mathrm{s}^1\left(\frac{\dot{M}}{10^{16}\mathrm{g}\mathrm{s}^1}\right),$$
(4)
and

$$\beta =\frac{B_p^2\mathrm{sin}^2\alpha R_{NS}^6}{6c^3},$$
(5)
where $`R_{NS}`$ is the neutron star radius and $`B_p`$ is the magnetic field strength at the magnetic pole, inclined at an angle $`\alpha `$ from the rotation axis (see, e.g., Shapiro & Teukolsky 1983).

In what follows, we assume that Eq. (3) is valid at all times. We note that this assumption is not entirely consistent with the evolutionary scenario of Chatterjee et al. (2000), in which a neutron star decouples completely from its surrounding fallback disk when it enters the radio pulsar phase. We expect this discrepancy to be minor given that magnetic dipole radiation still dominates the losses in the present model. More importantly, a young neutron star is not initially in the radio pulsar phase in the Chatterjee et al. model, but in a “dim propeller” phase (during which the magnetospheric radius has not yet reached the light cylinder radius). This phase can last up to several hundreds of years in the Chatterjee et al. model and it is not captured by our Eq. (3). This discrepancy suggests that the results derived below for the integrated neutron star spin history should be interpreted with caution (instantaneous quantities such as the braking index are not directly affected, however).

### 2.1 Steady mass inflow rate

From Eq. (3), we find that the braking index is

$$n\frac{\ddot{\mathrm{\Omega }}\mathrm{\Omega }}{\dot{\mathrm{\Omega }}^2}=3\frac{1}{\frac{1}{4}+\frac{\beta }{4\gamma }\mathrm{\Omega }^4},$$
(6)
if $`\gamma `$ (i.e. the disk mass inflow rate) is constant with time. The breaking index is therefore less than $`3`$, the value corresponding to pure magnetic dipole radiation losses (recovered in the limit $`\gamma 0`$).

Following the usual procedure, Eq. (3) can be separated and integrated to constrain the neutron star spin history. We rewrite Eq. (3) as

$$\frac{d\mathrm{\Omega }^2}{1+\frac{\beta }{\gamma }\mathrm{\Omega }^4}=\frac{2\gamma }{I}dt,$$
(7)
which can be integrated (assuming $`\gamma =`$ constant) by noting that the LHS is of the form $`dy/(1+y^2)=d(\mathrm{arctan}(y))`$. This yields, at any time $`t`$ after the system’s birth:

$$t=\frac{I}{2(\gamma \beta )^{1/2}}\left[\mathrm{arctan}\left[\left(\frac{\beta }{\gamma }\right)^{1/2}\mathrm{\Omega }_i^2\right]\mathrm{arctan}\left[\left(\frac{\beta }{\gamma }\right)^{1/2}\mathrm{\Omega }^2(t)\right]\right],$$
(8)
where $`\mathrm{\Omega }(t)`$ is the angular rotation speed at time $`t`$ and $`\mathrm{\Omega }_i`$ is the angular rotation speed at birth. In the limit of very large $`\mathrm{\Omega }_i`$, the left $`\mathrm{arctan}`$ is simply equal to $`\pi /2`$. There is no simple relation between the pulsar timing age and true age in this solution, contrary to the case of pure magnetic dipole radiation losses (Shapiro & Teukolsky 1983; Gunn & Ostriker 1969).

### 2.2 Time-dependent mass inflow rate

We generalize here the above calculations to the case of a mass inflow rate decreasing with time, as expected in the fallback disk scenario (Chatterjee et al. 2000; Menou, Perna & Hernquist 2001). Allowing the mass inflow rate to vary with time, the braking index becomes:

$$n=3\frac{1\frac{1}{4}\frac{\dot{\gamma }}{\gamma }\frac{\mathrm{\Omega }}{\dot{\mathrm{\Omega }}}}{\frac{1}{4}+\frac{\beta }{4\gamma }\mathrm{\Omega }^4},$$
(9)
where both $`\dot{\gamma }`$ and $`\dot{\mathrm{\Omega }}`$ are negative (assuming a decreasing mass inflow rate).222We note that this expression for the braking index can potentially yield values of $`n>3`$. The quantity $`\mathrm{\Omega }/\dot{\mathrm{\Omega }}`$ is twice the usual pulsar timing age ($`P/2\dot{P}`$), while the quantity $`\alpha \gamma /\dot{\gamma }`$ is actually the true age of the system if, as expected, the mass inflow rate obeys a power law evolution in time of the form $`\gamma (t)t^\alpha `$.

While the braking index and the second deceleration parameter (defined in §3.1) are instantaneous quantities which can easily be computed for any power-law solution (using, if known, the values of the true age and the timing age as explained above), the spin history cannot be computed analytically for an arbitrary power law index $`\alpha `$, but only in some special cases.

For simplicity, here we restrict our discussion to solutions for the neutron star spin history with a time-dependent mass inflow rate satisfying $`\gamma (t)=\gamma _0(t/t_0)^2`$. Let $`w=\mathrm{\hspace{0.17em}1}/\mathrm{\Omega }^2t`$, $`a=2\beta /I`$ and $`b=2\gamma _0t_0^2/I`$, then the differential equation (3) becomes

$$w+t\frac{dw}{dt}=a+bw^2$$
(10)
This can be separated and integrated:

$$\frac{dw}{bw^2w+a}\mathrm{ln}t=\mathrm{constant}$$
(11)
If $`4ab<1`$ then the integration yields

$$\frac{1}{c}\mathrm{ln}\left|\frac{2bw\mathrm{\hspace{0.17em}1}c}{2bw\mathrm{\hspace{0.17em}1}+c}\right|\mathrm{ln}t=\mathrm{constant}$$
(12)
where $`c=\sqrt{1\mathrm{\hspace{0.17em}4}ab}`$. The constant is fixed by the value of the variable $`w`$ at some initial time $`t_i`$. The solution for $`\mathrm{\Omega }`$ can then be easily derived:

$$\mathrm{\Omega }^2=\frac{2b}{t}\frac{(2bw_i1+c)(2bw_i1c)(t/t_i)^c}{(1+c)(2bw_i1+c)(1c)(2bw_i1c)(t/t_i)^c}.$$
(13)
At later times, $`tt_i`$, and the solution can be well approximated by

$$\mathrm{\Omega }^2\frac{1+c}{2at}.$$
(14)
Evaluating the above equation at the present time $`t_0`$ yields

$$t_0=\frac{I}{2\beta \mathrm{\Omega }_0^2\left(1+\frac{\gamma _0}{\beta \mathrm{\Omega }_0^4}\right)}.$$
(15)
Note that the factor $`I/(2\beta \mathrm{\Omega }_0^2)`$ is the pulsar timing age (in the case of pure magnetic dipole radiation losses). Although this solution corresponds to rather simple evolution laws for $`\mathrm{\Omega }`$ and $`\gamma `$, it is useful as an illustration of the effect of allowing for a time-varying mass inflow rate in the equations (as compared to the steady case considered in §2.1). Given that a fallback disk should have a mass inflow rate decreasing approximately as $`t^{1.2}`$ (Cannizzo et al. 1990; Chatterjee et al. 2000; Menou et al. 2001), the steady and simple time-varying cases considered here may be expected to bracket the real case ($`\gamma t^{1.2}`$), for which a numerical integration is required in full generality.

## 3 Applications

### 3.1 The Crab Pulsar

For the Crab pulsar (PSR B0531+21), we adopt $`I10^{45}`$ g cm-2 and $`\beta =3.85\times 10^{29}`$ erg s3. These values are derived from measurements of $`\mathrm{\Omega }`$, $`\dot{\mathrm{\Omega }}`$, the estimated value of $`\dot{E}5\times 10^{38}`$ erg s-1 and the assumption of pure magnetic dipole radiation losses (see Shapiro & Teukolsky 1983; Taylor, Manchester & Lyne 1993). The propeller losses are comparatively small in the applications considered below and for the Crab pulsar:

$$\frac{\beta \mathrm{\Omega }^4}{\gamma }25\left(\frac{\dot{M}}{10^{16}\mathrm{g}\mathrm{s}^1}\right)^1.$$
(16)
We recall here that, in 1972, the true age of the Crab pulsar was 918 yrs, while its timing age was measured to be 1243 yrs (which is half the characteristic age as defined by Shapiro & Teukolsky 1983). The braking index of the Crab pulsar is $`n=2.509`$ (e.g. Lyne et al. 1988).

In the case of steady mass inflow in the fallback disk outside the light cylinder, we find using Eq. (6) that the braking index is $`n=`$ 2.85, 2.57 and 1.86 for the mass inflow rates $`\dot{M}=10^{16}`$, $`3\times 10^{16}`$ and $`10^{17}`$ g s-1, respectively. For the same mass inflow rates, Eq. (8) predicts a system age of $`t_0=`$ 1130, 1101 and 1018 yrs, respectively (in the limit of very large initial angular speed $`\mathrm{\Omega }_i`$).

In the case of a mass inflow rate decreasing with time as $`t^{1.2}`$, we find using Eq. (9) that the braking index is $`n=`$ 2.97, 2.93 and 2.78 for the (current) mass inflow rates $`\dot{M}=10^{16}`$, $`3\times 10^{16}`$ and $`10^{17}`$ g s-1, respectively. The age cannot be simply computed with the power law index $`\alpha =1.2`$. However, it will be bracketed by the values obtained for steady mass inflow above and the $`\alpha =2`$ case, for which Eq. (15) holds. For the same (current) mass inflow rates as above, Eq. (15) predicts a system age of $`t_0=`$ 1099, 1020 and 817 yrs, respectively.

In each of the cases above (steady mass inflow and $`\gamma t^2`$), it is possible to find a mass inflow rate for which both the Crab pulsar true age and braking index are predicted to within $`20\%`$ or better. Consequently, it appears that a fallback disk with a current mass inflow rate in the range $`3\times 10^{16}10^{17}`$ g s-1 can account for the properties of the Crab pulsar with an accuracy of $`20\%`$ or so.

The model can be further tested with the third time derivative of the Crab pulsar frequency (Lyne et al. 1988; Blandford & Romani 1988). It is useful to define a second deceleration parameter:

$$p\frac{\stackrel{\mathrm{}}{\mathrm{\Omega }}\mathrm{\Omega }^2}{\dot{\mathrm{\Omega }}^3}=6+\left(3\frac{1}{\frac{1}{4}+\frac{\beta \mathrm{\Omega }^4}{4\gamma }}\right)^2\frac{1+\frac{\dot{\gamma }}{4\gamma }\frac{\mathrm{\Omega }}{\dot{\mathrm{\Omega }}}}{\frac{1}{4}+\frac{\beta \mathrm{\Omega }^4}{4\gamma }}+\frac{\frac{\ddot{\gamma }}{\gamma }\frac{\mathrm{\Omega }^2}{\dot{\mathrm{\Omega }}^2}\frac{\dot{\gamma }}{\gamma }\frac{\mathrm{\Omega }}{\dot{\mathrm{\Omega }}}}{1+\frac{\beta \mathrm{\Omega }^4}{\gamma }},$$
(17)
which has the correct limit $`p=15`$ for $`\gamma 0`$ (pure magnetic dipole radiation). The expression for steady mass inflow is found by forcing $`\ddot{\gamma }=\dot{\gamma }=0`$ in the above equation. Note again that $`\alpha (\alpha +1)\gamma /\ddot{\gamma }`$ is the square of the pulsar true age for a power law evolution of the disk mass inflow rate ($`\gamma (t)t^\alpha `$), while $`\mathrm{\Omega }/\dot{\mathrm{\Omega }}`$ is twice the pulsar timing age. The measured value of $`p`$ for the Crab is $`10.2`$ (Lyne et al. 1988; Blandford & Romani 1988; Melatos 1997). We find, for the case of steady mass inflow, that the equation above predicts $`p=12.6`$ for $`\dot{M}=3\times 10^{16}`$ g s-1 and $`p=9.46`$ for $`\dot{M}=10^{17}`$ g s-1. For a time dependent evolution with $`\alpha =1.2`$, we find $`p=13.6`$ and $`p=11.96`$ for the same mass inflow rates, respectively. The second deceleration parameter is therefore correctly predicted by the model (with an accuracy of $`30\%`$ or so) for fallback disk mass inflow rates in the range $`3\times 10^{16}10^{17}`$ g s-1.

### 3.2 Other Young Radio Pulsars

Braking indices have been measured for four other young radio pulsars. The measured values are $`n=2.2\pm 0.1`$, $`1.4\pm 0.2`$, $`2.91\pm 0.05`$ and $`2.837\pm 0.001`$ for PSR B0540-69, PSR B0833-45 (Vela), PSR J1119-6127 and PSR B1509-58, respectively (Lyne, Pritchard & Graham-Smith 1993; Deeter, Nagase & Boynton 1999; Lyne et al. 1996; Kaspi et al. 1994; Camilo et al. 2001).

It is possible to calculate, within the framework of our pulsar spindown model, the disk mass inflow rates required in each system to account for the measured breaking indices. In the absence of known true ages for these four pulsars, we calculate theoretical braking indices in two limits: (1) the steady mass inflow case (Eq. 6) and (2) the time-dependent case (Eq. 9) with the additional simplifying assumption that the pulsar timing age is a correct measure of the pulsar true age (the validity of this assumption, appropriate for the magnetic dipole model, was recently challenged by Gaensler & Frail 2000). The pulsar properties are taken from Taylor et al. (1993) and Camilo et al. (2001).

We find that the pulsar spindown model can account for the braking indices of PSR B0540-69, PSR J1119-6127 and PSR B1509-58 for disk mass inflow rates $`3\times 10^{16}`$, $`<\mathrm{\hspace{0.33em}10}^{14}`$ and $`<\mathrm{\hspace{0.33em}10}^{15}`$ g s-1, respectively. In each case, the propeller losses represent a reasonably small fraction of the total pulsar rotational losses. The very small value of $`n1.4`$ for the Vela pulsar cannot be accounted for by a model in which propeller losses are small. We note that the sub-Eddington mass inflow rates required for PSR B0540-69, PSR J1119-6127 and PSR B1509-58 (as well as the Crab) appear consistent with the values expected at late stages (before the fallback disk becomes neutral) in the evolutionary scenario described by Menou et al. (2001). One expects the magnetospheric interaction to be largely reduced or even suppressed when the gaseous disk becomes neutral (independently of the passive-disk outcome advocated in Menou et al. 2001).

### 3.3 Observational Constraints on Disk Sizes

The presence of disks around pulsars cannot be hidden: even a relatively compact disk will produce some optical emission, as radiation at these wavelengths is produced at relatively small radii.

Optical emission is a common feature of pulsars: since the first detection of the Crab in 1969, several more pulsars have been observed in the optical bandpass, and some of them have been seen to pulsate. In the cases where no pulsations are seen, the optical radiation is believed to be mainly of thermal origin, while pulsed optical emission is generally considered to be either a mixture of thermal and non-thermal (magnetospheric) radiation, or purely magnetospheric. Pacini (1971) and Pacini & Salvati (1987) proposed that the high energy emission comes from relativistic electrons radiating via synchrotron processes in the outer regions of the magnetosphere.

The Crab, Vela, and PSR B0540-69 pulsars belong to the above class, all possessing a pulsating optical counterpart. On the other hand, no pulsations have been detected for the candidate optical counterpart of the source PSR B1509-58 proposed333The interpretation of this optical counterpart is however a bit controversial (Mignani et al. 1998; Chakrabarty & Kaspi 1998; Gvaramadze 2001). by Caraveo et al. (1994).

In the following, we estimate, for each of the sources that we consider (except for PSR J1119-6127, for which no optical limits are available so far), the maximum size allowed for a disk to be compatible with the observed optical limits. For the three sources for which the optical emission is pulsed, we consider as an upper limit the lowest value of the emission (i.e. the minimum of the light curve), while for PSR B1509-58 we consider the detected value, assuming that it is indeed associated with the pulsar.

The inner radius of the disk is assumed to coincide with the light cylinder radius, $`R_{lc}`$. The disk model that we adopt is similar to that described by Perna, Hernquist & Narayan (2000). We include both the flux due to viscous dissipation and that due to reradiation. The irradiating X-ray luminosity is assumed to be isotropic, because, due to light bending in the vicinity of the star, beaming effects are smoothed out (Perna & Hernquist 2000).

Using the X-ray luminosities and periods from observations (see Becker & Trumper 1997 for a review), we find (assuming a disk inclination angle of $`60^o`$) that, for a mass inflow rate $`10^{17}`$ g s-1 a disk of radial extent $`6\times 10^9`$ cm is compatible with the optical limits for the Crab pulsar. A large disk ($`10^{13}`$ cm) is allowed for both PSR B0540-69 and PSR B1509-58 at the mass inflow rate of $`3\times 10^{16}`$ and $`10^{15}`$ g s-1 respectively, while for Vela only a very small disk ($`10^9`$ cm) at an high inclination ($`89^o`$) would be compatible with the observational constraints444Note however that, as the inner radii are $`10^810^9`$ cm for all the pulsars considered, we expect most of the intrinsic light from the disk to be produced in the UV, which is difficult to detect.. The allowed size becomes larger as $`\dot{M}`$ decreases. Again, we note that, with the exception of the Vela pulsar, these disk sizes appear consistent with the values expected at late stages (before the fallback disk becomes neutral) in the evolutionary scenario of Menou et al. (2001).

Support for this work was provided by NASA through Chandra Fellowship grant PF9-10006 awarded by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-39073.


text:
# An alternative to quintessence

## Abstract

We consider a FRW cosmological model with an exotic fluid known as Chaplygin gas. We show that the resulting evolution of the universe is not in disagreement with the current observation of cosmic acceleration. The model predict an increasing value for the effective cosmological constant.

The discovery that the expansion of the universe is accelerating  has promoted the search for new types of matter that can behave like a cosmological constant  by combining positive energy density and negative pressure. This type of matter is often called ”quintessence”.

Since in a variety of inflationary models scalar fields have been used in describing the transition from the quasi-exponential expansion of the early universe to a power law expansion, it is natural to try to understand the present acceleration of the universe, which has an exponential behaviour too, by constructing

models where the matter responsible for such behaviour is also represented by a scalar field . However, now we deal with the opposite task, i.e. we would like to describe the transition from a universe filled with dust-like matter to an exponentially expanding universe, and scalar fields are not the only possibility but there are (of course) alternatives. In particular, one can try to do it by using some perfect fluid but obeying “exotic” equations of state.

In this short letter we consider an example of this type: the so-called Chaplygin gas. Under this name we mean a perfect fluid having the following equation of state:

$$p=\frac{A}{\rho },$$
(1)
where $`p`$ and $`\rho `$ are respectively pressure and energy density in a comoving reference frame, with $`\rho >0`$; $`A`$ is a positive constant. This equation of state has raised recently a certain interest  because of its many interesting and, in some sense, intriguingly unique features. Indeed Eq. (1) has an amusing connection with string theory and it can be obtained from the Nambu-Goto action for $`d`$-branes moving in a $`(d+2)`$-dimensional spacetime in the light-cone parametrization . Also, the Chaplygin gas is the only fluid which, up to now, admits a supersymmetric generalization . We ourselves came across this fluid  when studying the stabilization of branes  in black hole bulks . We found that to obtain stabilization it is necessary to add matter on the branes which again obeys the equation of state (1).

For these reasons we have undertaken the simple exercise of stydying a FRW cosmology of a universe filled with a Chaplygin gas. The metric of a homogeneous and isotropic universe is usually written as follows

$$ds^2=dt^2a^2(t)dl^2,$$
(2)
where $`dl^2`$ is the metric of a 3-manifold of constant curvature ($`K=0,\pm 1`$), and the expansion factor $`a(t)`$ evolves according to the Friedmann equation

$$\frac{\dot{a}^2}{a^2}=\rho \frac{K}{a^2}.$$
(3)
Energy conservation

$$d(\rho a^3)=pd(a^3)$$
(4)
together with the equation of state (1) give the following relation:

$$\rho =\sqrt{A+\frac{B}{a^6}},$$
(5)
where $`B`$ is an integration constant.

By choosing a positive value for $`B`$ we see that for small $`a`$ (i.e. $`a^6B/A`$) the expression (5) is approximated by

$$\rho \frac{\sqrt{B}}{a^3}$$
(6)
that corresponds to a universe dominated by dust-like matter. For large values of the cosmological radius $`a`$ it follows that

$$\rho \sqrt{A},p\sqrt{A},$$
(7)
which, in turn, corresponds to an empty universe with a cosmological constant $`\sqrt{A}`$ (i.e a de Sitter universe). In the flat case it is possible also to find exact solutions as follows:

$`t={\displaystyle \frac{1}{6^4\sqrt{A}}}\left(\mathrm{ln}{\displaystyle \frac{{}_{}{}^{4}\sqrt{A+\frac{B}{a^6}}+^4\sqrt{A}}{{}_{}{}^{4}\sqrt{A+\frac{B}{a^6}}^4\sqrt{A}}}2\mathrm{arctan}^4\sqrt{1+{\displaystyle \frac{B}{Aa^6}}}\right)`$ (8)
Note that $`\sqrt{A}`$ solves the equation

$$\rho +p=\rho \frac{A}{\rho }=0.$$
(9)
The circumstance that this equation has a nonzero solution lies at the heart of the possibility of interpreting the model as a “quintessential” model. If this model were realistic we could estimate the constant $`A`$ by comparing our expressions for pressure and energy with observational data. An indirect and naive way to do it is to consider the nowadays accepted values for the contributions of matter and cosmological constant to the energy density of the universe. To use these data we decompose pressure and energy density as follows:

$`p=p_\mathrm{\Lambda }+p_M=\mathrm{\Lambda },`$ (10)
$`\rho =\rho _\mathrm{\Lambda }+\rho _M=\mathrm{\Lambda }+\rho _M.`$ (11)
An application of Eq. (1) gives

$$A=\mathrm{\Lambda }(\mathrm{\Lambda }+\rho _M).$$
(12)
If the cosmological constant contributes seventy percent of the energy we get $`\sqrt{A}1.2\mathrm{\Lambda }`$. We now observe that, in the context of a Chaplygin cosmology, once an expanding universe starts accelerating it cannot decelerate any more. Indeed eqs. (3) and (4) imply that

$$\frac{\ddot{a}}{a}=\frac{1}{2}(\rho +3p).$$
(13)
Condition $`\ddot{a}>0`$ is equivalent to

$$a^6>\frac{B}{2A},$$
(14)
which is obviously preserved by time evolution in an expanding universe. It thus follows that the observed value $`\mathrm{\Lambda }`$ of the (effective) cosmological constant will increase up to $`1.2\mathrm{\Lambda }`$.

Considering now the subleading terms in Eq. (5) at large values of $`a`$ (i.e. $`a^6B/A`$), one obtains the following expressions for the energy and pressure:

$`\rho \sqrt{A}+\sqrt{{\displaystyle \frac{B}{4A}}}a^6,`$ (15)
$`p\sqrt{A}+\sqrt{{\displaystyle \frac{B}{4A}}}a^6.`$ (16)
Eqs. (15) and (16) describe the mixture of a cosmological constant $`\sqrt{A}`$ with a type of matter known as “stiff” matter, and described by the following equation of state:

$$p=\rho .$$
(17)
Note that a massless scalar field is a particular instance of stiff matter. Therefore, in a generic situation, a Chaplygin cosmology can be looked at as interpolating between different phases of the universe: from a dust dominated universe to a de Sitter universe passing through an intermediate phase which is the mixture just mentioned above. The interesting point however is that such an evolution is accounted by using only one fluid.

In recent series of paper  a similar type of evolution has been described, where the universe passes from a dust dominated epoch to a de Sitter phase through an intermediate phase described as mixture of cosmological constant and radiation. In  the mechanism responsible for this behaviour is different from ours and is based on the quantum corrections to the effective action of a massive scalar field. However these corrections lead to a “standard” equation of state in the form of a mixture. We can reproduce this type of evolution by a slight modification to our “exotic” equation of state (1), namely

$$p=\frac{C}{{}_{}{}^{3}\sqrt{\rho }}$$
(18)
The obvious generalization $`p=C\rho ^\alpha `$ (with $`\alpha 1`$) gives a similar evolution from dust to cosmological constant with an intermediate epoch that can be seen as a mixture of a cosmological constant with a fluid obeying the state equation $`p=\alpha \rho `$.

For open or flat Chaplygin cosmologies ($`K=1,0`$), the universe always evolves from a decelerating to an accelerating epoch. For the closed Chaplygin cosmological models ($`K=1`$), the Friedmann equations (3) and (13) say that it is possible to have a static Einstein universe solution $`a_0=(3A)^{\frac{1}{4}}`$ provided the following condition holds:

$$B=\frac{2}{3\sqrt{3A}}$$
(19)
When $`B>\frac{2}{3\sqrt{3A}}`$ the cosmological radius $`a(t)`$ can take any value while if

$$B<\frac{2}{3\sqrt{3A}}$$

there are two possibilities: either

$$a<a_1=\frac{1}{\sqrt{3A}}\left(\sqrt{3}\mathrm{sin}\frac{\phi }{3}\mathrm{cos}\frac{\phi }{3}\right)$$
(20)
or

$$a>a_2=\frac{2}{\sqrt{3A}}\mathrm{cos}\frac{\phi }{3},$$
(21)
where $`\phi =\pi \mathrm{arccos}3\sqrt{3A}B/2`$. The region $`a_1<a<a_2`$ is not accessible.

Results connected to those presented in this letter have been obtained by Barrow  who has considered cosmologies with fluids admitting a bulk viscosity proportional to a power of the density.

In the flat $`K=0`$ case the FRW equations for Chaplygin fits in Barrow’s scheme as a special case  and indeed a transition from power law to exponential expansion is noticed already in . However since the state equation for the corresponding fluid is different from our Eq.(1) this coincidence of the solutions is destroyed by any small perturbation, for instance by a small spatial curvature or by adding another matter source. We mention also that the role of the Chaplygin-like behaviour in cosmology was also noticed in Ref. .

Following  we now try to find a homogeneous scalar field $`\varphi (t)`$ and a potential $`V(\varphi )`$ to describe the Chaplygin cosmology. Let us consider therefore the Lagrangian

$$L(\varphi )=\frac{1}{2}\dot{\varphi }^2V(\varphi )$$
(22)
and set the energy density of the field equal to that of the Chaplygin gas:

$$\rho _\varphi =\frac{1}{2}\dot{\varphi }^2+V(\varphi )=\sqrt{A+\frac{B}{a^6}}.$$
(23)
The corresponding ”pressure” coincides with the Lagrangian density:

$$p_\varphi =\frac{1}{2}\dot{\varphi }^2V(\varphi )=\frac{A}{\sqrt{A+\frac{B}{a^6}}}.$$
(24)
It immediately follows that

$$\dot{\varphi }^2=\frac{B}{a^6\sqrt{A+\frac{B}{a^6}}}$$
(25)
and

$$V(\varphi )=\frac{2a^6\left(A+\frac{B}{a^6}\right)B}{2a^6\sqrt{A+\frac{B}{a^6}}}.$$
(26)
Now let us restrict ourselves to the flat case $`K=0`$. In this case Eq. (25) also implies that

$$\varphi ^{}=\frac{\sqrt{B}}{a(Aa^6+B)^{1/2}},$$
(27)
where prime means differentiation w.r.t. $`a`$. This equation can be integrated and it follows that

$$a^6=\frac{4B\mathrm{exp}(6\varphi )}{A(1\mathrm{exp}(6\varphi ))^2}.$$
(28)
Finally, by substituting the latter expression for the cosmological radius in Eq. (26) one obtains the following potential, which has a surprisingly simple form:

$$V(\varphi )=\frac{1}{2}\sqrt{A}\left(\mathrm{cosh}3\varphi +\frac{1}{\mathrm{cosh}3\varphi }\right).$$
(29)
Note that the potential does not depend on the integration constant $`B`$ and therefore it reflects only the state equation (1) as it should.

In conclusion Chaplygin cosmology provides an interesting possibility to account for current observations about the expansion of the universe. It predicts also that the cosmological constant will increase (or that it was less in the past) and this could in principle be observed. Of course to take this model seriously one should have a good fundamental reason to believe in Eq. (1).

Acknowledgments: We are grateful to Prof. J.D. Barrow who has kindly attracted our attention to his papers . The work of A.K. was partially supported by RFBR via grants No 99-02-18409 and No 99-02-16224 and by the Cariplo Science Foundation.