text:
# Bose-Einstein condensate collapse: a comparison between theory and experiment

## Abstract

We solve the Gross-Pitaevskii equation numerically for the collapse induced by a switch from positive to negative scattering lengths. We compare our results with experiments performed at JILA with Bose-Einstein condensates of 85Rb, in which the scattering length was controlled using a Feshbach resonance. Building on previous theoretical work we identify quantitative differences between the predictions of mean-field theory and the results of the experiments. Besides the previously reported difference between the predicted and observed critical atom number for collapse, we also find that the predicted collapse times systematically exceed those observed experimentally. Quantum field effects, such as fragmentation, that might account for these discrepancies are discussed.

Introduction.– Most experiments on dilute gas Bose-Einstein condensates (BECs) are performed with atoms that have a repulsive two-body interaction. Exceptions are the experiments on 7Li Bradley ; Gerton  and, more recently, on 85Rb Donley ; Robertscollapse . For 85Rb a Feshbach resonance allows the two-body interaction strength to be tuned over a wide range of attractive and repulsive values. In particular, the scattering length has been rapidly switched from positive (repulsive interaction) to negative (attractive interaction) values, leading to the collapse and subsequent explosion of the condensate. Recently, the large positive scattering lengths attainable in this system have been used to produce atom-molecule condensates Donleymolecule .

In the following we report on our modelling of the 85Rb collapse experiments Donley , using the Gross-Pitaevskii (GP) equation for the expectation value of the field operator Dalfovo ; RMP ; MeystreAObook ; PethickSmith . Saito and Ueda Saitoexplosion  and Adhikari Adhikariexplosion  have also modelled these experiments by numerical solution of the cylindrically symmetric GP equation. Saito and Ueda conclude that this describes the collapsing and exploding dynamics at least qualitatively Saitoexplosion . Following their suggestion, we report a more quantitative comparison between the theoretical and experimental results, and find significant differences.

The series of experiments on the collapse and explosion of 85Rb BECs challenges theoretical models in a number of ways Claussen . A body of theoretical work based on the GP equation predicts the critical number of atoms $`N_{\mathrm{cr}}`$ for collapse to be significantly larger than is observed. The expression for the critical number is

$$N_{\mathrm{cr}}=k\frac{a_{\mathrm{ho}}}{|a|},$$
(1)
where $`a_{\mathrm{ho}}=\sqrt{\mathrm{}/(m\overline{\omega })}`$ is the harmonic oscillator scale length, with $`\overline{\omega }`$ the geometric mean of the trap frequencies in the three Cartesian directions, and $`a`$ the scattering length. Experimentally, $`k=0.46\pm 0.06`$ Robertscollapse , whereas $`k=0.57`$ according to various approximate solutions of the GP equation EleftheriouHuang ; Ruprecht .

We have confirmed this GP prediction for the specific cylindrically symmetric experimental case Robertscollapse  with cylindrically symmetric numerical solutions. We verified these with full three dimensional numerical solutions, and also confirmed that slight departures from cylindrical symmetry had no effect on the critical number 3dnumerics . Consequently there is a disagreement at the two standard deviations level, which should be regarded as significant.

We also report a new quantitative discrepancy between the predictions of the GP model and experiment. Under certain conditions, the GP predicted time to the initiation of collapse, $`t_{\mathrm{collapse}}`$, is systematically longer than that observed in the experiments Donley .

The GP model.– In the conclusion we will discuss the possiblility that these discrepancies result from quantum field effects beyond the GP approximation. We therefore now derive the GP equation from the quantum field theory.

The second-quantised Hamiltonian for a dilute gas, in terms of the field operator $`\widehat{\mathrm{\Psi }}(𝐫,t)`$, is

$`H`$ $`=`$ $`{\displaystyle 𝑑𝐫\widehat{\mathrm{\Psi }}^{}H_0\widehat{\mathrm{\Psi }}}`$ (2)
$`+{\displaystyle \frac{1}{2}}{\displaystyle 𝑑𝐫𝑑𝐫^{}\widehat{\mathrm{\Psi }}^{}\widehat{\mathrm{\Psi }}^{}V(𝐫𝐫^{})\widehat{\mathrm{\Psi }}^{}\widehat{\mathrm{\Psi }}},`$
where $`\widehat{\mathrm{\Psi }}^{}=\widehat{\mathrm{\Psi }}(𝐫^{},t)`$ and $`H_0`$ is the single particle Hamiltonian for the kinetic energy and trapping potential

$$H_0=\frac{\mathrm{}^2}{2m}^2+\frac{1}{2}m(\omega _x^2x^2+\omega _y^2y^2+\omega _z^2z^2),$$
(3)
where $`m`$ is the atomic mass ($`1.41\times 10^{25}`$ kg for 85Rb), and $`\omega _i`$ is the trap frequency along Cartesian axis $`i`$. In the limit of particles separated by distances much greater than the scattering length $`a`$ we approximate the two-body potential by a delta function interaction Dalfovo ; RMP ; MeystreAObook ; PethickSmith

$$V(𝐫𝐫^{})=g\delta (𝐫𝐫^{}),g=\frac{4\pi \mathrm{}^2a}{m}.$$
(4)
The Heisenberg dynamical equation for the field operator is then

$$i\mathrm{}\frac{}{t}\widehat{\mathrm{\Psi }}=H_0\widehat{\mathrm{\Psi }}+g\widehat{\mathrm{\Psi }}^{}\widehat{\mathrm{\Psi }}\widehat{\mathrm{\Psi }},$$
(5)
Taking the symmetry-breaking approach we assume that the field expectation value is not zero and define it as the GP wavefunction $`\widehat{\mathrm{\Psi }}(𝐫,t)=\mathrm{\Phi }(𝐫,t)`$, normalised to the number of particles $`N`$

$$N=|\mathrm{\Phi }(𝐫,t)|^2𝑑𝐫.$$
(6)
Then taking the expectation value of the Heisenberg equation (5) gives

$$i\mathrm{}\frac{}{t}\mathrm{\Phi }=H_0\mathrm{\Phi }+g\widehat{\mathrm{\Psi }}^{}\widehat{\mathrm{\Psi }}\widehat{\mathrm{\Psi }},$$
(7)
If we assume that the expectation value factorises, as it would, for example, if the system were in an eigenstate of the field operator,

$$\widehat{\mathrm{\Psi }}^{}\widehat{\mathrm{\Psi }}\widehat{\mathrm{\Psi }}=\mathrm{\Phi }^{}\mathrm{\Phi }\mathrm{\Phi },$$
(8)
then we obtain the GP equation

$$i\mathrm{}\frac{}{t}\mathrm{\Phi }=(H_0+g|\mathrm{\Phi }|^2)\mathrm{\Phi }.$$
(9)
In order to model atom loss due to three-body recombination we add a phenomenological term proportional to the density $`|\mathrm{\Phi }|^2`$ squared with rate coefficient $`K_3/2`$ Robertsparams ; PethickSmith

$$i\mathrm{}\frac{}{t}\mathrm{\Phi }=(H_0+g|\mathrm{\Phi }|^2i\frac{\mathrm{}}{2}K_3|\mathrm{\Phi }|^4)\mathrm{\Phi }.$$
(10)
We assume one-body and two-body loss are negligible, as was true for the relevant experiments. The number of atoms then decays as

$$\frac{dN}{dt}=K_3|\mathrm{\Phi }(𝐫,t)|^6𝑑𝐫.$$
(11)
GP Results.– As an example of the ability of the GP equation to correctly model the 85Rb Donley  experiments we present Fig. 1. It is the result of a numerical solution of the (two dimensional) cylindrically symmetric GP equation for $`\stackrel{~}{\mathrm{\Phi }}(r,z)`$

$`i\mathrm{}{\displaystyle \frac{}{t}}\stackrel{~}{\mathrm{\Phi }}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}(_r^2+r^1_r+_z^2)\stackrel{~}{\mathrm{\Phi }}`$ (12)
$`+{\displaystyle \frac{1}{2}}m(\omega _r^2r^2+\omega _z^2z^2)\stackrel{~}{\mathrm{\Phi }}+g|\stackrel{~}{\mathrm{\Phi }}|^2\stackrel{~}{\mathrm{\Phi }}`$
$`i{\displaystyle \frac{\mathrm{}}{2}}K_3|\stackrel{~}{\mathrm{\Phi }}|^4\stackrel{~}{\mathrm{\Phi }}.`$
Parameters are the same those of Fig. 1b of Donley et al. Donley . Specifically, the ground state of the GP equation for $`a=+7a_0`$ was switched in 1 ms to $`a=30a_0`$, where $`a_0=0.0529`$ nm is the Bohr radius. For the three-body recombination rate coefficient $`K_3=190\times 10^{28}`$ cm6s-1 the agreement with the experimental results is good. However it should be noted that the experimental points are the “remnant” atom number, while the numerical points are the total atom number, which overestimates the remnant atom number. A smaller value of $`K_3`$ agrees better with the earlier points, while overestimating the final atom number. The precise value of $`K_3`$ has little effect on the conclusions of this paper, which concern the initiation of collapse.

These results agree with those reported by Saito and Ueda Saitoexplosion  and Adhikari Adhikariexplosion . However the former authors used a much smaller value of the three-body recombination rate coefficient $`K_3=2\times 10^{28}`$ cm6s-1. This produces the collapses and revivals in condensate size that were observed in their simulations. These only become important for $`K_3`$ less than about $`10^{26}`$ cm6s-1. Adhikari Adhikariexplosion  used the much larger value $`K_3=13\times 10^{25}`$ cm6s-1. Since three-body recombination is responsible for the atom loss, it is remarkable that such a wide range of coefficients reproduces the experimental results.

The three-body recombination rate coefficient $`K_3`$ is expected to vary strongly near the Feshbach resonance Esry . Experimental determination of $`K_3`$ is difficult due to the low densities of 85Rb condensates. Upper bounds have been estimated to be $`5\times 10^{25}`$ cm6s-1, dropping to $`10^{26}`$ cm6s-1 nearer the Feshbach resonance Robertsparams .

The cylindrically symmetric numerical simulations were performed on a $`512\times 512`$ grid, 35.64 $`\mu `$m long in the axial ($`z`$) direction and with the radial coordinate extending to 11.88 $`\mu `$m. The corresponding spatial grid spacings were therefore 0.07 $`\mu `$m and 0.023 $`\mu `$m. The time steps were $`2.34`$ ns. All simulations were performed on a multiprocessor machine APAC , using up to 32 processors, and the RK4IP algorithm developed by the BEC theory group of R. Ballagh at the University of Otago RK4IP . This is a pseudo-spectral method with a Runge-Kutta time step. The cylindrically symmetric and full three-dimensional codes are independent and were cross checked. Grid spacings and time steps were varied to ensure convergence. Overall the results were found to be quite robust. As another test, we solved the GP equation for a half a radial period after the quenching of the collapse. As was observed experimentally, the condensate refocussed onto the axis, due to the oscillation in the harmonic radial potential. All this, together with the agreement of our results with those of Saito and Ueda Saitoexplosion  and Adhikari Adhikariexplosion , gives us confidence in their accuracy. Following Adhikari Adhikariexplosion , the initial condition for Fig. 1 was generated by adiabatically expanding the harmonic oscillator initial state $`a=0`$ to $`a=+7a_0`$ over 444 ms.

Figure 2 presents our calculations of the collapse times $`t_{\mathrm{collapse}}`$ for the conditions of Fig. 2 of Donley et al. Donley . The collapse times were determined by visually fitting plots of atom number versus time to the functional form

$$N=(N_0N_\mathrm{f})\mathrm{exp}[(tt_{\mathrm{collapse}})/\tau _{\mathrm{decay}}]+N_\mathrm{f},$$
(13)
where $`N_\mathrm{f}`$ is the long time atom number. An example is given in the inset to Fig. 2. We have also plotted the experimental results reported in Fig. 2 of Donley et al. Donley , and find a small, but significant, systematic disagreement with the GP results. Although the reported errors in the experimental collapse times are large, the GP values for $`t_{\mathrm{collapse}}`$ are consistently longer than the experimental ones. This is surprising as the GP model is expected to be valid for the low densities preceding the collapse. If it were to fail, it would be expected to do so at the high densities generated subsequently. Nevertheless, the disagreement is not unprecedented since, as we discussed earlier, the GP model also overestimates the critical number for collapse.

The estimates of $`t_{\mathrm{collapse}}`$ by Saito and Ueda Saitoexplosion  (their Fig. 3) are between five percent (low $`a`$) and ten percent (high $`a`$) smaller than ours. This is consistent with the smaller three-body recombination rate coefficient $`K_3`$ they used. However their results are still significantly longer than the experimentally measured times.

We have confirmed these cylindrically symmetric simulations by performing full three-dimensional simulations. In particular we broke the cylindrical symmetry by using trap frequencies of $`17.24\times 17.47\times 6.80`$ Hz Robertscollapse .

We were unable to substantially improve the agreement by either changing the initial condition to reflect the experimental uncertainty of $`a=\pm 2a_0`$, or by varying the three-body recombination rate coefficient. This suggests that some of the physics determining the collapse time is not captured by our GP model.

Discussion.– Both the collapse time and critical number discrepancies could be resolved by using a scattering length in the GP model larger in magnitude than the experimental value. This would reduce the collapse time and decrease the critical number, as required. The required increases in the scattering length magnitudes vary, ranging from a factor of $`0.57/0.46=1.2`$ for the critical number, up to a factor of about two for the collapse times for large $`a_{\mathrm{collapse}}`$. However, the scattering length is experimentally well calibrated Robertscalibrate , so any such change would reflect a deficiency of our GP model.

One possible origin of the discrepancy is the effect of thermal non-condensed atoms. Because of the quantum statistics of collisions betwen bosons, the scattering length between a condensed atom and an atom in another mode is twice that between two condensed atoms. Hence one might expect the presence of thermal uncondensed atoms to shorten the collapse time compared to the GP prediction, as observed. Furthermore, this might be approximately corrected for by using an increased magnitude effective scattering length in the GP model. However the uncondensed fraction is much less than 10% of the total number of atoms Claussen , so it seems unlikely that this effect is large enough to account for the discrepancy. Furthermore, Roberts et al. Robertscollapse  reported that the critical number for collapse $`N_{\mathrm{cr}}`$ was insensitive to varying the temperature. Therefore we do not expect finite temperature extensions of the GP theory to explain the discrepancy finite temperature .

Another possibile origin is quantised atom field effects. These might arise due a breakdown of the factorisation assumption Eq.(8). There have been several suggestions for how the quantised field might influence the collapse UedaLegget ; DuineStoof ; Yurovsky . Furthermore, Nozieres Nozieres  has emphasised that only for positive scattering lengths does an energy barrier protect BECs from fragmentation into many populated states. For negative scattering lengths, mean-field energy is released when atoms scatter from the condensate into other modes. Fragmentation could increase the effective scattering length by up to a factor of two.

In order to investigate the behaviour of a fully quantised atom field, we have used the gauge-P function approach recently developed by Deuar and Drummond DeuarDrummond . This method overcomes some of the problems that plague stochastic simulations based on the positive P-function quasi-probability distribution Drummond ; Steel . We were computationally limited to simulations in one spatial dimension and found agreement with the GP collapse times at the one percent level.

Although this preliminary work does not provide evidence for quantum field effects, it is important to extend the fully quantised field modelling to three spatial dimensions, and hence to use actual experimental parameters. As shown in Fig. 1, there are parameters for which the GP theory does agree with experiments. One approach is the recently developed perturbation theory which extends the GP model to include normal and anomalous densities of the quantum field HollandParkWalser . This method has recently been successfully applied to the formation of atom-molecule condensates in 85Rb KokkelmansHolland .

###### Acknowledgements.

This research was supported by the Australian Research Council and by an award under the Merit Allocation Scheme on the National Facility of the Australian Partnership for Advanced Computing APAC .


format: text

text:
# 1 Introduction

## 1 Introduction

Radio and optical polarimetry of the large scale jets in radio galaxies usually suggest strong velocity shearing effects at the jet edges, resulting in a magnetic field being parallel to the jet axis (see, e.g., Perlman et al. 1999 for M 87, Attridge et al. 1999 for 1055+018). Acceleration of radiating particles within a boundary region is required to explain the observed jet radial profile (e.g. in M87, Owen et al. 1989), or details of the spectral maps, like flattening of the radio spectrum near the jet surface in Mkn 501 (Edwards et al. 2000). One should mention, that the HST studies of the optical counterparts of large scale radio jets emphasize the necessity of the continuous electron reacceleration in a jet body (and not only in shock regions), as the lifetime of the synchrotron electrons is often much shorter than the light-travel time along the optical structure (cf. Jester et al. 2001 for the case of 3C 273). On the other hand, the broad band spectral properties of BL Lacs and FR I radio galaxies indicate a jet stratification. By comparing the multiwavelength observations of these two classes of AGNs in a framework of a unification scheme, Chiaberge et al. (2000) found evidences for a significant boundary layer emission, dominating the radiative jet output in FR Is. The inferred velocities of such boundary regions are significantly lower as compared to velocities of central spines, but are still relativistic, in order to explain anisotropic emission observed in the FR Is’ cores.

3D hydrodynamical simulations of relativistic jets reveal a highly turbulent cocoon and a shear layer surrounding the jet spine (Aloy et al. 1999). The reason for the presence of a turbulent medium effects in these simulations is numerical viscosity, which is not a real physical process. However, theoretical considerations support the model of a turbulent medium at the jet edges (cf. Bicknell & Melrose 1982). Such regions are therefore the promising places for particle acceleration, as discussed previously by Ostrowski (1998, 2000). The jet boundary layer acceleration, when applied to protons, was considered to provide ultra high energy cosmic rays, influence the jet dynamics, and significantly increase a pressure in the radio lobes of FR II sources (Ostrowski & Sikora 2001). Electrons accelerated at the jet boundary were considered to provide important contribution for the jet radiative output, including the large scale X-ray emission and its possible time (spatial) variation (Stawarz & Ostrowski 2002, 2001).

## 2 Radiating boundary layer

For illustration let us consider a relativistic large scale jet consisted of a spine surrounded by the boundary layer with a velocity shear. In order to specify the velocity structure, we assume a uniform flow Lorentz factor $`\mathrm{\Gamma }_j`$ inside the jet spine, and a linear radial profile $`1<\mathrm{\Gamma }<\mathrm{\Gamma }_j`$ within the boundary transition region. Below we put the thickness of the boundary layer to be of order of the jet radius, $`D1\mathrm{kpc}`$. Furthermore, we assume, that due to shearing effects the magnetic field is parallel on average to the jet axis inside the highly turbulent boundary layer, and that its intensity is the same as the one estimated for the spine, $`B=10^5\mathrm{G}`$.

The electrons injected into the jet boundary region can be accelerated due to cosmic ray viscosity connected with the flow velocity radial gradient within the shear layer and due to stochastic Fermi acceleration in the turbulent medium. For the assumed large scale jet parameters (with $`Dr_g`$, where $`r_g`$ is the electron gyroradius) the later mechanism acts more efficiently and, hence, the electrons gain energy mostly due to scattering on the long wavelength magnetic field irregularities moving with the velocity comparable to the Alfvén velocity $`V_A`$. The acceleration time scale for this process is roughly $`T_{acc}\zeta r_gcV_A^2`$, where $`\zeta `$ is the ratio of the energy density of the regular magnetic field to the turbulent one. For the jet dynamicaly dominated by nonrelativistic protons, $`V_A`$ is a few orders of magnitudes smaller than the flow velocity $`Uc`$. Below we assume $`V_A0.01c`$, what corresponds to the cold proton number density $`n_p10^4\mathrm{cm}^3`$ (cf. Sikora & Madejski 2000).

Radiative losses in the large scale jets are due to the synchrotron emission and comptonization of cosmic microwave background (CMB) photons and synchrotron core emission (AGN) illuminating the large scale jet from behind (cf. Celotti et al. 2001). In the Thomson regime one can write $`T_{loss}3m_ec/4\sigma _T\gamma u`$, where $`u=u_B+u_{cmb}+u_{agn}`$ is the energy density of the magnetic field and the seed photons in the source frame. The latter two represent anisotropic radiation fields which depend on the flow Lorentz factor of the emitting region, and hence – in our case of a shear layer – on the distance from the jet axis. For simplicity, we limit our calculation to the upper limits of both $`u_{cmb}`$ and $`u_{agn}`$, and neglect cosmological redshift corrections. Then the energy density of the CMB background is $`u_{cmb}=aT_{cmb}^4\mathrm{\Gamma }^2<aT_{cmb}^4\mathrm{\Gamma }_j^2`$, and the energy density of the core emission with intrinsic luminosity $`L_{agn}^{}`$ is $`u_{agn}=L_{agn}^{}\mathrm{\Gamma }_0^2/4\pi cz^2\mathrm{\Gamma }^2<L_{agn}^{}\mathrm{\Gamma }_0^2/4\pi cz^2`$, where $`\mathrm{\Gamma }_0`$ is the bulk Lorentz factor of the (sub-) parsec scale jet, and $`z`$ is the distance from the galactic nucleus. The electrons’ escape from the shear layer with the assumed magnetic field configuration proceeds through the cross-field diffusion process. The appropriate time scale for this process can be estimated as $`T_{esc}D^2\kappa _{}^1`$, where $`\kappa _{}\eta r_gc/3`$ is an effective cross-field diffusion coefficient, with the numerical scaling factor $`\eta 1`$. $`T_{esc}`$ is much longer than the time scale for radiative losses, $`T_{loss}`$. As a result, the electrons pile-up at highest energies below the maximum energy $`\gamma _{eq}mc^2`$, where $`T_{acc}T_{loss}`$. Below we will model the formed hard component in the electron spectrum with a monoenergetic bump at the highest energy ($`\delta (\gamma \gamma _{eq})`$, cf. Ostrowski 2000). For the parameters $`L_{agn}^{}10^{43}\mathrm{erg}/\mathrm{s}`$, $`\zeta =1`$ and $`\mathrm{\Gamma }_0\mathrm{\Gamma }_j10`$, the maximum energy the electrons can reach due to turbulent acceleration in a Bohm limit – at distances $`z10\mathrm{kpc}`$ where $`u_{agn}<u_{cmb}`$ – corresponds to the electron Lorentz factor $`\gamma _{eq}10^8`$.

In the case of efficient particle injection, the acceleration process acting continuously within the whole jet boundary layer leads to the distribution approximated here as a flat power-law spectrum $`n_e(\gamma )\gamma ^\sigma `$ finished with the pile-up bump at the maximum energies near $`\gamma _{eq}`$. The detailed evolution of the highest energy electrons needs careful analysis of the momentum diffusion of the accelerating particles within turbulent layer with velocity and density gradients (work in preparation). Below, in order to discuss main consequences of the boundary layer acceleration for the radiative jet output, we consider a case of forming the stationary two-component electron energy distribution with $`\sigma 2`$ and the high energy bump modeled as a monoenergetic peak (c.f. discussion in Stawarz & Ostrowski 2002). The observed multiwavelength radiation of such electrons, including synchrotron radiation, synchrotron self-Compton emission and external Compton scattering of CMB photons is plotted on figure 1 for different jet inclinations. In order to find the normalization constants for the particle distribution, we assumed the energy equipartition between the magnetic field and each one of the electron spectral components. The presented spectra are integrated over the assumed bulk flow Lorentz factor profile within the boundary layer, with the appropriate beaming pattern for synchrotron/self-Compton and external Compton radiation (Dermer 1995).

## 3 Discussion - X-ray emission of the large scale jets

After Chandra discovery that the large scale jets in quasars and radio galaxies are strong X-ray emitters, the question of the origin of such radiation became an important issue. For the jets observed at small angles, the inverse Compton scattering of the CMB photons by the low-energy tail of the non-thermal electron distribution flowing with large bulk Lorentz factor is the most likely explanation (Tavecchio et al. 2000). However, strong beaming effects exclude this process in case of the jets in radio galaxies. Celotti et al. (2001) proposed instead, that the X-ray radiation detected from such objects can be explained as the SSC emission of the slower (and therefore less beamed) boundary layer electrons, or as the inverse Compton scattering of the core emission. This, however, requires a departure from the energy equipartition between the magnetic field and radiating particles, at least if the X-ray flux is of the same order of magnitude as the synchrotron one. As illustrated on figure 1, the boundary layer acceleration at the tens-of-kpc scale jets can generate high energy synchrotron radiation peaking at the keV energy band. Luminosity of such radiation is high enough to account for the Chandra observations without departures from the equipartition condition (Stawarz & Ostrowski 2002). Continuous acceleration acting within the whole boundary layer volume compensates radiative losses of high energy electrons, and therefore the diffusive character of their X-ray emission is a natural consequence of the presented model (but of course any compressive perturbation – a shock – in the jet can disturb this smooth intensity variation). One should note that recent HST observations of the jet in 3C 273 (Jester et al. 2002) with a smooth spectral transition from optical, through UV, up to X-ray frequencies are consistent with our model.

The low energy (radio-to-optical) component of the boundary layer synchrotron emission can affect the jet-counterjet brightness asymmetry measurements, ‘hiding’ the highly relativistic spine. Thus the bulk Lorentz factor of the jet spine at large scales can be comparable to the one at VLBI scales. This has several consequences for the jet energetics (cf. Ghisellini & Celotti 2001).

## Acknowledgements

We are grateful to Marek Sikora for his help and discussions. The present work was supported by Komitet Badań Naukowych through the grant BP 258/P03/99/17.


text:
# Nonpartonic effects in pion electroproduction in the HERMES kinematical region

## 1 Introduction

Recently the interest in semi-inclusive processes has increased significantly. They appear to be a very useful tool for separating the flavour- and spin-dependent quark distributions in the nucleon. The interpretation of such experiments is based on the quark-parton model, often in its simplest form. Therefore, from a practical point of view, some of very important topics are: the presence of mechanisms beyond the QPM in the studied reactions, contributions of such mechanisms to measured cross sections and their influence on the analysed quantities.

The HERMES Collaboration has precisely measured multiplicities of charged and neutral pions produced in unpolarized DIS off the proton . Previously, charged pion production off proton and deutron targets was used for the extraction of the $`\overline{d}\overline{u}`$ asymmetry . Since the HERMES expe-riment corresponds to the kinematical region of relatively small momentum transfer where the applicability of the quark-parton model is not obvious a possible influence of mechanisms beyond QPM on the HERMES observables should be of particular interest.

## 2 Effects in the extraction of the $`\overline{d}\overline{u}`$ asymmetry

An extraction of the light quark asymmetry of the nucleon sea from charged pion multiplicities performed by HERMES  is based on the pure QPM description of the pion electroproduction. The pion production process is assumed to take place in two steps: first, a hard interaction of the incoming virtual photon with a quark in the nucleon, and second, a hadronisation of the stuck quark into final pion(s). Then, one can combine yields of positive and negative pions produced off proton and neutron, and using isospin invariance, isolate a quantity sensitive to the $`\overline{d}\overline{u}`$ asymmetry :

$$\frac{\overline{d}(x)\overline{u}(x)}{u(x)d(x)}=\frac{J(z)[1r(x,z)][1+r(x,z)]}{J(z)[1r(x,z)]+[1+r(x,z)]},$$
(1)
where $`J(z)=\frac{3}{5}\frac{1+D_{}(z)/D_+(z)}{1D_{}(z)/D_+(z)}`$ depends on the favoured $`D_+(z)`$ and unfavoured $`D_{}(z)`$ \[light quark\]-pion fragmentation functions, and $`r(x,z)=`$ $`\frac{N_p^\pi ^{}(x,z)N_n^\pi ^{}(x,z)}{N_p^{\pi ^+}(x,z)N_n^{\pi ^+}(x,z)}`$ is a ratio of differences of charged pion yields off proton and neutron.

In such an approach pions produced via other than partonic mechanisms are ignored. At the same time, experimental multiplicities $`N_{p,n}^{\pi ^+,\pi ^{}}`$, which one inputs into the right-hand side of Eq.(1) in order to obtain expe-rimental $`\overline{d}\overline{u}`$ asymmetry, do contain such nonpartonically produced pions. Thus, Eq.(1) is not an exact equation. In the $`Q^2\mathrm{}`$ limit the nonpartonic contributions vanish, but in the HERMES kinematical region where $`Q^22.5`$ GeV2 they cannot be neglected. Two mechanisms appear to be especially important: so-called VDM-like production and exclusive-$`\rho ^0`$ channel . In the VDM-like production the incoming virtual photon interacts with the nucleon strongly, via its intermediate vector meson state into which it fluctuates just before the interaction. The exclusive-$`\rho ^0`$ channel corresponds to the reaction $`\gamma ^{}+N\rho ^0+N`$ with a subsequent decay $`\rho ^0\pi ^{}\pi ^+`$.

Based on the world data for pion production in $`\pi N`$ collisions and for $`\rho ^0`$ electroproduction we have estimated contributions of these mechanisms111 For a detailed analysis see Ref..  and found them noteworthy.

The influence of the VMD and exclusive-$`\rho ^0`$ channels on the measured quantity $`\frac{\overline{d}\overline{u}}{ud}`$ can be seen in Fig.1.

The shown HERMES points were obtained with experimentally measured multiplicities $`N_{p,n}^{\pi ^+,\pi ^{}}`$. The solid curve independent of $`z`$ represents the pure QPM result calculable from both sides of Eq.(1). The two dashed curves correspond to calculations in which contributions of the nonpartonic mechanisms are included in two ways: either by simple adding them to the main partonic term, or with additional modification of the QPM term by a damping $`Q^2`$-dependent factor which is required for QPM to work properly at such a low $`Q^2`$ . A discrepancy between the solid and dashed lines gives us a theoretical error on the working formula (1). As one can see, it is comparable to or even larger then the given total (containing both statistical and systematic components) uncertainty of the final $`z`$-averaged HERMES experimental points (solid squares). The effect is stronger for the case of modified QPM term, i.e. for the more correct one. The averaging in $`z`$ also looks doubtful because of the clearly visible $`z`$-dependence.

## 3 Excess of charged over neutral pions at large $`z`$

A very interesting phenomenon has been recently observed by the HERMES Collaboration: different yields of charged and neutral pions produced off proton at large $`z`$: $`\frac{1}{2}(N_p^{\pi ^+}+N_p^\pi ^{})>N_p^{\pi ^0}`$ for $`z0.7`$ .

In the QPM picture with quark fragmentation there is no room for such an effect; there the neutral pion yield is exactly equal, due to isospin symmetry, to the charged pion yield. As a possible source for the difference, alongside the instrumental resolution, exclusive pion production channels (e.g. $`\gamma ^{}+p\pi ^++\mathrm{\Delta }^0`$) were suggested . In our opinion however, exclusive pion-production channels cannot be responsible for an effect in a broad region of $`z`$ since pions produced in this way carry almost the whole energy of the virtual photon and therefore contribute to measured multiplicities only in a narrow region near $`z1`$. In contrast, the discussed in the previous section the exclusive production of the $`\rho ^0`$ meson, which decays into two pions, gives a contribution to the charged pions multiplicity in a broad $`z`$-range. It is especially important at larger $`z`$ where the partonic rate is smaller. The contributions to the neutral pion multiplicity, caused e.g. by the exclusive-$`\rho ^\pm `$ channels, are much smaller due to smaller cross sections of these processes ($`\rho ^0`$ production, known as the dominant exclusive $`\gamma ^{}N`$ channel, is dominated by the Pomeron exchange which is absent in the case of charged-$`\rho `$ channels).

In the top panel of Fig.2 we show the pion multiplicities for $`z>0.4`$. The experimental points are from . The short-dashed curve represents the QPM prediction which is the same for both types of multiplicities. The solid curve includes also the contribution of the exclusive-$`\rho ^0`$ channel to the charged pion yield (this contribution is also shown separately). As one can see, the exclusive-$`\rho ^0`$ channel nicely explains the observed effect. This is especially well visible in the bottom panel where the ratio of charged-to-neutral pion multiplicities is shown.

## 4 Conclusions

We have discussed the role of some mechanisms beyond the quark-parton model in semi-inclusive production of pions in the kinematicl region relevant for the HERMES experiment.

We have shown that production of pions via VDM-like interaction and via exclusive-$`\rho ^0`$ channel have a strong impact on the extraction of the $`\overline{d}\overline{u}`$ asymmetry from the charged pion DIS multiplicities which was done by the HERMES Collaboration. The extracted asymmetry turns out to be very sensitive to the influence of the mentioned mechanisms and neglecting them leads to a significantly distorted result.

We have also shown that the exclusive production of $`\rho ^0`$ mesons accounts for the experimentally observed excess of charged over neutral pions produced at large $`z`$.

To conclude, nonpartonic mechanisms are essential for a good theore-tical description of the pion electroproduction in the HERMES kinematical region and cannot be omitted in a proper analysis. In addidition, we would like to stress that similar nonpartonic correction should be expected in polarised processes, as well.


text:
# A New Formation Mechanism for the Hottest Horizontal-Branch Stars

## 1. Evolution of Hot Flashers

Castellani & Castellani (1993) first showed that stars with large mass loss on the red-giant branch (RGB) can evolve to high effective temperatures before igniting helium in their cores. We have studied the properties of these “hot flashers” by following their evolution through the helium flash to the zero-age horizontal branch (ZAHB) for various amounts of mass loss (Brown et al. 2001; Sweigart et al. 2002). Two examples of our results are given in Figure 1. The “early hot flasher” in the left panel ignites helium between the tip of the RGB and the top of the white dwarf (WD) cooling curve. Mixing between the hydrogen envelope and helium core does not occur in such stars due to the high entropy of the hydrogen shell. Early hot flashers will thus have canonical, hydrogen-rich envelopes on the EHB.

The evolution is dramatically different for the “late hot flasher” in the right panel of Figure 1. Due to the weaker hydrogen shell the flash convection zone will then penetrate deeply into the stellar envelope. Protons from the envelope will be mixed into the helium core and rapidly burned, while helium and carbon from the core will be mixed outward to the stellar surface. Such flash mixing is a consequence of the basic properties of the stellar models and should therefore occur whenever a star ignites helium on the WD cooling curve.

Flash mixing has many implications for the properties of the EHB stars:

* Flash mixing produces hot EHB stars with greatly enhanced envelope helium and carbon abundances.
* Flash-mixed EHB stars will appear subluminous in the UV due to their redistributed far-UV flux and larger bolometric corrections, in agreement with the subluminous EHB stars observed in $`\omega `$ Cen and NGC 2808.
* Flash mixing leads to a sharp dichotomy in the properties of the EHB stars which may explain the gap within the EHB found in optical color-magnitude diagrams of NGC 2808.
* Flash mixing on the WD cooling curve provides a new evolutionary channel for producing hydrogen-deficient stars, particularly the high gravity, helium-rich sdO stars and the minority of sdB stars that are helium-rich.
* The high carbon abundance in the envelopes of the flash-mixed stars may be important for driving the pulsations of the variable sdB stars.

## References

Brown, T. M., Sweigart, A. V., Lanz, T., Landsman, W. B., & Hubeny, I. 2001, ApJ, 562, 368

Castellani, M., & Castellani, V. 1993, ApJ, 407, 649

Sweigart, A. V., Brown, T. M., Lanz, T., Landsman, W. B., & Hubeny, I. 2002, in $`\omega `$ Centauri: A Unique Window into Astrophysics, ed. F. van Leeuwen, J. D. Hughes, & G. Piotto (San Francisco: ASP), 261


text:
# 1 Introduction

## 1 Introduction

Though the isoscalar (IS) dipole resonance (a compression mode) was theoretically studied already more than 20 years ago  and various experimental efforts to pin down the IS giant dipole resonance (ISGDR) have been made , the distribution of the IS dipole strength obtained from the analysis of experimental data is still under dispute. The difficulty in extracting the experimental dipole strength comes from the ambiguity in the parameters of the optical potentials , which are needed in the analysis of the relevant $`(\alpha ,\alpha ^{^{}})`$ scattering data. In order to simplify our discussion, as a numerical example in the present article we discuss only the doubly magic nucleus $`{}_{82}{}^{}{}_{}{}^{208}`$Pb126. A current hot issue in both theory and experiments is whether or not a considerable amount of IS dipole strength is found in the energy region of $`E_x`$=8$``$17 MeV, which is much lower than the so-called ISGDR with a peak around $`E_x`$=23 MeV. Available calculations have a tendency to give 1-3 MeV higher peak-energy of ISGDR than that identified so far in experiments.

The calculation of IS dipole strength requires special care since the spurious translation mode has to be fully separated from intrinsic excitations. In the self-consistent random-phase-approximation (RPA) calculation the spurious mode should appear at zero energy with a finite contribution to the energy-weighted sum-rule (EWSR). The response function to the IS dipole operator

$$D=\underset{k}{}(r_k^3Y_{1\mu }(\theta _k))$$
(1)
contains the strengths of both the spurious mode and intrinsic excitations. In actual numerical calculations of self-consistent RPA the spurious mode appears not at zero energy but at either a finite energy or an imaginary energy, due to numerical inaccuracy as well as some possible practical approximations. This may lead to a question of whether or not the calculated IS dipole strength obtained in the energy region of intrinsic excitations is free from the spurious component. The difficulty in the present problem lies in the fact that a very small amount of the spurious component remaining in the energy region of intrinsic excitations can in turn produce an appreciable amount of the strength for the operator (1) due to the power of $`r^3`$. If the transition density coming from spurious component remaining in the region is proportional to $`\frac{d\rho _0}{dr}`$ where $`\rho _0`$ expresses the ground-state density, the IS dipole strength which is free from the spurious component is obtained by using the operator

$$\overline{D}=\underset{k}{}(r_k^3\eta r_k)Y_{10}(\theta _k)$$
(2)
where $`\eta `$=$`\frac{5}{3}r^2`$ , instead of using that in (1). However, if some spurious components remain at the energies of various unperturbed particle-hole (ph) excitations, the relevant transition densities are similar to those of respective original ph excitations rather than $`\frac{d\rho _0}{dr}`$. If so, using the operator in (2) may not help to take away all strength coming from the spurious component remaining in the region of intrinsic excitations.

In ref.  we studied both the IS and isovector (IV) dipole mode using the self-consistent RPA calculation with the SkM\* interaction, in which both the IS and IV correlations in RPA are taken into account solving both the Hartree-Fock (HF) and RPA equations in coordinate space using the Green’s function. Examining the IS dipole strength obtained by using the operator in (2), in the energy region of $`7<Ex<17`$ MeV we obtained an appreciable amount of IS dipole strength almost exactly at the energies of various unperturbed ph excitations. Since the residual interaction relevant to the IS dipole correlation is in general attractive, the strength increased by the RPA correlation is expected to be found at energies lower than the unperturbed ones. Therefore, using the procedure described in ref.  we tried to remove the strength at finite energies, which seemed to come from the remaining spurious component. In the continuum at a given energy the excitations of various modes with various quantum numbers may coexist. The separation of the strengths coming from modes with different quantum numbers is easy, however, we have found no established way of separating the strengths originating from different modes with the same quantum numbers. Since both the spurious mode and IS dipole resonance are isoscalar and have the same spin-parity, $`I^\pi =1^{}`$, our basic idea used in ref.  was the following. Near the isolated resonance, $`EE_n`$, the relation

$$Im\left(G_{RPA}(\stackrel{}{r},\stackrel{}{r^{}},E)\right)\rho _n^{tr}(\stackrel{}{r})\rho _n^{tr}(\stackrel{}{r^{}})$$
(3)
is known. At a high energy in the continuum, where resonances $`i`$ and $`j`$ with respective finite widths may overlap, we may write

$`Im\left(G_{RPA}(\stackrel{}{r},\stackrel{}{r^{}},E)\right)`$ $``$ $`a_i\rho _i^{tr}(\stackrel{}{r})\rho _i^{tr}(\stackrel{}{r^{}})+b_j\delta \rho _j^{tr}(\stackrel{}{r})\delta \rho _j^{tr}(\stackrel{}{r^{}})`$ (4)
$``$ $`F(\stackrel{}{r},\stackrel{}{r^{}})`$
where the modes, $`i`$ and $`j`$, have the same spin-parity and are diagonal in RPA. Choosing some $`f(\stackrel{}{r})`$ so that

$$𝑑\stackrel{}{r^{}}\rho _i^{tr}(\stackrel{}{r^{}})f(\stackrel{}{r^{}})=0$$
(5)
one obtains

$$\delta \rho _j^{tr}(\stackrel{}{r})𝑑\stackrel{}{r^{}}F(\stackrel{}{r},\stackrel{}{r^{}})f(\stackrel{}{r^{}})$$
(6)
Choosing $`f(\stackrel{}{r})=rY_{1\mu }`$ in the IS dipole case, we obtained $`\delta \rho _j^{tr}(\stackrel{}{r})`$, which was supposed to be the transition density at the energy $`E`$, that carried the spurious $`rY_{1\mu }`$ strength. Then, we subtracted $`\delta \rho _j^{tr}(\stackrel{}{r})`$ from the transition density $`\rho _{tr}(\stackrel{}{r})`$ at $`E`$ for the operator $`D`$ in (1),

$$\rho _{tr}(\stackrel{}{r})𝑑\stackrel{}{r^{}}Im\left(G_{RPA}(\stackrel{}{r},\stackrel{}{r^{}},E)\right)(r^{})^3Y_{1\mu }(\theta _k^{})$$
(7)
Using the above subtraction procedure, it seemed that in ref.  we subtracted a bit too much strength as was criticized, for example, in ref. . The oversubtraction happened, since the remaining ph excitations which carry some spurious strength were not really RPA solutions.

In all available calculations of ISGDR except those in refs.  the states in the continuum are approximated using the expansion in terms of discrete basis. Namely, the continuum wave functions are expanded in terms of harmonic oscillator basis with a maximum principal-quantum-number $`N_{max}`$. In those discrete-basis calculations the operator in (2) is commonly used in order to subtract the spurious components. The IS dipole strength obtained in the energy region of $`E_x=817`$ MeV depends very much on publications. The relativistic RPA calculation in ref.  produced especially a large amount of the IS dipole strength in the low-energy region, for all effective interactions used. In none of the discrete-basis calculations the EWSR carried by the spurious state and intrinsic excitations is carefully examined.

In the present paper we use EWSR as a measure of the validity of numerical calculations. The EWSR for the IS dipole operator in (1) is written as

$`S(r^3Y_{10})`$ $``$ $`{\displaystyle \underset{n}{}}E_nB(D;01^{}n)`$ (8)
$`=`$ $`({\displaystyle \frac{3}{4\pi }})\left({\displaystyle \frac{\mathrm{}^2A}{2M}}\right)11r^4`$
where $`E_n`$ expresses the energy of RPA states, while $`r^4`$ on the r.h.s. is evaluated with the HF ground state . The total EWSR in (8) consists of the contribution by the spurious state $`S_{spr}`$ and that by the intrinsic excitations $`S_{intr}`$, $`S=S_{spr}+S_{intr}`$. The latter can be written as

$`S_{intr}(r^3Y_{10})`$ $`=`$ $`{\displaystyle \underset{intrinsic,i}{}}E_iB(D;01^{}i)`$ (9)
$`=`$ $`({\displaystyle \frac{3}{4\pi }})\left({\displaystyle \frac{\mathrm{}^2A}{2M}}\right)\left(11r^4{\displaystyle \frac{25}{3}}r^2^2\right)`$
which is indeed equal to the EWSR for the IS dipole operator in (2).

In sect. 2 we use a harmonic oscillator model to evaluate the distribution of the EWSR of IS dipole strength among the spurious state, the lower-lying and higher-lying intrinsic excitations. This schematic model is very much simplified, nevertheless, we expect that the semi-quantitative feature of the result obtained from the model should remain unchanged in realistic calculations. In sect. 3 we present the result of the self-consistent RPA calculation with the Skyrme SkM\* interaction on 208Pb, which is the continuum calculation solving both HF and RPA in coordinate space. The result is compared with discrete-basis calculations, in which the wave functions in the continuum are expanded in terms of harmonic oscillator basis. In sect.4 conclusion and discussions are given.

## 2 Harmonic oscillator model

The distribution of IS dipole strength may be studied in terms of a simplified model  in which the particle motion is described by a harmonic oscillator potential. The particle excitations produced by the field (1) are governed by the selection rule $`\mathrm{\Delta }N`$=1 or 3 and have the energies $`\mathrm{}\omega _0`$ and $`3\mathrm{}\omega _0`$, respectively. For a single completed shell with principal quantum number $`N`$, the transition strength is given by

$$\underset{\nu _N,\nu _{N+1}}{}|\nu _{N+1}|r^3Y_{10}|\nu _N|^2=\frac{3}{4\pi }\left(\frac{\mathrm{}}{2M\omega _0}\right)^3\frac{1}{12}(N+1)(N+2)(N+3)(13N^2+52N+50)$$
(10)
$$\underset{\nu _N,\nu _{N+3}}{}|\nu _{N+3}|r^3Y_{10}|\nu _N|^2=\frac{3}{4\pi }\left(\frac{\mathrm{}}{2M\omega _0}\right)^3\frac{1}{12}(N+1)(N+2)(N+3)(N+4)(N+5)$$
(11)
where $`\nu _N`$ represents the quantum numbers needed to specify the single-particle states in the shell $`N`$. For simplicity, we do not include the spin-isospin degeneracy in the present model. Since the last three occupied shells contribute to the total $`\mathrm{\Delta }N`$=3 strength we obtain

$$\underset{N=N_F2}{\overset{N_F}{}}\underset{\nu _N,\nu _{N+3}}{}|\nu _{N+3}|r^3Y_{10}|\nu _N|^2=\frac{3}{4\pi }\left(\frac{\mathrm{}}{2M\omega _0}\right)^3\frac{1}{12}(N_F+1)(N_F+2)(N_F+3)(3N_F^2+12N_F+20)$$
(12)
while only the last-filled $`N=N_F`$ shell contributes to the $`\mathrm{\Delta }N`$=1 strength.

Now, the major part of the $`\mathrm{\Delta }N`$=1 strength, which is obtained by setting $`N=N_F`$ in (10), corresponds to the spurious excitation, namely center of mass motion, while in the present harmonic oscillator model the $`\mathrm{\Delta }N`$=3 strength contains no spurious component. Subtracting the strength coming from the spurious excitation

$$\frac{3}{4\pi }\left(\frac{\mathrm{}}{2M\omega _0}\right)^3\frac{25}{24}(N_F+1)(N_F+2)^3(N_F+3),$$
(13)
from the expression (10) with $`N=N_F`$, we obtain the strength of the intrinsic $`\mathrm{\Delta }N`$=1 excitation

$$\frac{3}{4\pi }\left(\frac{\mathrm{}}{2M\omega _0}\right)^3\frac{1}{24}N_F(N_F+1)(N_F+2)(N_F+3)(N_F+4)$$
(14)
We note that the total number of particles in the present model is written as

$$\frac{1}{6}(N_F+1)(N_F+2)(N_F+3)$$
(15)
while

$$\underset{N=0}{\overset{N_F}{}}r^2=\left(\frac{\mathrm{}}{2M\omega _0}\right)\frac{1}{4}(N_F+1)(N_F+2)^2(N_F+3)$$
(16)
and

$$\underset{N=0}{\overset{N_F}{}}r^4=\left(\frac{\mathrm{}}{2M\omega _0}\right)^2\frac{1}{2}(N_F+1)(N_F+2)(N_F+3)(N_F^2+4N_F+5)$$
(17)
The spurious contribution (13) is estimated by recognizing that the spurious component is found only at the $`\mathrm{\Delta }N`$=1 excitations and that the transition density of the spurious excitation is proportional to the radial derivative of the ground-state density, $`d\rho _0/dr`$.

From the expressions, (12) and (14), one obtains the ratio of the unperturbed transition strength of $`\mathrm{\Delta }N`$=1 excitations to that of $`\mathrm{\Delta }N`$=3 ones to be 1:6 for $`N_F1`$. Correspondingly, the ratio of the unperturbed EWSR is 1:18.

Introducing the residual interaction of a separable form

$$V=cD^{}D$$
(18)
where $`c`$ expresses an attractive coupling constant, we solve the RPA equation and estimate the distribution of the transition strength together with the RPA energies. The following results of cases (a) and (b) are summarized in Table I.
(a) If we choose $`c`$ so that the higher-lying RPA solution is obtained at (2.5)$`\mathrm{}\omega _0`$, the lower-lying RPA solution is found at (0.872)$`\mathrm{}\omega _0`$. Then, the RPA transition strength is distributed with 1.09:2.79 over the lower-lying and higher-lying RPA solutions, while the distribution of the EWSR is 0.95:6.98 .
(b) If we choose $`c`$ so that the higher-lying RPA solution is obtained at (2.8)$`\mathrm{}\omega _0`$, the lower-lying RPA solution is found at (0.961)$`\mathrm{}\omega _0`$. Then, the RPA transition strength is distributed with 0.60:2.62 between the lower-lying and higher-lying RPA solutions, while the distribution of the EWSR is 0.58:7.34 .

Comparing the peak energies of the RPA solutions with those obtained by available self-consistent HF plus RPA calculations, one expects that the realistic situation may lie somewhere between the case (a) and (b). Before the spurious component $`\delta \rho _j^{tr}(\stackrel{}{r})`$ in (6) was subtracted in our self-consistent RPA calculations of ref. , a number of lower-energy peaks at Ex $`\stackrel{<}{}`$ 17 MeV for the operator (2) were found nearly at respective unperturbed ph energies. In contrast, the higher-energy peak around 25 MeV was collectively constructed by shifting the unperturbed strength coherently from the region of higher energy. In other available self-consistent RPA calculations, for example in refs.  and  in which the expansion in terms of harmonic oscillator basis was used, or in ref.  which was a continuum calculation of relativistic Hartree plus RPA calculation, a considerable amount of low-energy strength appeared also around unperturbed ph energies. Compared with those HF plus RPA calculations, in the case (a) of Table I the higher-lying RPA peak is pushed down by about the same order of magnitude and the lower-lying peak is too strongly pushed down, while in (b) the higher-lying peak is appreciably less pushed down and the lower-lying peak is still slightly lower.

On the other hand, one might wonder whether or not the ratio of the coupling between the $`\mathrm{\Delta }N`$=1 and 3 excitations to that within the excitations of a given $`\mathrm{\Delta }N`$ is realistic in the present schematic interaction of (18). In the present harmonic-oscillator model the effective strength of the coupling between the $`\mathrm{\Delta }N`$=1 and 3 excitations can be measured from the increase of EWSR in the lower-lying RPA solution compared with that in the unperturbed $`\mathrm{\Delta }N`$=1 excitation. The portion of the EWSR in (8) carried by the lower-energy peak is increased from 2.3 to 5.2 percent in the case (a) and to 3.2 percent in the case (b), as shown in Table I.

## 3 Hartree-Fock and RPA calculations with Skyrme interactions

In the HF calculation with the SkM\* interaction for 208Pb we obtain

$`r^2`$ $`=`$ $`30.85fm^2`$
$`r^4`$ $`=`$ $`1218.2fm^4`$ (19)
Thus, using eqs. (8) and (9) the EWSR becomes

$$S(r^3Y_{10})=\mathrm{\hspace{0.33em}13.79}\times \mathrm{\hspace{0.17em}10}^6fm^6MeV$$
(20)
$$S_{intr}(r^3Y_{10})=\mathrm{\hspace{0.33em}5.632}\times \mathrm{\hspace{0.17em}10}^6fm^6MeV$$
(21)
From eqs. (20) and (21) the contribution by the spurious state to the EWSR for the operator $`r^3Y_{1\mu }`$ is 8.16 $`\times `$ 106 $`fm^6MeV`$, which is 59.2 percent of the total EWSR in (20).

In the present work we perform the self-consistent RPA calculation with the SkM\* interaction, in which both the IS and IV correlations in RPA are taken into account solving both HF and RPA equations in coordinate space using the Green’s function. The difference of the present calculation from that in ref.  is that we use the radial mesh ($`\mathrm{\Delta }r`$) of 0.1 $`fm`$ in both the HF and RPA calculations and, moreover, we carefully check the sum-rules consumed by both intrinsic and spurious excitations. We present the numerical result without adopting the subtraction procedure, which was used in ref.  and was also explained in the Introduction. In ref.  $`\mathrm{\Delta }r`$=0.1 $`fm`$ was used in HF, while we adopted $`\mathrm{\Delta }r`$=0.3 $`fm`$ in RPA. In Fig. 1 we show the RPA response functions, which are smeared out using the width of 0.5 MeV, for the operator (1) obtained with $`\mathrm{\Delta }r`$=0.1 $`fm`$ and 0.3 $`fm`$ and for the operator (2) calculated with $`\mathrm{\Delta }r`$=0.1 $`fm`$. A great, somewhat unexpected improvement in the result of the present calculation with a finer radial mesh in RPA is : the appreciable difference between the response functions to the operators (1) and (2) in the energy region of the higher peak around 25 MeV is drastically decreased. Moreover, the spurious state appears now slightly below 0 MeV (namely, at an imaginary energy) and consumes a much larger portion of EWSR, $`S_{spr}`$.

As in ref. , our effective interaction used in RPA consists only of the part in the spin-independent channels. Using discrete-basis calculations, in which the inclusion of the full Skyrme effective interaction is straightforward and easy, we have checked that the RPA response function of the present IS dipole operator changes very little when we use the full effective interaction, as seen in Tables II and III.

In Table II we show the transition strength for the operator (1), the contribution to EWSR, and the percent of the total EWSR, $`S(r^3Y_{10})`$, at the spurious state and in the region of $`2<E_x<17`$ MeV, $`17<E_x<70`$ MeV, and $`E_x>70`$ MeV. In Table III the same quantities are tabulated for the operator (2). The borderline 17 MeV is chosen, since almost all $`\mathrm{\Delta }N=1,I^\pi =1^{}`$ ph excitations in the HF potential lie below 17 MeV and no unperturbed 1- ph excitations are found for $`18<E_x<24`$ MeV.

As is previously known, it is seen from Table II that the portion of the EWSR found in the energy region of intrinsic excitations ($`E_x>2`$ MeV) clearly shows the presence of the spurious components remaining in the region. From Table III it is seen that the ratio of the EWSR lying in the region of $`E_x<17`$ MeV to that for $`E_x>17`$ MeV is 1:4.9. More importantly, the portion of the EWSR obtained for $`E_x>2`$ MeV is 39.1 percent of the total EWSR, compared with the full value, 40.8 percent. It should be mentioned that in the continuum calculation it is technically very difficult to obtain accurate values of the transition strength carried by very sharp peaks including the spurious state, in contrast to the discrete-basis calculations. Furthermore, in the present continuum calculation the spurious state appears slightly below 0 MeV. Thus, reducing the effective ph interactions by multiplying a factor of 0.9902, we obtain the spurious state at $`E_x`$=0.107 MeV, from which we estimate the carried transition strength. Though it is clear that the continuum calculations performed using a much finer radial mesh than 0.1 $`fm`$ would further improve the accuracy, here we do not pursue it any further.

In Fig. 2 we show the calculated RPA response functions to the operators, (1) and (2), as a function of excitation energy. Though the peaks at $`E_x>`$ 7.5 MeV lie in the continuum, calculated lower-lying peaks are too sharp to be plotted. Since the spreading width is anyhow not included in the RPA response function, in Fig. 2 we plot the quantities, which are obtained by smearing out the calculated RPA response function with the width of 1 MeV. It is seen that an appreciable amount of IS dipole strength is present in the energy region of $`E_x`$=7-17 MeV, which is much lower than the so-called ISGDR around 25 MeV. The presence of the lower-lying IS dipole strength indicates that the ISGDR appearing around 25 MeV is not really a good collective mode in the nuclear system.

For reference, in Tables II and III we show also the result of discrete-basis calculations for the operators, (1) and (2). How to truncate the space of harmonic oscillator basis seems to depend on publications. What we have adopted in the present work, which is along the line of ref. , is the following. First we construct the HF potential, which is calculated in coordinate space and, then, calculate the one-particle energies and wave-functions for the potential using the harmonic oscillator basis with the radial node $`n^{HO}n_{max}^{HO}`$, where $`n_{max}^{HO}`$ is taken to be independent of $`\mathrm{}`$. Then, the RPA equation is solved including the lowest-lying $`\frac{1}{2}(N_{max}^{HF}\mathrm{})`$ one-particle levels for a given $`\mathrm{}`$, where $`N_{max}^{HF}`$ is taken to be independent of $`\mathrm{}`$. This truncation of the ph space included in RPA is approximately the same as that limiting the ph space by the $`\mathrm{}`$-independent maximum energy of particle configurations. It is noted that, as $`n_{max}^{HO}`$ increases, the bound one-particle energies and wave-functions approach the eigen energies and eigen functions of the HF potential, respectively, while the properties of the positive-energy one-particle levels can be far away from the correct ones unless $`n_{max}^{HO}`$ is taken to be $`\mathrm{}`$. In actual calculations one must take finite values of $`n_{max}^{HO}`$ and $`N_{max}^{HF}`$ and, thus, it is not guaranteed that using finite larger values of $`n_{max}^{HO}`$ and $`N_{max}^{HF}`$ produces a more accurate result of IS dipole strength.

For the parameter set, $`n_{max}^{HO}=12`$ and $`N_{max}^{HF}=24`$, in Tables II and III we show both the result obtained by using the full effective interaction in RPA and the one by employing the same effective interaction used in our continuum RPA calculation. The fact that very little difference is found between these two results indicates that the IS dipole strength obtained from our continuum RPA calculation is reliable. Comparing the numbers in the last column of Table II with those in Table III, in the discrete-basis calculation the difference between the strengths calculated with the operators (1) and (2) is relatively minor in the region of $`E_x>17`$ MeV. On the other hand, from Table II we notice that the calculated energy of the spurious state is pretty high compared with the continuum calculation and, furthermore, only 92 percent of the correct portion of EWSR is consumed by the calculated spurious state. The properties of the spurious state converge surprisingly slowly to the right ones, as we increase $`n_{max}^{HO}`$ and $`N_{max}^{HF}`$. On the other hand, if we use a value of $`n_{max}^{HO}`$ larger than 12, the calculated response function in the continuum starts to behave strangely instead of converging to some reasonable result. This strange behavior of discrete-basis calculations seems to be known for experts in the fields .

## 4 Conclusion and discussions

First, the distribution of IS dipole strength for heavy nuclei is studied using the schematic harmonic-oscillator model, in which the strength is distributed among the spurious mode, the low-energy and high-energy peaks. Though 96 percent of the unperturbed $`\mathrm{\Delta }N=1`$ strength for the operator $`r^3Y_{1\mu }`$ corresponds to the spurious excitation, the remaining 4 percent contributes to the intrinsic excitations. The spurious state consumes 56.8 percent of the total EWSR for the operator $`r^3Y_{1\mu }`$, while the low-energy RPA peak may carry up till 5 percent of the EWSR.

Secondly, the result of IS dipole strength in 208Pb is presented, which is calculated using the self-consistent HF plus RPA solved in coordinate space using the Green’s function method, with better accuracy than in ref. . If one performs numerical calculations with very good accuracy, the IS dipole strength of intrinsic excitations should not depend on whether the operator (1) or (2) is used. At the same time, the spurious mode should appear close to zero energy and consume the whole portion of EWSR expected. In our present continuum calculation we could not attain such a satisfactory accuracy. On the other hand, if the calculated strengths with the operators (1) and (2) are appreciably different from each other in the energy region of intrinsic excitations, the spurious component remaining in the region may not have the transition density proportional to $`\frac{d\rho _0(r)}{dr}`$. In the present work we have not attempted to take away this kind of possible spurious component from the region of intrinsic excitations, since we find no reliable way of excluding it.

Using the radial mesh of 0.1 $`fm`$ in both the HF and RPA calculations, the spurious state is obtained slightly below 0 MeV. Multiplying the effective ph interaction by a factor of 0.9902, we push up the energy to 0.107 MeV and estimate the EWSR carried by the state, which turns out to be the major portion of the expected EWSR. On the other hand, the EWSR carried by the calculated states for $`E_x>2`$ MeV is 114 and 96 percent of the expected EWSR for the operators (1) and (2), respectively. The accuracy of the continuum calculations is appreciably improved by decreasing $`\mathrm{\Delta }r`$ used in solving RPA from 0.3 $`fm`$ in ref.  to 0.1 $`fm`$, however, it is still away from being very satisfactory. In the continuum calculations the accuracy of separating the spurious component from intrinsic excitations is surprisingly sensitive to the size of radial mesh employed in RPA. Thus, the response function of IS dipole strength has to be calculated using a very fine radial mesh, in order to obtain a reliable numerical result. From the present study we find that slightly more than 5 percent of the total EWSR is expected as intrinsic excitations lying in the low-energy region below $`E_x<17`$ MeV.

The expansion of HF and RPA solutions in terms of harmonic oscillator basis can work for bound states if one includes a sufficiently large number of bases . Thus, in the present study of IS dipole strength the properties of the spurious mode become better and better as the number of bases increases, though we find that the convergence of the properties is unexpectedly slow as $`n_{max}^{HO}`$ and $`N_{max}^{HF}`$ increase. On the portion of EWSR found in the continuum the discrete-basis calculation presented in Table II and Fig. 1 produces 9 percent larger in the $`2<E_x<17`$ MeV region and 6 percent smaller for $`E_x>17`$ MeV, compared with our continuum calculation. However, the real problem is : If we further enlarge the space of discrete-basis, the calculation starts to produce more than 100 percent of EWSR, which are certainly unreasonable. Thus, it is not clear whether or not the discrete-basis calculations can be reliably used for the estimate of IS dipole strength.

The authors are grateful to Drs. N. Van Giai and G. Colo for the generous advice on the discrete-basis calculation. I.H. acknowledges the financial supports provided by Crafoordska Stiftelsen and RIKEN. H.S. acknowledges financial support by the Japanese Ministry of Education, Science, Sports and Culture by Grant-In-Aid for Scientific Research under the program number C(2) 12640284.

Figure captions

RPA response functions to the IS dipole operators, (1) and (2), as a function of excitation energy, which are obtained from the self-consistent HF plus RPA calculations solved in coordinate space. A radial mesh $`\mathrm{\Delta }r`$=0.1 $`fm`$ is used in HF, while the results obtained by using $`\mathrm{\Delta }r`$=0.1 and 0.3 $`fm`$ in RPA are compared. We show the calculated response functions, which are smeared out using the width of 0.5 MeV.

RPA response functions to the IS dipole operators, (1) and (2), as a function of excitation energy, which are obtained from the self-consistent continuum calculation (CRPA) with $`\mathrm{\Delta }r`$=0.1 $`fm`$ and the discrete-basis calculation with $`n_{max}^{HO}=12`$ and $`N_{max}^{HF}=24`$ (DRPA). In all calculations the (1, $`\tau \tau `$) effective interaction (see the caption to Table II) is used. In order to see more easily the total distribution of the strength, we plot both CRPA and DRPA response functions, which are smeared out using the width of 1 MeV.