text:
# Neutrinos Associated With Cosmic Rays of Top-Down Origin

## I Introduction

The discovery of cosmic rays with energy exceeding the GZK cutoff gzk  presents an interesting challenge to astrophysics, particle physics, or both agasa ; hires . Numerous scenarios have been proposed to solve the problem. These include exotic particles exo , neutrinos with QCD scale cross sections qcdneu , semi-local astrophysical sources loc  and top-down models top .

Top-down models can be motivated by a variety of arguments. For example, the recent measurements of the cosmic microwave background and of supernova redshifts have dramatically confirmed that our universe contains a large fraction of cold dark matter cmb . A top-down model in which annihilating or decaying superheavy relic particles produce the highest energy cosmic rays could potentially solve both of these problems wim1 ; wim . Topological defects could also solve the ultra-high energy cosmic ray problem in a similar way top . In this paper, we reevaluate the implications of generic top-down models producing parton jets of $`10^{21}`$ or $`10^{25}`$ eV that fragment into the observed super GZK cosmic rays.

Established particle physics implies that such ultra high-energy jets fragment predominantly into pions and kaons, with a small admixture of protons lep . The mesons will eventually decay into photons or electrons plus neutrinos. A typical QCD jet therefore produces more photons than protons. This is true in particular at relatively low values of $`x=E_{\mathrm{particle}}/E_{\mathrm{jet}}`$, but even at large $`x`$ the photon flux is at least as large as the proton flux in a jet. This seems to be in disagreement with mounting evidence that the highest energy cosmic rays are not photons pro1 . The observed shower profile of the original Fly’s Eye event, with energy exceeding $`10^{20}`$ eV, fits the assumption of a primary proton, or, possibly, that of a nucleus. The shower profile information is sufficient to conclude that the event is unlikely to be of photon origin. The same conclusion is reached for the Yakutsk event that is characterized by a large number of secondary muons, inconsistent with a purely electromagnetic cascade initiated by a gamma ray. A recent reanalysis of Haverah Park data further reinforces this conclusion park . In light of this information, it seems likely that protons, and not gamma rays, dominate the highest energy cosmic ray spectrum. This does not necessarily rule out superheavy particles as the source of the highest energy cosmic rays. The uncertainties associated with the cascading of the jets in the universal and galactic radio backgrounds and with the strength of intergalactic magnetic fields leave open the possibility that ultra high-energy photons may be depleted from the cosmic ray spectrum near $`10^{20}`$ eV, leaving a dominant proton component at GZK energies lee ; radio . With this in mind, we will choose to normalize the proton spectrum from top-down scenarios to the observed ultra high-energy cosmic ray flux.

Neutrinos are produced more numerously than protons and travel much greater distances. The main point of this paper is to point out that this “renormalization” of the observed cosmic ray flux to protons generically predicts observable neutrino signals in operating experiments such as AMANDA II, RICE and AGASA. Top-down models, if not revealed, will be severely constrained by high-energy neutrino observations in the near future.

## II Nucleons from Ultra-high energy Jets

The assumption that nucleons from the decay (or annihilation) of very massive $`X`$ particles are the source of the highest energy cosmic rays normalizes the decay or annihilation rate of their sources, once the shape of the spectrum of the produced nucleons is known. One needs mass $`M_X10^{21}`$ eV in order to explain the observed UHECR events. The presence of such very massive particles strongly indicates the existence of superparticles with masses at or below the TeV scale, since otherwise it would be difficult to keep the weak energy scale ten or more orders of magnitude below $`M_X`$ in the presence of quantum corrections. Moreover, we know that all gauge interactions are of comparable strength at energies near $`M_X`$. These two facts together imply that the evolution of a jet with energy $`10^{21}`$ eV shows some new features not present in jets produced at current particle collider experiments.

First of all, primary $`X`$ decays are likely to produce approximately equal numbers of particles and superparticles, since $`M_X`$ is much larger than the scale $`M_{\mathrm{SUSY}}1`$ TeV of typical superparticle masses. Even if the primary $`X`$ decay only produces ordinary particles, superparticles will be produced in the subsequent shower evolution. Note also that (at least at high energies) electroweak interactions should be included when modeling the parton shower. Both effects taken together imply that the jet will include many massive particles – superparticles, electroweak gauge and Higgs bosons, and also top quarks. The decays of these massive particles increase the overall particle multiplicity of the jet, and also produce quite energetic neutrinos, charged leptons and lightest supersymmetric particles (LSP). Eventually the quarks and gluons in the jet will hadronize into baryons and mesons, many of which will in turn decay.

We model these jets at the point of their origin using the program described in Ref.bd1 . This program allows us to calculate spectra for different $`X`$ decay modes. It then follows the supersymmetric parton cascade down to virtuality (or inverse time) of the order of $`M_{\mathrm{SUSY}}`$, including all gauge interactions as well as third generation Yukawa interactions. At $`M_{\mathrm{SUSY}}`$ all massive particles are decoupled from the parton shower, and decay. Supersymmetric cascade decays are fully taken into account; the results presented below have been obtained using the same spectrum of superparticles as in ref.bd1 . At virtualities below $`M_{\mathrm{SUSY}}`$ only ordinary QCD interactions contribute significantly to the development at the jet; $`b`$ and $`c`$ quarks are decoupled at their respective masses, hadronize, and decay. At a virtuality near $`1`$ GeV the light quarks and gluons hadronize, with a meson to baryon ratio of roughly thirty to one (five to one) at small (large) $`x`$. All baryons will eventually decay into protons, while the mesons (mostly pions) decay into photons, electrons footnote  and neutrinos (plus their antiparticles). The heavier charged leptons (muons and taus) also decay. The final output of the code is the spectra of seven types of particles which are sufficiently long–lived to reach the Earth: protons, electrons, photons, three flavors of neutrinos, and LSPs. We assume that $`X`$ decays are CP–symmetric, i.e. we assume equal fluxes of particles and antiparticles of a given species.

The calculation of Ref.bd1  was based on conventional one–loop evolution equations for the relevant fragmentation functions. These may not be reliable in the region of very small $`x`$. We wish to calculate neutrino fluxes at energies down to $`10^{15}`$ eV (1 PeV), which corresponds to $`x10^6(10^{10})`$ for $`M_X=10^{21}(10^{25})`$ eV. At these very small $`x`$ values color coherence effects are expected to suppress the shower evolution mlla . We try to estimate the size of these effects by matching our spectra computed using conventional evolution equations to the so–called asymptotic MLLA spectra; details of this procedure will be described elsewhere bd2 . The effect of this modification on the neutrino event rate is relatively modest for primary jet energy near $`10^{21}`$ eV, but becomes significant at $`10^{25}`$ eV. However, even at this higher energy the proton flux, which we only need at $`x10^5`$, is not affected significantly.

This calculation gives us the shape of the spectra of the stable particles at source. The spectra on Earth might differ significantly due to propagation effects. As stated in the introduction, we will assume that (almost) all UHE photons get absorbed. This is actually expected to be true for a homogeneous source distribution. However, according to current estimates of the strengths of the magnetic fields and of the radio wave background in (the halo of) our own galaxy most UHE photons produced in the halo of our galaxy are expected to reach the Earth. As stated in the Introduction, this seems to be in conflict with observation. We will therefore assume that the interaction length of UHE photons in our galaxy has been greatly over–estimated, and explore the consequences of this assumption for neutrino signals.

As well known, (anti)protons lose energy when traveling through the intergalactic medium, mostly through scattering off photons of the ubiquitous cosmic microwave background (CMB). We calculate the observed spectrum of protons taking into account scattering off the CMB at the $`\mathrm{\Delta }`$resonance and scattering by $`e^+e^{}`$ pair production; energy losses through the Hubble expansion of the Universe are also included lee ; prop . Note that the photoproduction of charged pions contributes to the observed neutrino flux on Earth. In order to solve the ultra high-energy cosmic ray problem, the (anti)proton flux must accommodate the events above the GZK cutoff. Observations indicate on the order of a few times $`10^{27}`$ events $`\mathrm{m}^2\mathrm{s}^1\mathrm{sr}^1\mathrm{GeV}^1`$ in the energy range above the GZK cutoff ($`5\times 10^{19}`$ eV to $`2\times 10^{20}`$ eV)hires ; agasa .

The formalism of a generic top-down scenario is sufficiently flexible to explain the data from either the HIRES hires  or AGASA agasa  experiments. Figure 1 compares HIRES and AGASA data to the proton spectrum predicted for a galactic distribution of decaying particles with mass $`M_X=210^{21}`$ eV. The drop near a few times $`10^{19}`$ eV is a manifestation of the GZK cutoff. Note, however, that there are sufficient semi-local events to explain all observed super GZK events. Similarly, figure 2 compares HIRES and AGASA data to the spectrum predicted for $`M_X=210^{25}`$ eV, rather than $`210^{21}`$ eV, decaying particles for the same distribution. Although HIRES and AGASA data differ at face value, especially above the GZK cutoff, top-down scenarios can accommodate all events observed above the GZK cutoff in either experiment.

If the cosmic ray sources are not distributed with a large overdensity in the galaxy, the resulting cosmic ray and neutrino spectrum will be modified. For example, using a homogeneous distribution, the GZK cutoff will again be manifest and the observed cosmic ray spectrum will be difficult to explain. A galactic overdensity of $`10^3`$ to $`10^4`$ or more seems necessary to fit the data. The figure 1 shows a $`10^5`$ overdensity, which is the overall overdensity of matter in our galaxy at the location of the Sun. Note that for less extreme overdensities, the average distance at which a proton is produced will be larger. This implies larger energy losses, and hence a reduced proton flux on Earth for a given number of sources. Conversely, if we fix the proton flux to the observed flux of UHECR events, models with lower overdensity require more sources. Since neutrino fluxes are not degraded by propagation through the intergalactic medium, the number of neutrinos increases proportionally to the number of sources, with additional contributions to the neutrino flux coming from pion production on the CMB background. Thus, the neutrino event rates and spectrum shown in the figures reflect the most conservative choice of distributions. The table shows results for both homogeneous and galactic distributions.

## III Neutrinos from Ultra high-energy Jets

As discussed earlier, the program computing the proton flux at source also gives the neutrino flux at source. Neutrinos, not being limited by scattering, travel up to the age of the universe at the speed of light ($``$ 3000 Mpc in an Euclidean approximation). The only nontrivial effect of neutrino propagation is due to oscillations. In our case the propagation distance of neutrinos amounts to many oscillation lengths, if oscillation parameters are fixed by the currently most plausible solutions of the atmospheric and solar neutrino deficits bahcall . As a result, the UHE neutrino flux on Earth is the same for all three flavors, and amounts to the average of the fluxes of the three neutrinos flavors at source.

The predicted neutrino flux is shown in figures 3 and 4. At $`E_\nu E_{\mathrm{jet}}`$ the main contribution comes from $`\pi ^\pm \mu ^\pm \nu _\mu e^\pm \nu _e\nu _\mu `$ decays, but at larger $`E_\nu `$ there can be significant contributions from the decays of heavy (s)particles. The peak in the dotted curves at $`E_\nu =E_{\mathrm{jet}}`$ results from our assumption that in this scenario $`X`$ decays directly into first or second generation $`SU(2)`$ doublet (s)leptons, which implies that 50% of all $`X`$ decays give rise to a primary neutrino; in this case the ratio of neutrino and proton fluxes has a maximum at high energy. On the other hand, if primary $`X`$ decays are purely hadronic, the neutrino flux at the largest energy is only slightly above the proton flux at that energy. The reason is that neutrinos from meson decays only carry a fraction of the energy of the meson, so a five to one meson to proton ratio at large $`x`$ leads to a nearly one to one neutrino to proton ratio. We see that the neutrino flux at the highest energy depends quite strongly on how the $`X`$ particles decay; there is also some dependence on the parameters of the SUSY model bd1 ; bd2 . For given proton flux the neutrino flux at smaller $`x`$ is much less model dependent. At very small $`x`$ a new uncertainty appears due to coherence effects. These have so far only been studied in a pure QCD parton shower; our treatment of these effects is therefore of necessity rather crude.

## IV Event Rates in High-energy Neutrino Telescopes and Air Shower Experiments

We will discuss two classes of experiments capable of observing high energy cosmic neutrinos: neutrino telescopes and air shower experiments.

Optical Cerenkov neutrino telescopes such as the operating AMANDA II and next generation IceCube are designed to observe muon tracks from charged current interactions as well as showers which occur in the detector. The probability of detecting a neutrino passing through the detector from its muon track is given by

$$P_{\nu \mu }(E_\nu ,\theta _{\mathrm{zenith}})=\sigma _{\nu N}(E_\nu )n_{\mathrm{H}_2\mathrm{O}}R_\mu (E_\mu ,\theta _{\mathrm{zenith}})$$
(1)
where $`n_{\mathrm{H}_2\mathrm{O}}`$ is the number density of nucleons in the detector medium (water or ice), and the muon range $`R_\mu (E_\mu ,\theta _{\mathrm{zenith}})`$ is the average distance traveled by a muon of energy $`E_\mu `$ before falling below some threshold energy (we have used 100 TeV). This quantity depends on the zenith angle of the incoming neutrino because for a detector depth of $`2`$km, only quasi-horizontal or upgoing events can benefit from longer muon ranges. At the energies we are most concerned with, the majority of muon events will be quasi-horizontal. The number of muon events observed is then given by

$$N_{\mathrm{events}}=𝑑E_\nu 𝑑\mathrm{\Omega }\frac{d\varphi _\nu }{dE_\nu }P_{\nu \mu }(E_\nu ,\theta _{\mathrm{zenith}})A_{\mathrm{eff}}T,$$
(2)
where $`T`$ is the time observed and $`A_{\mathrm{eff}}`$ is the effective area of the detector: one twentieth square kilometers for AMANDA II and one square kilometer for IceCube.

AMANDA II and IceCube can also observe showers generated in charged or neutral current interactions within the detector volume. The event rate from showers is not enhanced by long muon ranges, but can be generated by all three flavors of neutrinos and with greater cross section (neutral + charged current). We use a shower energy threshold of 100 TeV. The energy threshold imposed effectively removes any background events from atmospheric neutrino events. For a review of Optical Cerenkov neutrino telescopes see Ref.ice2 .

The operating radio Cerenkov experiment, RICE, is capable of observing showers generated in charged current electron neutrino events. RICE’s effective volume increases with energy. At 100 TeV, RICE has an effective volume less than one hundredth of a cubic kilometer. By 10 PeV, however, it increases to about ten cubic kilometers rice . Again, we use a hard 100 TeV shower threshold.

Air shower experiments can also observe very high energy cosmic neutrinos. We consider AGASA, the largest ground array currently in operation agasa2 , and the next generation AUGER array auger .

To determine that an air shower was initiated by a neutrino, rather than a proton or other cosmic ray, we require a slant depth greater than 4000 $`\mathrm{g}/\mathrm{cm}^2`$. This corresponds to a zenith angle very near 75 degrees. Therefore, only quasi-horizontal air shower events can be identified as neutrinos. Additionally, unlike showers generated in the upper atmosphere, deeply penetrating showers provide both muon and electromagnetic shower components which help them be differentiated from showers with hadronic primaries. The probability of detecting and identifying a neutrino initiated air shower is described in terms of the array’s acceptance, $`A`$, in units of volume times water equivalent steradians (we sr). The detector’s acceptance increases with energy. For AGASA, the acceptance is about $`0.01\mathrm{km}^3\mathrm{we}\mathrm{st}`$ at $`10^7`$ GeV but increases to $`1.0\mathrm{km}^3\mathrm{we}\mathrm{st}`$ at $`10^{10}`$ GeV and above. For AUGER, the acceptance is about $`0.1\mathrm{km}^3\mathrm{we}\mathrm{st}`$ at $`10^7`$ GeV, $`10.0\mathrm{km}^3\mathrm{we}\mathrm{st}`$ at $`10^9`$ GeV and $`50.0\mathrm{km}^3\mathrm{we}\mathrm{st}`$ at $`10^{12}`$ GeV. The number of events observed is then

$$N_{\mathrm{events}}=𝑑E_\nu 𝑑\mathrm{\Omega }n_{\mathrm{H}_2\mathrm{O}}\frac{d\varphi _\nu }{dE_\nu }\sigma _{\nu N}(E_\nu )A(E_\nu )T,$$
(3)
where $`T`$ is again the time observed, $`n_{\mathrm{H}_2\mathrm{O}}`$ is the number density of nucleons in water and $`A(E_\nu )`$ is the detector’s acceptance. AGASA presently has about five years of effective running time between 1995 and 2000 analyzed. A useful treatment of air shower event rates from neutrinos can be found in Ref. feng .

Table 1 shows the event rates expected for a variety of models, and for several experiments. AMANDA-II, with an effective area of $``$50,000 square meters can place the strongest limits on high energy neutrino flux presently. Furthermore, AGASA, with five years of effective observing time, has similar sensitivity. RICE, just beginning to release results, will be capable of raising the level to which top-down scenarios can be tested, perhaps being capable of testing all 16 models shown in the table. Even if no events are observed with operating experiments, next generation experiments, especially IceCube, will be able to test all models with adequate sensitivity.

Event rates shown in table 1 include only events below $`10^{12}`$ GeV. Above this energy, uncertainties in the neutrino-nucleon cross sections and in detector performance make such calculations difficult and unreliable. Our most reasonable extrapolations into this energy range indicate about a 20% enhancement to the event rate if all energies are considered for $`10^{25}`$ eV jets. There is no effect for the $`10^{21}`$ eV jet case.

High-energy neutrino event rates have been calculated in Ref. jaime  for a similar model. Their calculation used the model of reference bere  which normalized the ultra high-energy cosmic ray flux to the photons and protons generated in superheavy particle decay rather than the proton flux alone. For this reason, their results show only two events per year in a square kilometer neutrino telescope, a smaller rate than we predict for most models. Another recent estimate of neutrino fluxes on Earth in top–down models kalashev  finds broadly similar results as our’s. However, there the ‘MLLA’ form for the fragmentation functions was used for all energies, which (incorrectly) predicts nearly energy–independent ratios of neutrino, photon and proton fluxes.

## V Conclusions

If a top-down scenario, such as the decay or the annihilation of superheavy relics, is the source of the highest energy cosmic rays, then a high-energy neutrino flux should accompany the observed cosmic ray flux. This neutrino flux will be much higher than the flux of nucleons due to the much greater mean free path of neutrinos and greater multiplicity of neutrinos produced in high-energy hadronic jets.

The high-energy neutrino flux generated in such a scenario can be calculated by normalizing the flux of appropriate particles to the ultra–high energy cosmic ray flux. With mounting evidence that the highest energy cosmic rays are protons or nuclei and not photons, we have assumed that the ultra high-energy photons are degraded by the universal and/or galactic radio background, leaving protons to dominate the highest energy cosmic ray flux. The neutrino flux must then be normalized to the proton flux resulting in significantly improved prospects for its detection.

A word about the uncertainties in our calculation might be in order. First of all, the uncertainty of the measured UHECR flux, and in particular the discrepancy between the HIRES and AGASA results, leads to an overall uncertainty of a factor of $`23`$. On the theoretical side, the main uncertainty probably comes from the calculation of the particle spectra at “small” energies, where currently not very well understood coherence effects can play a role. This effect is bigger for higher primary jet energy, and can change the event rate by up to a factor of about 7 (see table). Relaxing our assumption that all UHE photons are absorbed would lead to a corresponding reduction of the fitted source density, and hence of the neutrino flux. In this context it is worth mentioning that in the scenario which seems to fit the data best, with primary jet energy near $`10^{21}`$ eV and a galactic source overdensity of about $`10^5`$ (see Fig. 1 and ref.frag2 ), including the photon flux fully would only reduce the predicted event rate by a factor of two to three, since in this case the flux of $`10^{20}`$ eV photons at source is only slightly larger than the corresponding proton flux. This would still give a neutrino flux in easy striking range of km<sup>2</sup> scale detectors.

This paper shows that the neutrino flux accompanying the highest energy cosmic rays in top-down scenarios is of order of the limits placed by operating experiments such as AMANDA II, RICE and AGASA. Further data from these experiments, or next generation experiments IceCube and AUGER, can test the viability of top-down scenarios which generate the highest energy cosmic rays. If a signal is found soon, future high statistics experiments should be able to map out the neutrino spectrum, thereby allowing us direct experimental access to physics at energy scales many orders of magnitude beyond the scope of any conceivable particle collider on Earth.

###### Acknowledgements.

We would like to thank Michael Kachelriess, John Learned and Tom Weiler for valuable discussions. This work was supported in part by a DOE grant No. DE-FG02-95ER40896 and in part by the Wisconsin Alumni Research Foundation. The work of M.D. was partially supported by the SFB375 of the Deutsche Forschungsgemeinschaft.


format: text

text:
# Stirring Bose-Einstein condensate

## Abstract

By shining a tightly focused laser light on the condensate and moving the center of the beam along the spiral line one may stir the condensate and create vortices. It is shown that one can induce rotation of the condensate in the direction opposite to the direction of the stirring.

When a spoon (stick, or another object of a similar shape) is used to stir a liquid, the latter rotates in the direction induced by the stirring object. Is it possible to make the liquid rotating counter-clockwise while stirring it clockwise? The aim of this paper is to show that this counter-intuitive scenario may be realized in a quantum fluid, or more precisely a Bose-Einstein condensate (BEC) when stirring it with the help of a tightly focused laser beam.

Creation of vortices in a BEC and study of their properties has been a subject of quite intensive research for last couple of years (an extensive list of references may be found in ). Let us mention that vortices have been created in the BEC experimentally using various methods. Paris group  used the rotating anisotropic potential (created by a detuned broad laser beam) to make a direct analog of the rotating bucket experiments . The formation of a vortex is then the result of dynamical instabilities appearing in the course of the experiment . Similar method was used by Ketterle group . Boulder group  created vortices in two component condensate, where one fraction was made to rotate with respect to the other by means of the phase engineering technique. The latter technique attempts to create directly the desired vortex state in the condensate.

Various “stirring” propositions have been discussed theoretically  for creation of vortices. In particular  used a localized potential moving on a circular path (with an appropriate smooth turn-on and turn-off of the stirrer). Such a stirring produces the condensate state which may be approximately described as a time dependent combination of the ground state and the vortex state. As time evolves the system undergoes a generalized Rabi oscillation between the ground state and the vortex state.

Similar in spirit is our recent proposition for creation of vortices in a BEC . It relies on an appropriate deformation of a harmonic trapping potential by means of an additional, tightly focused laser beam. The beam approaches the center of the trap moving along a spiral line. The effective interaction of the detuned laser beam with atoms results in an additional effective potential seen by the atomic external degrees of freedom. Neglecting the interaction between atoms, the effective two-dimensional Hamiltonian, in the frame rotating with the center of the laser beam, reads

$`\widehat{H}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2}{x^2}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2}{y^2}}\mathrm{\Omega }\widehat{L}_z+{\displaystyle \frac{x^2+y^2}{2}}`$ (2)
$`u_0\mathrm{arctan}(|x_0|)\mathrm{exp}\left({\displaystyle \frac{(xx_0)^2+y^2}{2\sigma ^2}}\right),`$
where $`\mathrm{\Omega }`$ is the frequency of the rotation of the laser beam around the center of the trap while $`u_0>0`$, $`\sigma `$ and $`x_0`$ stand for the parameters of the beam. In Eq. (2) and in the following we work in units defined by the harmonic trap. Changing $`x_0`$ from an initial negative value to zero, according to $`x_0(t)=x_0(0)+vt`$, corresponds to the motion of the laser beam along a spiral line in the laboratory frame. When the beam reaches the center of the trap its intensity is reduced to zero \[see $`\mathrm{arctan}(|x_0|)`$ function in Eq. (2)\] and we end up with the harmonic trapping potential only (for details of the method see ).

We have shown, on the other hand,  that sweeping the laser beam across the condensate along the spiral line may serve as an efficient and stable way to create vortices in the system. This can be easily understood by looking at the energy levels of the Hamiltonian (2) for different (fixed) values of $`x_0`$, see Fig. 1. Narrow avoided crossings between neighboring energy levels indicate that starting with the system in the ground state and changing $`x_0`$ from some negative value to zero, one may pass the avoided crossing diabatically and end up (with a high efficiency) in the first excited state of the trap that possesses the angular momentum $`L_z=1`$. Moreover, the process may be repeated, Fig. 1 suggests that having the system in the first excited state of the harmonic trap with $`L_z=1`$ (after a single sweep) it is possible to employ the second similar process and transfer the population to $`L_z=2`$ state with a high efficiency.

To look more quantitatively at the stirring process, we write the wave function of the system, in the hydrodynamical approach , as $`\mathrm{\Psi }(\stackrel{}{r},t)=\sqrt{\rho (\stackrel{}{r},t)}\mathrm{exp}(i\chi (\stackrel{}{r},t))`$, where $`\rho (\stackrel{}{r},t)`$ is the density of a probability fluid. The velocity field is defined as

$$\stackrel{}{v}=\stackrel{}{}\chi (\stackrel{}{r},t).$$
(3)
The single valuedness of the wave function requires that the circulation of the velocity field $`\mathrm{\Gamma }_C`$ around any closed contour $`C`$ is quantized (Feynman-Onsager quantization condition )

$$\mathrm{\Gamma }_C=_C\stackrel{}{v}𝑑\stackrel{}{l}=2\pi n,$$
(4)
where $`n=0,\pm 1,\pm 2,\mathrm{}`$. The value of $`n`$ characterizes vortices in the wave function. We say that we have vortex with unit charge at a given point, when calculation of $`\mathrm{\Gamma }_C`$ gives $`n=1`$ as contour $`C`$, encircling that point, shrinks down to this point.

Creation of vortices by our method, which is nothing but a smooth time-dependent modification of the potential, requires a sudden (due to the quantization) appearance of a non-zero circulation. This is necessarily accompanied by an appearance of a singularity in the velocity field. It is interesting to find out how this process occurs since we know that at the beginning of the laser sweep there is no circulation in the velocity field, at the end there is a vortex approximately at the center of the trap. Integrating the time-dependent Schrödinger equation (recall that we discuss noninteracting particles first) we have looked for the wave function modulus minimum and calculated the circulation around small contour encircling it. If $`n`$ is equal to $`1\pm 0.04`$ ($`1\pm 0.04`$) we assume that vortex (antivortex) with unit charge is located at such a minimum. As almost non-interacting condensates have been realized in laboratories already , it is perfectly legitimate to consider non-interacting particle case first.

Let us inspect the first sweep of the laser beam through the system initially in the ground state. We have found (compare Fig. 2) that the vortex moves in from the border of the trap (i.e. the range of the configuration space where we are able to control the velocity field numerically). Please note that we are able to observe the vortex after some time since the beginning of the simulation. Indeed, it crosses the trap border instead of being suddenly created at $`t=145`$—see Fig. 2. The position of the vortex follows quite closely the center of the focused laser beam. At the end of the excitation process the vortex lands very close to the trap center.

Similarly for a second laser sweep aiming at increasing $`L_z`$ to 2 an additional vortex with the topological charge $`n=1`$ comes from the border of the trap along a spiral line (similar to the line depicted in Fig. 2) and collides with the first vortex which, during the whole time evolution, is situated in the vicinity of the trap center. In the numerical implementation the final wave function consists mainly of the eigenstate with $`L_z=2`$ (the square overlap on this state is $`p_20.9997`$). However, there is also a slight contribution from the $`L_z=0`$ eigenstate ($`p_00.0003`$). A simple calculation immediately shows that instead of a single vortex with $`n=2`$ we get two separate vortices with $`n=1`$ in this case. This observation confirms that vortex with $`n=2`$ is unstable. The two vortices are situated symmetrically with respect to the trap center at a distance $`2(2p_0/p_2)^{1/4}`$. Plot of the phase of the final wave function in the vicinity of the trap center confirms such prediction, see Fig. 3.

Energy levels of the Hamiltonian (2) as a function of $`x_0`$ have been calculated in Fig. 1 for $`\mathrm{\Omega }=0.6`$. For $`x_0=0`$ the ground state corresponds to $`L_z=0`$, the first excited state corresponds to $`L_z=1`$ and the second one to $`L_z=2`$. However, the order can be different if we decrease $`\mathrm{\Omega }`$. Indeed for $`\mathrm{\Omega }<1/3`$ the second excited state (for $`x_0=0`$) corresponds to $`L_z=1`$. It offers an opportunity for the following counterintuitive situation which is of main interest for our study. Suppose, we start with the $`L_z=0`$ state. After a potential sweeping we end up with a very high probability in the state with $`L_z=1`$ where the rotation of the probability fluid coincides with the rotation of the applied laser beam. Then another, identical stirring by our “laser spoon” results in probability fluid rotating in the opposite direction (a state with $`L_z=1`$)! Needless to say such a situation is quite surprising and no analogy to some process in a classical fluid appears.

The prediction based on Hamiltonian levels behavior can again be tested by a direct integration of the time-dependent Schrödinger equation and indeed the $`L_z=1`$ state is excited with very high accuracy. Analyzing the process of such change of the angular momentum from $`L_z=1`$ to $`L_z=1`$ by looking at the time dependent motion of vortices we find that the vortex with $`n=1`$ initially situated at the center moves out to the border of the trap while the other vortex (born at the border) with an opposite $`n=1`$ circulation arrives at the center along a complex trajectory shown in Fig. 4. The latter vortex, before reaching the center, experiences a sequence of collisions with another $`n=1`$ vortex that affects its trajectory. Therefore a transition from $`n=1`$ to $`n=1`$ case is a result of (a bit complicated as seen in Fig. 4) dynamics of vortices.

It remains to be seen whether the counterintuitive stirring scheme is feasible also in the presence of atom-atom interactions since so far we have presented a creation of vortices for a non-interacting BEC. It is known, however, that the stability of vortices may be strongly affected by the atom-atom interactions . To analyze the effect of interactions we have performed numerical integration of the Gross-Pitaevskii equation

$$i\frac{\mathrm{\Psi }}{t}=(\widehat{H}+g|\mathrm{\Psi }|^2)\mathrm{\Psi },$$
(5)
with $`\widehat{H}`$ given by (2). The interaction parameter $`g`$ is proportional to the number of atoms in the system and to the $`s`$-wave scattering length. In an experiment, $`g`$ can be easily of order of thousands but it can be also reduced to a very small value exploring Feshbach resonances . In the present work, we have chosen $`g=100`$ for the numerical calculations.

If the ground state of the system is chosen as an initial state, applying the potential sweeping allows one to obtain the $`L_z=1`$ state with a high efficiency as described in Ref. . We performed such numerical simulation taking $`\mathrm{\Omega }=0.1`$. Now, we apply the second similar laser sweep on the state obtained after the first one. It creates a vortex with the topological charge $`n=1`$ similarly as it takes place for a noninteracting BEC if $`\mathrm{\Omega }<1/3`$. However, contrary to the linear case, the initial vortex with $`n=1`$ does not disappear — the interaction between atoms makes the initial vortex more robust to the perturbation. The vortex with $`n=1`$ lands close to the center of the trap while the original $`n=1`$ moves to the edge of the trap. In effect the total angular momentum per particle is $`\widehat{L}_z=0.42`$ with the dispersion $`\sigma _L=\sqrt{\widehat{L}_z^2\widehat{L}_z^2}=1.13`$.

The position of vortices may be observed using the interference approach . In left frame of Fig. 5 the square modulus of the final wave function superimposed with a plane wave traveling vertically in the figure’s plane is presented. A vortex–antivortex pair, clearly visible in the figure, can be observed experimentally as the interference technique has been applied in a laboratory already . The appearance of such a vortex-antivortex pair might be interesting from an experimental point of view, since interactions between such pairs in BEC confined in a harmonic trap, are still an unexplored topic experimentally.

Increasing the frequency $`\mathrm{\Omega }`$ of the stirring to $`\mathrm{\Omega }=0.25`$ we come back to the case of two $`n=1`$ vortices discussed previously for noninteracting case. Please note, that now the energy spacing between eigenstates with different value $`L_z`$ has decreased, so it is possible to have a $`L_z=2`$ state as a second excited eigenstate for $`\mathrm{\Omega }<1/3`$. In the presence of atom-atom interaction ($`g=100`$) we again sweep the laser across the condensate twice, first stirring creates a single vortex, a second stirring process adds an additional vortex with the topological charge $`n=1`$. This is again in a qualitative agreement with the noninteracting case considered previously. Quantitatively, the final state is characterized by $`\widehat{L}_z=1.69`$ with $`\sigma _L=1.92`$. The interaction between atoms leads now to a much larger separation between the two vortices, see the right frame in Fig. 5. Indeed, the separation between them is now comparable with the size of the entire condensate (note different scales in Fig. 3 and Fig. 5).

It is interesting to ask what is the critical stirring frequency for a transition from the regime of ‘vortex-antivortex’ to that of ‘vortex-vortex’ production during the second laser sweep. We estimate the critical frequency $`\omega _c`$ as satisfying the following equation:

$`\mu (L_z=1)+\omega _c=\mu (L_z=2)2\omega _c`$

where $`\mu (L_z=1)`$ and $`\mu (L_z=2)`$ are chemical potentials of two lowest eigenstates of the time-independent GP equation. The latter are found solving the 2D equation $`(\widehat{H}+g|\mathrm{\Psi }|^2)\mathrm{\Psi }=\mu \mathrm{\Psi }`$, with Hamiltonian $`\widehat{H}=\frac{1}{2}\stackrel{}{}^2+\frac{1}{2}(x^2+y^2)`$, i.e., the Hamiltonian (2) in the laboratory frame without laser beam. For stirring frequencies $`\mathrm{\Omega }`$ lower than $`\omega _c`$ three lowest GP eigenstates, in the frame rotating with stirrer, possess angular momentum $`L_z=0,1,1`$ while in the case of frequencies higher (but not too high) than $`\omega _c`$ the order is $`L_z=0,1,2`$. The frequency $`\omega _c`$ is an upper bound for the real critical frequency since its definition is based solely on the ordering of eigenstates in the frame rotating with stirrer. Indeed, an efficient transfer requires also that the distance in energy between the level that we would like to populate and the next one should be sufficient to assure adiabaticity, which is by definition not the case when we stir the BEC with $`\mathrm{\Omega }=\omega _c`$. Therefore, one might expect that the optimal realization should require lower frequency, probably in the middle between $`\omega _c`$ and the lowest estimate for a creation of vortex-antivortex pairs (equal to $`0`$). That gives $`\omega _c/2`$ as a good guess. We expect that the critical frequency should be somewhere between these two estimates, namely between $`\omega _c/2`$ and $`\omega _c`$. Calculation for $`g=100`$ gives $`\omega _c=0.18`$ which interestingly compares with $`0.125\pm 0.01`$ determined from a direct integration of Eq. (5) for $`g=100`$, $`20u_025`$ and a duration of a single laser sweep between $`40`$ and $`60`$ (compare Fig. 5 and Eq. (2)). Therefore, a numerical calculation gives a value which is higher than $`\omega _c/2`$ and lower than $`\omega _c`$, even though the lower bound is just a rough estimate. Similar calculations of upper bounds for the critical frequency ($`\omega _c`$) yield $`0.237`$ for $`g=30`$ and $`0.195`$ for $`g=70`$.

Finally, we would like to comment on an influence of the stirring scheme’s details on final results in the interacting case. First of all, we have observed that the width $`\sigma `$ \[compare (2)\] should be small, of the order of $`0.2`$; two times bigger widths lead to a significant decrease in the stirring process’ efficiency. Secondly, the parameter $`u_0`$ has to be high enough, of the order of $`20`$, for an efficient transfer of atoms from the ground state to the vortex state(s). These two conditions provide non-trivial restrictions on laser beam width and intensity, respectively. Thirdly, changes of the switching time within about $`20\%`$ of a given time scale ($`8`$ periods of harmonic trap for $`g=100`$) do not affect the dynamics qualitatively. Further details can be found in .

To summarize we have investigated the details of the creation of vortices in BEC when the laser sweep scheme  is applied. Especially, we have shown that rotating the probability fluid by means of the “laser spoon” may introduce a circulation with the opposite direction with respect to the steering one.

We are grateful to an anonymous referee for bringing references  to our attention and several useful hints. Support of KBN under projects 5 P03B 088 21 (K.S. and J.Z.), 2 P03B 124 22 (B.D.) is acknowledged. B. D. is grateful for hospitality extended to him in Hannover during preparation of the final form of manuscript.


text:
# Evolution of the resistivity anisotropy in Bi2Sr2-xLaxCuO6+δ single crystals for a wide range of hole doping

## I INTRODUCTION

It has long been recognized that the peculiar $`c`$-axis transport properties Gray  of the high-$`T_c`$ cuprates are strong manifestations of the unusual electronic state in these materials. Anderson  In particular, the following two features have been considered to be most unusual: (i) The magnitude of the $`c`$-axis resistivity $`\rho _c`$ is orders of magnitude larger than that expected from band calculations, Pickett  leading to a huge resistivity anisotropy where $`\rho _c`$ is up to 105 times larger Martin  than the in-plane resistivity $`\rho _{ab}`$. (ii) The temperature dependence of $`\rho _c`$ is in most cases semiconducting or insulating (i.e., $`d\rho _c/dT<0`$), while that of $`\rho _{ab}`$ is metallic ($`d\rho _{ab}/dT>0`$); such a contrasting behavior Nakamura ; Takenaka ; Ando  is not expected in ordinary anisotropic metal. Gray ; Anderson

The feature (i) indicates that the strong correlations in the cuprate materials give rise to some unconventional mechanism which renormalizes the $`c`$-axis transfer-matrix element to a value much smaller than what is expected for an uncorrelated system. This renormalization of the $`c`$-axis transfer-matrix element is usually called “charge confinement”. There have been many theoretical proposals to explain the charge confinement in the cuprates, but it is fair to say that the confinement mechanism is not yet understood on the fundamental level. Since many of the theoretical models for the charge confinement, such as the resonating-valence-bond theories Zou ; Nagaosa  or self-organized stripe theories, Carlson ; Zaanen  are closely tied to the mechanism of the high-$`T_c`$ superconductivity, it remains to be important to study the $`c`$-axis transport in the cuprates.

The feature (ii) is now better understood than the feature (i), both experimentally and theoretically, though the understanding is still far from complete. To appreciate the current understanding of the feature (ii), it should first be recognized Ioffe  that the $`c`$-axis transport in the normal-state of the cuprates is essentially an incoherent tunneling process; this fact gives rise to the intrinsic Josephson effect Muller  as well as the Josephson-plasma resonance JPR  in the superconducting state. It should at the same time be noted that the tunneling conductance for ordinary metal-insulator-metal junctions should be independent of temperature, according to the Fermi golden rule; this means that the “semiconducting” behavior of $`\rho _c(T)`$ does not automatically arise from the tunneling nature of the $`c`$-axis transport. (In this sense, the expression “semiconducting behavior”, which implies that the transport is essentially an activation process over a gap, is inappropriate, though this expression is often naively used in this context.) Thus, there must be some unconventional mechanisms which cause the peculiar “insulating” temperature dependence of $`\rho _c`$ in the cuprates.

Experimentally, there seem to be two additive mechanisms that are both responsible for the “insulating” $`\rho _c(T)`$. The first is the charge-confinement mechanism which appears to become increasingly more effective with lowering temperature; this causes the tunneling matrix element to be reduced and keeps $`\rho _c`$ increasing with decreasing $`T`$. This point was most clearly demonstrated Ando  by the low-temperature normal-state resistivity measurement of Bi2Sr2-xLaxCuO6+δ (BSLCO) under 60 T, which showed that $`\rho _c`$ of Bi-2201 keeps increasing even down to 0.6 K. The more recently recognized mechanism is the effect of the pseudogap Timusk  which causes destruction of the Fermi surface starting from the ($`\pm \pi `$, 0), (0, $`\pm \pi `$) points, as has been demonstrated by angle-resolved photoemission spectroscopy (ARPES) measurements Norman ; this causes the available density of states (DOS) for tunneling to be reduced, and $`\rho _c`$ gets steeply increased as the pseudogap develops. Note that the $`c`$-axis matrix element of the cuprates has a particular $`𝐤`$-dependence Xiang  ($`𝐤`$ is the in-plane wave vector of the conduction electrons) that tends to amplify the contribution of the electrons on these gapped portion of the Fermi surface, and this $`𝐤`$-dependence of the matrix element causes $`\rho _c`$ to be very sensitive to the opening of the pseudogap. Ioffe  Recent tunneling-spectroscopic studies Suzuki ; Krasnov  of the intrinsic junctions in Bi2Sr2CaCu2O8+δ (Bi-2212) have made strong cases for the interpretation that the steep upturn in $`\rho _c(T)`$ below the pseudogap temperature $`T^{}`$ is largely due to this second mechanism. It is expected that these two independent mechanisms, confinement and pseudogap, together cause the rather complicated, “insulating” $`\rho _c(T)`$ behavior in the cuprates, though their relative roles in determining the actual temperature dependence is not clarified yet.

Given the above-mentioned recent understanding of the $`c`$-axis transport, it would be instructive to see how the two mechanisms (confinement and pseudogap) are actually working in a single system. For such a purpose, it is useful to elucidate the evolution of the behavior of $`\rho _c(T)`$, as well as the anisotropy ratio $`\rho _c/\rho _{ab}`$, over a wide range of hole doping, because $`T^{}`$ is known to decreases with increasing doping. In the past, there were many studies of the anisotropic resistivities in La2-xSrxCuO4 (LSCO), Nakamura ; Kimura ; logT ; Boebinger  YBa2Cu3O7-δ (YBCO), Takenaka ; Terasaki  Bi-2212, Watanabe97  and BSLCO, Martin ; Wang  but none of them cover a wide enough doping range from underdoped insulators to overdoped superconductors to yield a coherent picture.

Very recently, Komiya et al. measured Komiya  the anisotropic resistivity of high-quality LSCO single crystals for $`x`$ = 0.01 – 0.10 and found that in this system the effect of pseudogap is not apparent (i.e., the additional mechanism to cause the insulating $`\rho _c(T)`$ is absent), probably because of the unique Fermi surface topology Zhou ; InoFS  of LSCO. Therefore, it is useful to study another cuprate system in which one can change the hole doping over a really wide range. Obviously, BSLCO is an ideal candidate, because we have demonstrated Murayama ; Ono ; Hanaki  that the hole doping can be varied (by changing the La concentration) from heavily-underdoped insulator to overdoped superconductor; moreover, ARPES measurements have found Harris ; Chuang ; Sato  that the Fermi surface topology and the pseudogap structure of BSLCO are very similar to those of Bi-2212. We can thus expect that it may be possible to substantiate the effects of charge confinement and the pseudogap on the $`c`$-axis transport by tracing its evolution with the hole doping in BSLCO.

The purpose of this paper is to provide high-quality data, obtained from well-controlled experiments, on the evolution of the temperature dependences of $`\rho _c`$ and $`\rho _c/\rho _{ab}`$ of BSLCO for a wide range of doping, which showcases the complicated origin of the unusual $`c`$-axis transport in the cuprates. Our data exemplifies that there are three different mechanisms that govern the behavior of the anisotropy ratio $`\rho _c/\rho _{ab}`$: charge confinement, pseudogap, and localization. In particular, the temperature dependence of $`\rho _c`$ of our underdoped samples makes a good case for the two different origins of the insulating behavior below $`T^{}`$.

This paper is organized as follows: Details of our crystal-growth technique to produce the high-quality BSLCO single crystals, as well as the details of our measurement technique, are presented in Section II. Section III presents our data for $`\rho _{ab}(T)`$, $`\rho _c(T)`$, and $`\rho _c/\rho _{ab}`$, putting emphasis on how they evolve with changing hole doping. The implication of the systematics of our observation is discussed in Section IV, and Section V summarizes our findings.

## II EXPERIMENTS

### II.1 Bi2Sr2-xLaxCuO6+δ Crystals

The Bi2Sr2-xLaxCuO6+δ system is a Bi-based single-layer cuprate where the hole doping is controlled by replacing Sr with La. Until a few years ago, this material had been considered to be a particularly “dirty” system among the cuprates, because the in-plane resistivity of even the best crystal at those times showed the residual resistivities of around 70 $`\mu \mathrm{\Omega }`$cm (Refs. Martin, ; Ando, ) and the Hall coefficient $`R_H`$ displayed only a weak temperature dependence, Mackenzie  both of which are characteristic of dirty cuprate samples.

In 1999, Ando and Murayama demonstrated Murayama  that it is possible to grow high-quality single crystals of BSLCO that show much smaller residual resistivity and a strongly temperature-dependent $`R_H`$, which conform to the “standard” behavior of cuprates. Also, though the BSLCO system was believed to possess the lowest optimum $`T_c`$ among the well-studied cuprate systems, we have demonstrated Ono  that the optimum $`T_c`$ of this system can be raised to 38 K, which is almost equal to the optimum $`T_c`$ of LSCO.

In BSLCO, larger La concentration $`x`$ corresponds to smaller hole doping per Cu, $`p`$, and the relation between $`x`$ and $`p`$, as well as the phase diagram of $`T_c`$ vs $`p`$, has been sorted out Hanaki ; it was found that optimum doping occurs at $`p`$ 0.16 (which is achieved with $`x`$ 0.4), but the superconductivity disappears at $`p`$ 0.10 ($`x`$ 0.8) upon underdoping. Low-temperature normal-state resistivity measurement of this system using 60-T pulsed magnetic field found Ono  that the metal-to-insulator crossover occurs at $`p`$ 1/8, below which $`\rho _{ab}(T)`$ shows a peculiar $`\mathrm{log}(1/T)`$ divergence.

Single crystals of Bi2Sr2-xLaxCuO6+δ are grown by the floating-zone (FZ) technique for a wide range of La concentration (0.23 $`x`$ 1.02). Before the FZ operation, we first prepare polycrystalline rods of BSLCO. Raw powders of Bi2O3, SrCO3, La2O3, and CuO, with purities of 99.9% or higher, are dried, weighed, mixed into the nominal molar ratio of the target composition, and well ground in an agate mortar; they are then calcined at 750 – 850C for 20 hours to form the BSLCO phase. The resulting powders are reground and calcined again, and this process is repeated twice. The x-ray diffraction analysis reveals that the powders become 100%-pure BSLCO phase after the third calcination. The resulting BSLCO powders are isostatically pressed into a rod shape ($``$ 5 mm$`\varphi `$ $`\times `$ 100 mm) and finally sintered at 850C for 20 hours, to form a rigid polycrystalline feed rod to be used for the FZ operation.

Single-crystal growth is carried out using an infrared image furnace (NEC Machinery SC K-15HD) with two halogen lamps and double ellipsoidal mirrors. We grow BSLCO crystals with a “simple” FZ technique, which means that we do not use any solvent to reduce the temperature of the molten zone. One important technique we use is what we call “fast scan”, with which we quickly melt and quench the whole feed rod using the FZ furnace with the feed speed of 30 mm/h. This process is necessary for making the feed rod dense, which is important for avoiding the emergence of bubbles during the crystal growth; actually, the rod becomes appreciably thinner than the as-sintered rod after the fast scan (it becomes roughly 3 mm$`\varphi `$). The “scanned” rod is then remounted to the FZ furnace and is used for the actual crystal growth. The growth rate is kept constant at 0.5 mm/h, which is two times slower than the growth rate for LSCO. Komiya  Another important technique we use during the growth is what we call “necking”; since many independent domains start to crystallize at the initial stage of the growth, after roughly 1 cm of the rod is grown, the diameter of the molten zone is gradually reduced to approximately 2 mm$`\varphi `$, which causes fewer domains to remain in the crystal. The molten zone is then gradually fattened, which leads to the growth of the fewer domains in the resulting rod. After a few times of the necking, when the growth is successful, there remains only one domain and the whole chunk becomes a large single-domain crystal. The space around the molten zone is filled with 1 atm of flowing oxygen during the growth of underdoped crystals, while for the growth of overdoped crystals the space is filled with 1 atm of flowing air. The size of the single-domain crystals that can be obtained from the grown rod is typically 20 $`\times `$ 3 mm2.

The actual La concentrations in the crystals are determined by the inductively-coupled plasma atomic-emission spectroscopy (ICP-AES). The crystals tend to contain more La than the nominal composition, but the difference is small; here we report crystals with actual $`x`$ values of 0.23, 0.39, 0.49, 0.66, 0.73, 0.84, 0.92, and 1.02, and these crystals were obtained from the nominal $`x`$ values of 0.20, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, and 1.00, respectively. In all these crystals, the ICP-AES analysis confirms that there is essentially no excess Bi in the crystals; rather, Bi tends to be deficient in highly La-doped crystals, as shown in Table I. On the other hand, the ratio of (Sr+La):Cu is always near 2:1.

Uniformity of the distribution of the La dopants in the crystals is confirmed by the electron-probe microanalysis (EPMA), with which the variation of the $`x`$ value is estimated to be within 5%. The superconducting transition temperature $`T_c`$ is determined both by the zero resistivity in the transport measurement and by the onset of the Meissner signal in the SQUID magnetization measurement; for all the crystals reported here, these two definitions give $`T_c`$ values that match within $`\pm `$0.5 K, and the variation among different pieces of crystals with the same composition is also within $`\pm `$0.5 K. (Thus we refer to the $`T_c`$ values with two significant digits in this paper.) The $`T_c`$ values for each of the compositions are listed in Table II.

### II.2 Measurements

The crystals are cut into dimensions of typically 2 $`\times `$ 1 $`\times `$ 0.05 mm3 for the $`\rho _{ab}`$ measurements and 1 $`\times `$ 1 $`\times `$ 0.05 mm3 for the $`\rho _c`$ measurements. Accurate determination of the crystal thickness is done by measuring the weight with 0.1 $`\mu `$g resolution and converting it into volume using the nominal density for each composition. We employ the straightforward four-terminal method to measure $`\rho _{ab}`$ and $`\rho _c`$: For the $`\rho _{ab}`$ measurements, the current contacts are painted to cover two opposing side faces of the platelet-shaped crystals to ensure a uniform current flow and the voltage contacts are painted on the two remaining side faces of the crystals. For the $`\rho _c`$ measurements, two current contacts are painted to almost completely cover the opposing $`ab`$ faces of the crystal and two voltage contacts are placed in the small window reserved in the center of the current contacts. In both cases, the contact pads are tactfully hand-drawn with gold paint, followed by a heat treatment at 400C for 30 minutes in the air, which makes the gold particles to well adhere to the sample surface. After this heat-treatment to cure the gold contact pads, the samples are annealed at higher temperatures for longer time to control the oxygen concentration; samples with $`x`$ larger than 0.3 are annealed in flowing air at 650C for 48 hours, while those with $`x0.3`$ are annealed in flowing oxygen at 400C for 60 hours. In both cases, the crystals are annealed together with a sufficient amount of BSLCO powders of the same composition, and they are quenched to room temperature at the end of the annealing. Finally, gold lead wires are attached to the contact pads using silver epoxy, which is cured at relatively low temperature, 130C. The contact resistance we achieve with this technique is less than 1 $`\mathrm{\Omega }`$. We note that the above annealing conditions are chosen to minimize the superconducting transition width $`\mathrm{\Delta }T_c`$ of our crystals; in all the superconducting samples, $`\mathrm{\Delta }T_c`$ in $`\rho _c(T)`$ data is around 2 K after the above annealing.

Before finishing the description of our sample preparation technique, we would like to emphasize the importance of the quenching at the end of the high-temperature annealing. Figure 1 shows the data of $`\rho _{ab}(T)`$ and $`\rho _c(T)`$ of our BSLCO crystals with $`x`$ = 0.66, for both quenched and slowly-cooled samples. Although both sets of the samples are annealed at 650C for 48 hours, the data show non-negligible differences; the difference in the $`\rho _{ab}(T)`$ are comparatively small (the magnitude of $`\rho _{ab}`$ of quenched sample is 10% larger than that of slowly-cooled sample), but the difference in $`\rho _c(T)`$ is significant, differing by about 50%. Also, $`\mathrm{\Delta }T_c`$ is about 2 K for the quenched samples, while it is about 6 K for the slowly-cooled sample. Those differences are consistent with the interpretation that the slow cooling causes the oxygen concentration in the samples to be inhomogeneous and there are more oxygen near the surface of the slowly-cooled samples; remember that the equilibrium oxygen concentration of Bi-2201 during the annealing is dependent on both the temperature and the oxygen partial pressure, and the equilibrium concentration becomes higher as the temperature is lowered. Thus, to obtain uniformly oxygenated samples and reproducible data, quenching is indispensable.

In this work, we measure $`\rho _{ab}`$ and $`\rho _c`$ on more than 20 samples for each La concentration to check for the reproducibility; the variation in the absolute magnitude of the resistivity of the quenched sample is less than $`\pm `$10%, which is consistent with the size of the errors due to uncertainties in the geometrical factors.

We have previously determined the actual hole concentration per Cu, $`p`$, for various La concentrations $`x`$, and have obtained an empirical relation. Hanaki  In this paper we use this empirical relation to convert $`x`$ into $`p`$, which are listed in Table I. Hereafter we mostly use the $`p`$ value in referring to the samples, because $`p`$ is much more intuitively understandable than $`x`$.

## III RESULTS

### III.1 In-plane resistivity

Figure 2 shows the temperature dependences of $`\rho _{ab}`$ for all the $`p`$ values studied here. As expected (and has been reported in Refs. Murayama, ; Ono, ; Hanaki, ), the magnitude of $`\rho _{ab}(T)`$ shows a systematic increase with decreasing carrier concentration from $`p`$ = 0.18 to 0.03. The $`\rho _{ab}(T)`$ curve of the overdoped sample ($`p`$ = 0.18) shows a behavior that can be described Murayama  by $`a+bT^n`$ with $`n>1`$, while a strictly $`T`$-linear behavior is apparent in the optimally-doped sample ($`p`$ = 0.16). In the most underdoped non-superconducting sample ($`p`$ = 0.03), $`\rho _{ab}(T)`$ shows an insulating behavior ($`d\rho _{ab}/dT<0`$) from room temperature, which is quite different from the $`\rho _{ab}(T)`$ behavior of lightly hole-doped LSCO that shows a metallic behavior ($`d\rho _{ab}/dT>0`$) even at $`x`$ = 0.01. mobility  We assert this difference is not intrinsic but is rather due to disorder, because the room-temperature sheet resistance per CuO2 plane of our BSLCO at $`p`$ = 0.03 is 110 k$`\mathrm{\Omega }`$, while that of LSCO at $`x`$ = 0.03 (Ref. mobility, ) is only 60 k$`\mathrm{\Omega }`$; this indicates that a significant disorder is present in the heavily-underdoped BSLCO samples. On the other hand, the residual resistivity at optimum-doping ($`p`$ = 0.16) is only 20 $`\mu \mathrm{\Omega }`$cm, indicating that the electron transport can be quite clean in this system; since we find that the crystals are morphologically very clean for both $`p`$ = 0.03 and 0.16, the source of the electronic disorder is apparently not any crystalline defects but is likely to be the local disordering potentials due to the La dopants.

It is useful to note that the barely superconducting sample, $`p`$ = 0.10, shows $`\rho _{ab}`$ value of $``$0.7 m$`\mathrm{\Omega }`$cm just above $`T_c`$; this value corresponds to the sheet resistance per CuO2 plane of 5.8 k$`\mathrm{\Omega }`$, which is very close to the quantum sheet resistance $`h/4e^2`$ (= 6.45 k$`\mathrm{\Omega }`$). This suggests that the disappearance of the superconductivity in BSLCO at unusually high doping of $`p0.10`$ is caused by the disorder-driven superconductor-insulator transition in two-dimensional (2D) superconductors. Fukuzumi

### III.2 Out-of-plane resistivity

Figure 3 shows the temperature dependences of $`\rho _c`$ for the same compositions as those in Fig. 2. As is the case with $`\rho _{ab}(T)`$, the magnitude of $`\rho _c(T)`$ shows a systematic increase with decreasing carrier concentration. Since the magnitude of $`\rho _c`$ changes by 2 orders of magnitude from $`p`$ = 0.18 to 0.03, the data for the overdoped and optimally-doped samples are shown in Fig. 3(a), while those of the underdoped samples are shown in Fig. 3(b), both in linear scale; Fig. 3(c) plots $`\rho _c(T)`$ in semi-logarithmic scale to equally show the behavior of all the samples. In Fig. 3(a), $`\rho _c(T)`$ data for La-free Bi2.13Sr1.89CuO6+δ (“pure” sample) are also shown for comparison (dotted line). It is worthwhile to note that, for each composition, the $`\rho _c`$ value of our crystals is generally smaller than that reported for BSLCO in the past. Martin ; Wang

We note that a “metallic” temperature dependence of $`\rho _c`$ is often observed, as is the case with the present experiment, in overdoped samples at high enough temperatures, Nakamura ; Takenaka ; Watanabe ; Lavrov  which probably indicates that the in-plane scattering of the electrons is involved in determining the temperature dependence of $`\rho _c`$. Jayannavar  In fact, recent $`c`$-axis magnetoresistance study by Hussey et al. showed that $`\rho _c`$ of LSCO is strongly affected by the scattering within the planes. Hussey  This suggests that looking at $`\rho _c/\rho _{ab}`$ is useful for investigating the intrinsic mechanisms of the unusual $`c`$-axis transport, because the in-plane scattering rate is roughly cancelled out in $`\rho _c/\rho _{ab}`$.

### III.3 Anisotropy ratio $`\rho _c`$/$`\rho _{ab}`$

Figure 4 shows the temperature dependences of the anisotropy ratio $`\rho _c`$/$`\rho _{ab}`$ for all the concentrations studied. At moderate temperatures, the magnitude of $`\rho _c`$/$`\rho _{ab}`$ monotonically increases with decreasing $`p`$, except for the most underdoped concentration ($`p`$ = 0.03) at which $`\rho _c`$/$`\rho _{ab}`$ suddenly drops; it is probably the case that this drop is caused by the heightened $`\rho _{ab}`$ value of the $`p`$ = 0.03 sample due to localization, and therefore if there were no enhanced disorder at this composition we would expect $`\rho _c`$/$`\rho _{ab}`$ to be not much different from that for $`p`$ = 0.07. Note that the $`\rho _c`$/$`\rho _{ab}`$ value generally increases with lowering temperature, and almost reaches 106 just above $`T_c`$ for $`p=0.100.14`$.

Before our FZ-grown crystals became available, the anisotropic resistivity of BSLCO was measured by Wang et al. using flux-grown crystals that showed rather large residual resistivities. Wang  Though the data by Wang et al. demonstrated reasonably systematic behavior and are mostly consistent with our data reported here, their data are different from ours in one important aspect: they observed that the anisotropy ratio $`\rho _c/\rho _{ab}`$ is almost independent of doping at 300 K, while we observed that $`\rho _c/\rho _{ab}`$ at 300 K increases systematically with decreasing $`p`$, as can be seen in Fig. 4.

We note that in the anisotropic resistivity measurements using the flux-grown crystals, it was already recognized that the BSLCO system shows a very large anisotropy ratio $`\rho _c/\rho _{ab}`$ of more than 105, which is the largest among the cuprates; for example, Martin et al. reported Martin  that $`\rho _c/\rho _{ab}`$ can be as large as 3$`\times `$105 in Bi-2201, and the work by Wang et al. mentioned above reported Wang  maximum $`\rho _c/\rho _{ab}`$ of 2.5$`\times `$105. In our crystals, perhaps because of the smaller $`\rho _{ab}`$ values, the observed $`\rho _c/\rho _{ab}`$ is even larger: the $`\rho _c/\rho _{ab}`$ value exceeds 5$`\times `$105 in four of the concentrations ($`p`$ = 0.10, 0.11, 0.12, and 0.14), and the maximum $`\rho _c/\rho _{ab}`$ we observe is 8$`\times `$105.

## IV DISCUSSIONS

### IV.1 La-free (“pure”) sample

In Fig. 3, one can infer that the onset of steeply insulating behavior in $`\rho _c(T)`$ moves to higher temperature with decreasing $`p`$ in the La-doped series of samples, which is consistent with the expected behavior of the pseudogap-opening temperature $`T^{}`$. Although the data for pure sample in Fig. 3(a) appears to conform nicely to this trend, there is one issue that we need to elaborated on: When we estimate the actual hole concentration of the pure sample using the method we reported in Ref. Hanaki, , the $`p`$ value of this pure sample is 0.17 \[Fig. 5(a)\]; the slope of $`\rho _{ab}(T)`$ also indicate that the $`p`$ value of the pure samples should be around 0.17 \[Fig. 5(b)\]. On the other hand, compared to the La-doped $`p`$ = 0.18 sample, the $`\rho _c`$ value of the pure sample is smaller and the onset temperature of the insulating behavior is lower, both of which would normally be expected for more overdoped sample. Although this apparent puzzle looks problematic, it can be understood to be a result of a combination of extrinsic effects: First, since the constituent atoms of the apical-oxygen layer are different \[(Sr,La)O for BSLCO and (Sr,Bi)O for pure samples\], it is expected that the $`c`$-axis tunneling matrix element $`t_{}`$ is also different; it appears that $`t_{}`$ is larger in the pure samples, though the microscopic reason for the difference is not very clear. Second, since the $`c`$-axis magnetoresistance study of BSLCO has revealed Lavrov  that in the overdoped region the pseudogap is largely due to superconducting fluctuations and thus $`T^{}`$ is governed primarily by $`T_c`$, the low $`T_c`$ value of the pure sample is the reason for the low onset temperature of the insulating behavior. note

In Fig. 5(b), one notices that residual resistivity of the pure sample is higher than that of the La-doped samples, which may be counter-intuitive for a “pure” sample. It is important to recognize, however, that the composition of the “pure” sample is Bi2.13Sr1.89CuO6+δ, which means that there is sizable disorder in the Sr-O block due to excess Bi. (It is known that stoichiometric Bi2.0Sr2.0CuO6+δ is unstable and cannot be grown.) The larger residual resistivity and the lower $`T_c`$ of the pure sample suggest that the disorder due to (Sr,Bi) substitution is more harmful to the electronic system in this material than the disorder due to (Sr,La) substitution.

### IV.2 Two origins of the pseudogap

It is worthwhile to mention that previous measurements of the $`c`$-axis magnetoresistance of BSLCO have uncovered two additive mechanisms for the pseudogap.Lavrov  In underdoped samples, a magnetic-field-insensitive pseudogap emerges at $`T^{}`$, and a distinct one is formed above $`T_c`$ as a precursor to superconductivity (which is essentially the superconducting fluctuations); it was found that the latter is mostly responsible for large negative magnetoresistance in $`\rho _c`$. In overdoped samples, on the other hand, it is only the precursor to superconductivity that appears above $`T_c`$ and serves as pseudogap.

Recently, scanning tunneling spectroscopy (STS) measurements have been done on our pure samples (which is overdoped) and the pseudogap was found to open below $`68\pm 2`$ K. Kugler  In Fig. 3(a), the minimum in $`\rho _c(T)`$ of the pure sample is observed at 70 K, which actually agrees with $`T^{}`$ determined by STS and gives strong support to the interpretation that the steeply insulating behavior in $`\rho _c(T)`$ is largely due to the suppression of the DOS in BSLCO. Moreover, the STS study has concluded Kugler  that the pseudogap and the superconductivity has a common origin in our pure sample, which is in agreement with the conclusion drawn from the magnetoresistance study Lavrov  that the pseudogap in the overdoped region is in fact a precursor to superconductivity.

### IV.3 Temperature dependence of $`\rho _c`$

During the course of the analysis of our $`\rho _c(T)`$ data, we found that the log-log plot of the data (Fig. 6) uncovers intriguing systematics. In Fig. 6, one can see that for most of the compositions there is a finite range of temperature where the $`\rho _c(T)`$ data can be roughly described by $`T^\alpha `$ with $`\alpha 1`$ (thin solid lines in Fig. 6). Since the pseudogap is known to cause a rapid increase in $`\rho _c`$ with decreasing $`T`$, it is natural to interpret that the observed $`T^\alpha `$ dependence of $`\rho _c`$ is caused by the development of the pseudogap; this interpretation motivates us to evaluate the temperature $`T_1`$ below which the $`\rho _c(T)`$ obeys the $`T^\alpha `$ dependence. As shown in Fig. 7, the $`p`$-dependence of $`T_1`$, as well as the value of $`T_1`$ itself, is consistent with what one expects for the pseudogap temperature $`T^{}`$.Takenaka ; Timusk ; Watanabe ; Lavrov

An interesting feature in the data in Fig. 6 is that the $`\rho _c(T)`$ data tend to deviate downwardly from the $`T^\alpha `$ ($`\alpha `$ 1) dependence at lower temperatures in underdoped samples; moreover, the characteristic temperature for this deviation, $`T_2`$, systematically moves to higher temperature with decreasing $`p`$. In the light of the current understanding of the $`c`$-axis transport phenomena in cuprates, it is most reasonable to interpret this deviation to come from a saturation of the rapid pseudogap opening, and the more slower increase in $`\rho _c`$ below $`T_2`$ is likely to be due to the temperature dependence in the strength of the confinement. Therefore, the temperature dependence of $`\rho _c`$ as revealed in the log-log plot appears to reflect the fact that the “insulating” behavior is caused by the two phenomena, the pseudogap and the confinement. This is new information which the systematics of our data allows us to infer. (The low-temperature data of the non-superconducting samples ($`p`$ = 0.03 and 0.07), which show steeper divergence below $`T_2`$ compared to the superconducting samples, are probably affected by the additional effect of localization.)

We note that the effect of superconducting fluctuations, which cause additional pseudo-gapping near $`T_c`$,Lavrov  appears to be inconsequential in the overall behavior of the underdoped samples shown in Fig. 6, and the evolution of $`T_2`$ with $`p`$ is unrelated to $`T_c`$. This is probably because the additional reduction in the DOS near $`T_c`$ due to superconducting fluctuations is relatively small in underdoped samples, which can be inferred by the small size of the magnetoresistanceLavrov  observed in underdoped samples.

It is useful to mention that several empirical temperature dependences have been proposed to describe the $`\rho _c(T)`$ data of cuprates. In particular, $`\rho _c(T)=(a/T)e^{(\mathrm{\Delta }/T)}+bT+c`$ has been reported to fit the $`\rho _c(T)`$ data of Bi-2212 very well, Watanabe97 ; Yan  while $`\rho _c(T)=a_0+aT+1/(bT+c)`$ has been reported to fit the $`\rho _c(T)`$ data of BSLCO. Wang  We have tried to fit our data using these formulae, and found that the latter formula fits well to the data for $`p`$ = 0.16 and 0.18; however, neither of them fits well to the data for underdoped samples, which is apparently due to the existence of “$`T_2`$” below which the divergence of $`\rho _c(T)`$ weakens. Therefore, the extended range of doping of our study unveils the limitation of empirical temperature dependences.

### IV.4 $`\rho _c`$/$`\rho _{ab}`$

Figure 8 replots the $`\rho _c/\rho _{ab}`$ data of Fig. 4 using a linear scale, which makes it easier to infer the quantitative change of the anisotropy with doping. In Fig. 8, one can see that at high temperatures the $`\rho _c/\rho _{ab}`$ value drops rapidly with increasing $`p`$ for $`p0.07`$, but this rapid change tends to saturate in highly doped samples ($`p0.14`$); this situation can be understood more graphically by plotting the $`p`$ dependence of $`\rho _c`$/$`\rho _{ab}`$ at 200 K and 300 K, as is shown in Fig. 9. It can be seen in Fig. 9 that at 200 K the $`p`$ dependence of $`\rho _c/\rho _{ab}`$ for high doping ( $`p0.14`$) is weak, while the $`\rho _c/\rho _{ab}`$ value steeply increases with decreasing $`p`$ for lower doping ($`p0.12`$), except for $`p`$ = 0.03. It is most likely that the weak $`p`$ dependence at high doping reflects the change in the strength of the charge confinement, while the steep increase in the underdoped samples reflects the enhancement of the anisotropy due to the pseudogap. In fact, the $`T_1`$ value we determined in Fig. 6 is near 200 K for $`p`$ = 0.12. Therefore, the peculiar hole-doping dependence of the anisotropy is understood to be a combined result of the confinement and the pseudogap.

In Fig. 8, one notices that for non-superconducting concentrations ($`p`$ = 0.03 and 0.07) there is a peak in the temperature dependence of $`\rho _c/\rho _{ab}`$ (at 50 and 100 K, respectively), and below the peak temperature the anisotropy is diminished. In particular, $`\rho _c/\rho _{ab}`$ for $`p`$ = 0.03 decreases dramatically with decreasing $`T`$, suggesting that the system is heading towards an anisotropic three-dimensional (3D) state for $`T0`$. Such diminishment of $`\rho _c/\rho _{ab}`$ has also been observed in insulating samples of Bi-2212 (Ref. Kitajima, ) and LSCO (Ref. Komiya, ), and has been discussed to be due to localization of carriers. In fact, when the carriers are completely localized and becomes immobile in any direction, it is natural to see the system as a “3D” insulator; note that the system is expected to retain some amount of anisotropy in the 3D-insulating state, reflecting the anisotropy of the “bare” matrix elements.

Using the resistivity data for $`p`$ = 0.03, where the superconductivity is lost and carriers are localized due to disorder, we can actually confirm that the system is approaching a 3D insulator. Figure 10 shows the 3D variable range hopping (VRH) plots of $`\rho _{ab}(T)`$ and $`\rho _c(T)`$ for $`p`$ = 0.03. The dotted lines are fits to the expected temperature dependence for the 3D VRH, $`\rho \mathrm{exp}[(T_0/T)^{1/4}]`$. Since the 3D-VRH formula describes the temperature dependences of both $`\rho _{ab}`$ and $`\rho _c`$ very well, one can conclude that the charge transport along all directions is governed by the same mechanism, giving a rationale to call the system to be an anisotropic 3D insulator. The fits give $`T_0`$ values of $`2.6\times 10^6`$ K and $`8.6\times 10^4`$ K for $`\rho _{ab}(T)`$ and $`\rho _c(T)`$, respectively. These values are in the same range as the $`T_0`$ values obtained for flux-grown LSCO crystals in the lightly-doped regime, Chen  but the source of its anisotropy is not clear at this stage.

## V SUMMARY

Systematic measurements of the $`c`$-axis resistivity $`\rho _c`$ and the resistivity anisotropy ratio $`\rho _c/\rho _{ab}`$ of BSLCO for a wide doping range showcase that the anisotropy is determined by three mechanisms: charge confinement, pseudogap, and localization. At high temperature and in highly doped samples, the charge confinement is the sole determining factor, and thus one can measure the intrinsic strength of the charge confinement in this corner of the phase diagram; our data indicate that without the effect of the pseudogap the $`\rho _c/\rho _{ab}`$ value is order of 104 for BSLCO. At low temperature and in low doped samples, the localization governs the charge transport and the $`\rho _c/\rho _{ab}`$ ratio is diminished. In underdoped superconducting samples, the pseudogap significantly enhances the $`\rho _c/\rho _{ab}`$ value with decreasing temperature, causing $`\rho _c/\rho _{ab}`$ to reach $``$106. Close examination of the temperature dependence of $`\rho _c`$ of underdoped samples finds that below $`T^{}`$ there are two distinct temperature regions where the divergence of $`\rho _c`$ appears to be governed by different mechanisms; the data suggest that the development of the pseudogap leads to a $`T^\alpha `$ dependence of $`\rho _c`$ with $`\alpha `$ 1 below $`T^{}`$, while at lower temperature $`\rho _c`$ shows a weaker temperature dependence where we argue that the charge confinement dictates the $`T`$ dependence. Since the effect of the pseudogap is apparently lacking in the anisotropic behavior of LSCO, Komiya  the BSLCO system offers a complementary testing ground to study the source of the peculiar $`c`$-axis transport in the cuprate superconductors.

We thank A. N. Lavrov for helpful discussions.


text:
# Multiband model for penetration depth in MgB2

## Abstract

The results of first principles calculations of the electronic structure and the electron-phonon interaction in MgB2 are used to study theoretically the temperature dependence and anisotropy of the magnetic field penetration depth. The effects of impurity scattering are essential for a proper description of the experimental results. We compare our results with experimental data and we argue that the two-band model describes the data rather well.

The electronic structure of the recently discovered superconductorakimitsu  MgB2 is now rather well understood and the superconductivity may be ascribed to the conventional electron-phonon mechanismKortus ; An ; Kong ; Bohnen . The Fermi surface consists of two three-dimensional sheets, from the $`\pi `$ bonding and antibonding bands, and two nearly cylindrical sheets from the two-dimensional $`\sigma `$ bands. The qualitative difference between the 2-D $`\sigma `$\- and the 3-D $`\pi `$-bands in connection with the large disparity of the electron-phonon interaction (EPI) for the different Fermi surface sheets suggested a multiband description of superconductivity Liu ; Shulga ; Brink ; Golubov ; Choi . Recent reports on quantum oscillations Yelland  provided not only important information on the electronic structure near the Fermi level but it also probed directly the disparity of the EPI in the $`\pi `$\- and $`\sigma `$-band systems. The excellent agreement of the calculated EPI with the de Haas - van Alphen mass renormalization clearly confirms the basic assumption of the two-gap model Mazin-dhva . Further experimental support of this model comes from STM and point-contact spectroscopy Giub ; Bugoslavsky ; Szabo , high-resolution photo-emission spectroscopy Tsuda , Raman spectroscopy Chen , specific heat measurements sph ; Golubov  and studies of the magnetic penetration depth Nieder ; Klein ; Kim .

There is still some debate concerning the applicability of the multiband description to MgB2, in particular since some tunneling measurements Gonelli  show only a single gap with a magnitude smaller than the BCS value of $`\mathrm{\Delta }=1.76`$ $`T_c`$. A recently proposed multiband scenario for tunneling Brink  in MgB2 explains the reason for the differences in the observed tunneling spectra and thus helps to settle this debate.

A similar discussion has been going on concerning the penetration depth. The measured magnetic penetration depth shows a large variety of behavior (see Table 1). In order to interprete these results, a microscopic model is required. In this paper we shall use the multiband model Liu ; Brink ; Golubov  to calculate the temperature dependence and the anisotropy of the penetration depth using the Eliashberg formalism and the results of first principles electronic structure calculations.

Generalization of the BCS theory to the multiband model was first suggested in Refs. Suhl,  and Moskal,  and it has been observed experimentally in Nb doped SrTiO3 Bednorz . More recently, Kresin and Wolf Kre  suggested a two-band model in the strong-coupling regime for describing the properties of high Tc superconductors. Strong-coupling two-band-model calculations of the microwave response, and in particular the penetration depth, were performed in Refs. Adrian,  and Zhu, .

The penetration depth of the magnetic field $`\lambda _{L,\alpha \beta }`$ in the local (London) limit is related to the imaginary part of the optical conductivity by

$$1/\lambda _{L,\alpha \beta }^2=\underset{\omega 0}{lim}4\pi \omega Im\sigma ^{\alpha \beta }(\omega ,𝐪=0)/c^2,$$
(1)
where $`\alpha ,\beta `$ denote Cartesian coordinates and $`c`$ is the velocity of light. If we neglect strong-coupling effects (or, more generally, Fermi-liquid effects) then for a clean uniform superconductor at $`T=0`$ we have the relation $`\lambda _{L,\alpha \beta }=c/\omega _p^{\alpha \beta }`$, where $`(\omega _p^{\alpha \beta })^2=8\pi e^2_{𝐤j}\delta (\epsilon ^{𝐤j})v_F^\alpha v_F^\beta `$ is the squared plasma frequency and $`v_F^\alpha `$ the $`\alpha `$-component of the Fermi velocity. Impurities and interaction effects drastically enhance the penetration depth, and it is therefore suitable to introduce a so called ’superfluid plasma frequency’ $`\omega _{p,\alpha \beta }^{sf}`$ by the relation $`\omega _{p,\alpha \beta }^{sf}=c/\lambda _{L,\alpha \beta }`$. It has often been mentioned that this function corresponds to the charge density of the superfluid condensate, but we would like to point out that this is only the case for noninteracting clean systems at $`T=0`$.

In the two-band model we have the standard expression (neglecting vertex corrections)Adrian ; Zhu

$`1/\lambda _{L,\alpha \beta }^2(T)(\omega _{p,\alpha \beta }^{sf}(T)/c)^2=`$
$`{\displaystyle \underset{i=\sigma ,\pi }{}}\left({\displaystyle \frac{\omega _{p,i}^{\alpha \beta }}{c}}\right)^2\pi T{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}_i^2(n)}{[\stackrel{~}{\omega }_i^2(n)+\stackrel{~}{\mathrm{\Delta }}_i^2(n)]^{3/2}}},`$ (2)
where $`\stackrel{~}{\omega }(n)=\omega _nZ(\omega _n)`$ and $`\stackrel{~}{\mathrm{\Delta }}(\omega _n)=\mathrm{\Delta }(\omega _n)Z(\omega _n)`$ are the solutions of the Eliashberg equations Brink ; Golubov  and the calculated plasma frequencies for the $`\sigma `$ and $`\pi `$ bands are given in Ref. Brink, . The Eliashberg equations were solved numerically as described in Refs. Brink,  and Golubov, . The influence of impurities is incorporated into the model by including shifts of the gap function $`\mathrm{\Delta }_i^0(\omega _n)`$ and the renormalization factor $`Z_i^0(\omega _n)`$

$`\mathrm{\Delta }_i`$ $``$ $`\mathrm{\Delta }_i^0+{\displaystyle \underset{j}{}}\gamma _{ij}\mathrm{\Delta }_j^0/2\sqrt{\omega _n^2+(\mathrm{\Delta }_j^0)^2},`$
$`Z_i(\omega _n)`$ $``$ $`Z_i^0(\omega _n)+{\displaystyle \underset{j}{}}\gamma _{ij}/2\sqrt{\omega _n^2+(\mathrm{\Delta }_j^0)^2}`$
in the Eliashberg equations. Intraband scattering does not change $`T_c`$ and the gap values (Anderson’s theorem), but influences strongly the penetration depth.

Before we start discussing the exact solutions to the Eliashberg equations we present a simplified model consisting of two independent BCS superconducting bands with different plasma frequencies and different gaps (and consequently different $`T_c`$’s). In spite of the fact that this model is clearly an oversimplification, it captures qualitatively most of the observed behavior. In this model the band with the larger $`T_c`$ has the smaller plasma frequency (see Fig. 1).

For a clean system the resulting inverse squared penetration depth is the sum of the ’superfluid plasma frequencies’ (solid line). The kink in Fig. 1 is an artefact of the simplified model which will be smoothed out by interband scattering, but an inflection point in the temperature dependence of the penetration depth remains which has been observed experimentally. The low temperature dependence is determined by the band with the smallest gap, whereas the high temperature behavior results from the band with the larger gap. This is in accordance with the temperature dependence for low temperatures observed in some experiments. If the superconducting band with the smaller gap will be ’overdamped’ due to impurities, then the penetration depth is only determined by the other band and it will show a BCS temperature dependence, which has also been observed in some experiments.

For a proper understanding of the observed physical behavior of MgB2 it is important that impurities are taken into account properly. Recently, the influence of impurities on the two-gap superconductivity has been discussedmazin . Using the same arguments we shall discuss two cases: (a) The clean case with scattering rates $`\gamma _\sigma =\gamma _\pi `$=2meV as realized in low-resistivity dense wires wires  and (b) the dirty case with $`\gamma _\sigma `$=54meV and $`\gamma _\pi `$=1.2eV. The values for the scattering rates in the dirty case, as well as negligibly small interband scattering rates mazin , are in accordance with the results on high resistivity films Tu .

Exact calculations, i.e. solving the Eliashberg equations for the effective two-band model with parameters derived from first-principle electronic structure calculations, have been carried out for the clean and dirty cases for a magnetic field along the $`c`$-direction or in the $`ab`$-plane. The results obtained can be presented in the form of the effective superfluid plasma frequency, $`\omega _p^{sf}`$. Fig. 2 displays the calculated $`\left(\omega _p^{sf}(T)/\omega _p^{sf}(T0)\right)^2`$=$`\left(\lambda _L(T0)/\lambda _L(T)\right)^2`$ as a function of reduced temperature. First, we shall discuss the temperature dependence of $`\lambda _L^{ab}`$, when the magnetic field is oriented exactly along the $`c`$-axis (this means that screening currents run in the $`ab`$-plane). In the clean case the situation is similar to the model discussed above (Fig. 1). $`\lambda _L^{ab}(T)^2`$ has an inflection point and the low temperature behavior is determined by the band with the small gap $`\mathrm{\Delta }_\pi `$. In the dirty case the conductivity in the $`\pi `$-band is strongly suppressed. This means that the screening currents in the $`ab`$-plane are determined by $`\sigma `$-band with a BCS-like temperature dependence with a large gap $`\mathrm{\Delta }_\sigma `$. For the intermediate case the temperature dependence of $`\lambda _L^{ab}(T)`$ is between these limiting cases. One can even have situations with a nearly linear dependence in some temperature interval, as may be seen in Fig. 2. Experimental data from microwave experiments on single crystals Manzano  and oriented films Jin  as well as $`\mu SR`$ data on polycrystals Nieder  are shown for comparison in Fig. 2.

As may be seen from Eq. (2), the penetration depth in the $`c`$-direction in the clean case is only determined by the $`\pi `$-bands because of the very small plasma frequency of the $`\sigma `$-band in this direction. It is interesting to notice that the temperature dependence of the penetration depth in the dirty $`c`$-axis case follows the clean $`ab`$-direction case. This is because both the $`\sigma `$-channel is blocked by the small plasma frequency and the $`\pi `$-channel is blocked by impurity scattering. However, the absolute values of the penetration depths differ by a factor of about 8 in these two cases. One may see that an inflection point is present in the temperature dependence even for this single-band contribution, because of the induced superconductivity at higher temperatures in the $`\pi `$-band.

The corresponding London penetration depths at $`T`$=4 K have the values $`\lambda _{L,ab}^{clean}`$=39.2 nm, $`\lambda _{L,c}^{clean}`$=39.7 nm, $`\lambda _{L,ab}^{dirty}=105.7`$ nm, and $`\lambda _{L,c}^{dirty}`$=316.5 nm. One may observe that even in the clean case the value of $`\omega _{p,\alpha \beta }^{sf}(T0)5`$eV from Eq. 2 differs from the total plasma frequency $`7`$ eV as a consequence of strong-coupling effects due to electron-phonon interaction.

Table 1 summarizes the experimental information on the penetration depth $`\lambda _L(T0)`$, the temperature dependence of the penetration depth $`\mathrm{\Delta }\lambda _L(T)`$ and the estimated superconducting gap obtained by different experimental methods and groups. Our theoretical values of the penetration depths for the clean case are smaller than the smallest experimental value. On the other hand the values for the dirty case are in reasonable agreement with experiment. Nearly all measured penetration depths fall within our limiting cases (clean and dirty), and especially the BCS-like behavior at lower temperatures observed in experiment, reflecting the $`\pi `$-band contribution, is in agreement with our theoretical calculations.

It is well known that in a clean superconductor the low temperature penetration depth is independent of the superconducting gap. In this case, the anisotropy is only determined by the ratio of the plasma frequencies in the $`ab`$-plane and $`c`$-directions. Hence, if the $`\pi `$-band is very clean the penetration depth is nearly isotropic due to the small difference in the effective plasma frequencies. This is shown in Fig. 3, where one may also see that impurities in the $`\pi `$-band drastically enhance the anisotropy. In the inset of Fig. 3 the temperature dependence of the anisotropy is shown for the clean case. The reason for the strong variation with temperature is that for high temperatures the difference of gaps also contributes. A similar observation has recently been made in a weak coupling modelkogan .

A final remark concerns the orientation of the magnetic field. According to estimates in Ref. Brink,  the $`\sigma `$-band does not contribute to the electronic transport for angles with the $`c`$-axis larger than of the order of 1. This implies that for larger angles the effective penetration depth is determined by the $`\pi `$-band only. Only for angles approaching zero, the $`\sigma `$-band contributes and the penetration depth decreases towards the minimal value corresponding to the screening current flowing in the $`ab`$ plane, namely 39 nm and 106 nm for the clean and dirty case respectively.

The above mentioned considerations must be taken into account when interpreting the experimental data. In polycrystalline samples the penetration depth is mostly determined by the $`\pi `$-band and therefore practically isotropic and it is similar to the $`c`$-axis penetration depth in Fig. 2. Our calculations for the dirty case describes qualitatively well the data in Ref. Nieder, . On the other hand, the data for single crystals and oriented films correspond to our calculations of the $`ab`$-plane penetration depth, provided the magnetic field is oriented with an accuracy better than 1 along the $`c`$-axis. Therefore the data from Ref. Manzano,  and Jin,  are described by our calculation for the clean case. Quantitative deviations can be attributed to a different impurity content and a possible admixture of $`c`$-axis contribution. The temperature dependence of the specific electrical resistivity is not provided in these papers however, which would be needed in order to estimate the impurity scattering rates.

In conclusion, we have used the results of first principles calculations of the electronic structure and electron-phonon interaction in MgB2 to calculate the magnetic field penetration depth. The measured temperature dependence of the penetration depth is qualitatively well reproduced in a two band model with the same set of parameters which was used to fit DC resistivity wires ; Tu . We predict strong dependence of the anisotropy of the penetration depth on impurity scattering. This anisotropy increases with increasing temperature.

We acknowledge many useful discussions with I.I. Mazin and A.D. Caplin. J.K. would like to thank the Schloeßmann Foundation for financial support. This work was supported by the Dutch Foundation for Research on Matter (FOM).