text:
# Effect of realistic astrophysical inputs on the phase and shape of the WIMP annual modulation signal

## I Introduction

Arguably the best motivated non-baryonic dark matter candidate is the neutralino (the lightest supersymmetric particle), and direct detection experiments are just reaching the sensitivity required to probe the relevant region of parameter space . Since the expected event rates are so small ( $`𝒪(10^510)`$ counts $`\mathrm{kg}^1\mathrm{day}^1`$ see e.g. Refs. ) distinguishing a putative Weakly Interacting Massive Particle (WIMP) signal from backgrounds due to, for instance, neutrons from cosmic-ray induced muons or natural radioactivity, is crucial. The Earth’s motion about the Sun provides two potential WIMP smoking guns: i) an annual modulation  and ii) a strong direction dependence  of the event rate. In principle the dependence of the differential event rate on the atomic mass of the detector (see e.g. Refs. ) is a third possibility, however this would require good control of systematics for detectors composed of different materials. While direction sensitive detectors probably offer the best long-term prospects for the unambiguous detection of WIMPs, and development of such a detector is underway , annual-modulation searches are already feasible using large detector masses .

If the local WIMP velocity distribution is isotropic then the annual modulation is roughly sinusoidal with a maximum in early June (when the component of the Earth’s velocity in the direction of the Sun’s motion is largest) and amplitude of order a few per-cent. The DAMA collaboration, using a detector consisting of NaI crystal scintillators, have detected an annual modulation with roughly these properties, which they interpret as a WIMP signal . Assuming a standard halo model with a Maxwellian velocity distribution and circular velocity $`v_\mathrm{c}=220\mathrm{km}\mathrm{s}^1`$, they find a best fit WIMP mass $`m_\chi =52`$ GeV and cross-section $`\zeta \sigma _\mathrm{p}=7.2\times 10^9\mathrm{nb}`$ ***Here $`\zeta =\rho _\chi /(0.3\mathrm{GeV}\mathrm{cm}^3)`$ parameterizes the uncertainty in the local WIMP density, $`\rho _\chi `$. with the 3-$`\sigma `$ allowed region encompassing masses and cross-sections in the range $`30\mathrm{GeV}<m_\chi <100\mathrm{GeV}`$ The lower limit $`m_\chi >30`$ GeV is imposed by hand and represents the, somewhat model dependent, limit on the neutralino mass from accelerator sparticle searches (e.g. ). and $`10^9\mathrm{nb}<\zeta \sigma _\mathrm{p}<10^8\mathrm{nb}`$ .

Taken at face value this allowed region is incompatible with the exclusion limits from the Cryogenic Dark Matter Search (CDMS) , Edelweiss  and Zeplin I  experiments. While these exclusion limits depend relatively weakly on the WIMP velocity distribution (varying by of order tens of per-cent for fixed $`v_\mathrm{c}`$ ) the annual modulation signal depends sensitively on the halo model assumed  and the region of WIMP mass–cross-section parameter space corresponding to the DAMA signal may be significantly enlarged if non-standard halo models are considered . In particular Belli et. al.  carried out an extended analysis of the DAMA data, for a large range of halo models and parameters, and found a significant enlargement of the allowed region to encompasses masses in the range $`30\mathrm{GeV}<m_\chi <270\mathrm{GeV}`$ and cross-sections in the range $`10^{10}\mathrm{nb}<\zeta \sigma _\mathrm{p}<6\times 10^8\mathrm{nb}`$. However, as recently pointed out by Copi and Krauss  and Fornengo and Scopel , the phase of the annual modulation can be significantly different if the WIMP velocity distribution is anisotropic and the time variation may not be close to sinusoidal . In fact, in apparent contradiction to the results of Belli et. al. , Copi and Krauss  found, considering a similar range of halo models, that the DAMA annual modulation signal could not be made compatible with the exclusion limits from the CDMS and Edelweiss experiments.

Motivated by this apparent conflict, and also the possibility that if a WIMP annual modulation signal were detected its phase and shape  could allow us to probe the structure of the Milky Way (MW) halo, we study in detail the astrophysical uncertainties in the calculation of the annual modulation signal. In Section II we concentrate on the motion of the Earth relative to the Galactic rest frame. In Section III we discuss the local WIMP velocity distribution, in particular the importance of using physically and observationally reasonable models. In Section IV we examine the effects of the modeling of the Earth’s motion and the local velocity distribution on the phase and amplitude of the signal and the extent to which it can be approximated by a sinusoidal variation, before concluding with discussion in Section V.

## II Motion of the Earth

The differential WIMP elastic scattering rate due to scalar interactions is given by :

$$\frac{\mathrm{d}R}{\mathrm{d}E}(E,t)=\zeta \sigma _\mathrm{p}\left[\frac{\rho _{0.3}}{\sqrt{\pi }v_c}\frac{(m_\mathrm{p}+m_\chi )^2}{m_\mathrm{p}^2m_\chi ^3}A^2T(E,t)F^2(q)\right],$$
(1)
where $`\rho _{0.3}=0.3\mathrm{GeV}\mathrm{cm}^3`$ (and $`\zeta `$ is defined so that the local WIMP density is $`\zeta \rho _{0.3}`$), $`\sigma _\mathrm{p}`$ is the WIMP scattering cross section on the proton, $`v_\mathrm{c}`$ is the local circular velocity, $`A`$ and $`F(q)`$ are the mass number and form factor of the target nuclei respectively, $`E`$ is the recoil energy of the detector nucleus, and $`T(E)`$ is defined, so as to be dimensionless, as

$$T(E,t)=\frac{\sqrt{\pi }v_c}{2}_{v_{\mathrm{min}}}^{\mathrm{}}\frac{f_v(t)}{v}dv,$$
(2)
where $`f_v(t)`$ is the WIMP speed distribution in the rest frame of the detector, normalized to unity, and $`v_{\mathrm{min}}`$ is the minimum WIMP speed that can cause a recoil of energy $`E`$:

$$v_{\mathrm{min}}=\left(\frac{E(m_\chi +m_A)^2}{2m_\chi ^2m_A}\right)^{1/2}.$$
(3)
where $`m_A`$ is the Atomic mass of the detector nuclei.

The WIMP velocity distribution in the rest frame of the detector is found by making a, time dependent, Galilean transformation $`𝐯\stackrel{~}{𝐯}=𝐯+𝐯_\mathrm{e}(t)`$, where $`𝐯_\mathrm{e}(t)`$ is the Earth’s velocity relative to the Galactic rest frame, which has two components: the Earth’s orbit about the Sun and the Sun’s motion with respect to the Galactic rest frame.

### A Orbit about the Sun

The Earth moves in a close to circular orbit, inclined at an angle of roughly $`60^{}`$ to the Galactic plane, with orbital speed $`v_\mathrm{e}=29.79\mathrm{km}\mathrm{s}^1`$. The annual modulation is mainly determined by the component of the Earth’s motion in the direction of the Sun’s orbit For simplicity we ignore the Sun’s motion with respect to the local standard of rest in this sub-section.:

$$v_{\mathrm{e},}(t)=v_{}+v_\mathrm{e}\mathrm{sin}\gamma \mathrm{cos}\alpha (t),$$
(4)
where $`\alpha (t)=2\pi (tt_0)/T`$, $`T=1`$ year, $`t_0153`$ days (June 2nd), $`\gamma 30^{}`$ is the angle between the axis of the ecliptic and the galactic plane and $`t`$ is measured in days. Since $`v_{}v_e`$, if the velocity distribution is close to isotropic, the differential event rate can be expanded in a Taylor series in $`\mathrm{cos}\alpha (t)`$ so that, to first order, :

$$\frac{\mathrm{d}R}{\mathrm{d}E}(E,t)R_0(E)+R_1(E)\mathrm{cos}\alpha (t).$$
(5)
with $`R_0(E)R_1(E)`$, i.e. the modulation is roughly sinusoidal.

A commonly used  expression for the Earth’s motion, in Galactic co-ordinates (X,Y,Z) where X is toward the Galactic center, Y in the direction of rotation and Z toward the north Galactic pole, is:

$$𝐯_\mathrm{e}(t)=v_\mathrm{e}[\mathrm{sin}\alpha (t),\mathrm{cos}\alpha (t)\mathrm{sin}\gamma ,\mathrm{cos}\alpha (t)\mathrm{cos}\gamma ].$$
(6)
This expression assumes that the Earth’s orbit is circular and the axis of the ecliptic lies in the Y-Z plane .

A more accurate expression for the components of the Earth’s velocity relative to the Sun can be found by using the expression for the Sun’s ecliptic longitude (from e.g. p77 of Ref. ), $`\lambda =L+1^{}.915\mathrm{sin}g+0.020^{}\mathrm{sin}2g`$ where $`L=280^{}.460+0^{}.9856003t`$ is the mean longitude of the sun, corrected for aberration, $`g=357^{}.528+0^{}.9856003t`$ is the mean anomaly (polar angle of orbit) and $`t`$ is the time in days from 1200 Universal Time §§§Universal Time is equivalent to Greenwich Mean Time to within an accuracy of seconds. on Dec 31st 1999, and transforming from ecliptic to Galactic co-ordinates (e.g. p13 of Ref. ). Lewin and Smith  carried out his procedure and found:

$`𝐯_\mathrm{e}(t)`$ $`=`$ $`v_\mathrm{e}(\lambda )[\mathrm{cos}\beta _x\mathrm{sin}(\lambda \lambda _x),`$ (8)
$`\mathrm{cos}\beta _y\mathrm{sin}(\lambda \lambda _y),\mathrm{cos}\beta _z\mathrm{sin}(\lambda \lambda _z)],`$
where $`v_\mathrm{e}(\lambda )=v_\mathrm{e}\left[1e\mathrm{sin}(\lambda \lambda _0)\right]`$, $`e=0.016722`$ and $`\lambda _0=13^{}\pm 1^{}`$ are the ellipticity of the Earth’s orbit and the ecliptic longitude of the orbit’s minor axis respectively, and $`\beta _i=(5^{}.5303,59^{}.575,29^{}.812)`$ and $`\lambda _i=(266^{}.141,13^{}.3485,179^{}.3212)`$ are the ecliptic latitudes and longitudes of the of the (X,Y,Z) axes respectively.

Fornengo and Scopel  used a similar expression where the the non-uniform motion of the Sun in right ascension is neglected in the co-ordinate transformation so that the $`\mathrm{sin}(\lambda \lambda _i)`$ terms in eq. (8) above are replaced by $`\mathrm{cos}\left[\omega (tt_i)\right]`$ where $`t_i=(76.1,156.3,352.4)`$ days and $`\omega =2\pi /(1\mathrm{year})`$. Gelmini and Gondolo  used a slightly simpler expression, which neglects the ellipticity of the Earth’s orbit and the non-uniform motion of the Sun in right ascension. In terms of the Sun’s ecliptic longitude, which under these approximations is given by $`\stackrel{~}{\lambda }(t)=\omega (tt_1)`$ where $`t_1`$ is the spring (vernal) equinox, $`t_1=79.55`$ days This figure is for 2003 and is accurate to within about 0.04 days. The time at which the Spring equinox occurs increases by 0.24 days a year (see e.g. Ref. ).,:

$$𝐯_\mathrm{e}(t)=v_\mathrm{e}\left[\widehat{𝐞}_1\mathrm{sin}\stackrel{~}{\lambda }(t)\widehat{𝐞}_2\mathrm{cos}\stackrel{~}{\lambda }(t)\right],$$
(9)
where $`\widehat{𝐞}_1`$ and $`\widehat{𝐞}_1`$ are in the Sun’s direction at the Spring equinox and Summer solstice respectively:

$`\widehat{𝐞}_1`$ $`=`$ $`(0.0670,\mathrm{\hspace{0.17em}0.4927},0.8676),`$ (10)
$`\widehat{𝐞}_2`$ $`=`$ $`(0.9931,0.1170,\mathrm{\hspace{0.17em}0.01032}).`$ (11)
These expressions for the components of the Earth’s orbital velocity are displayed in Fig. 1, and the deviations between the expressions in Fig. 2. In the plots time is measured in days from noon on Dec 31st 2002, but the correct time (measured in days from noon on Dec 31st 1999) is used in the calculations. If time is erroneously measured from the beginning of the year rather than the beginning of 2000, then this leads to a horizontal shift in the curves of 0.27 days for every year elapsed since 2000. The deviations in the Y component are reassuring small (of order a few per-cent), however assuming that the Earth’s orbit is circular and the axis of the ecliptic lies in the Y-Z plane (eq. (6)) leads to large errors in the X and Z components. The other three expressions are in relatively good agreement for all three components.

### B Motion of the Sun

We now consider the motion of the Sun, which can be divided into two components: the motion of the local standard of rest (LSR) and the Sun’s peculiar motion with respect to the LSR, $`v_{,\mathrm{pec}}`$. The “standard values” often used in WIMP event rate calculations date back to Kerr and Lynden-Bell’s 1986 “Review of Galactic constants” , however since then the Hipparcos satellite has provided accurate measurements of the motions of large numbers of nearby stars, allowing more accurate determinations of the Galactic constants, so the “standard values” should be updated.

The standard value used for the Sun’s peculiar motion, with respect to the LSR, is $`v_{,\mathrm{pec}}=(10,15,8)\mathrm{km}\mathrm{s}^1`$ , while Lewin and Smith  use $`v_{,\mathrm{pec}}=(9,12,7)\mathrm{km}\mathrm{s}^1`$, with errors of order a few $`\mathrm{km}\mathrm{s}^1`$ in each component leading to an uncertainty of several days in the phase of the modulation . The value determined more recently using data from Hipparcos is $`v_{,\mathrm{pec}}=(10.0,5.2,7.2)\mathrm{km}\mathrm{s}^1`$ , with the stated errors in each component being roughly $`0.5\mathrm{km}\mathrm{s}^1`$.

Assuming that the MW is axisymmetric, then the motion of the LSR is, by definition, $`(0,v_\mathrm{c}(R_0),0)`$ where $`v_\mathrm{c}(R_0)`$ is the circular velocity at the solar radius. Kerr and Lynden-Bell combined a large number of independent determinations of the circular velocity and found $`v_\mathrm{c}=222.2\mathrm{km}\mathrm{s}^1`$, with standard deviation $`20\mathrm{km}\mathrm{s}^1`$ . The proper motions of Cepheids measured by Hipparcos allow accurate determinations of the Oort constants $`A`$ and $`B`$ , which lead to a value for the circular velocity, in terms of the Galactocentric distance $`R_0`$: $`v_\mathrm{c}(R_0)=(27.2\pm 0.9)(R_0/\mathrm{kpc})\mathrm{km}\mathrm{s}^1`$ . One possible caution is that this calculation assumes that stars move on circular orbits, which may lead to a systematic error as the MW is not exactly axisymmetric . Accurate very long baseline interferometry measurements of the proper motion of the $`\mathrm{SgrA}^{}`$ radio source at the Galactic center provides an alternative, marginally inconsistent, determination of $`v_\mathrm{c}`$: $`(30.1\pm 0.9)(R_0/\mathrm{kpc})\mathrm{km}\mathrm{s}^1`$ . In this case non-zero peculiar motion of $`\mathrm{SgrA}^{}`$ produces a systematic error. Using $`R_0=8.0\pm 0.5\mathrm{kpc}`$ from the most recent compilation of determinations of $`R_0`$ , leads to $`v_\mathrm{c}(R_0)=(218\pm 15)\mathrm{km}\mathrm{s}^1`$ from Cepheids and $`v_\mathrm{c}(R_0)=(241\pm 17)\mathrm{km}\mathrm{s}^1`$ from $`\mathrm{SgrA}^{}`$.

## III Local WIMP velocity distribution

Data analyzes nearly always assume a standard halo model with an isotropic Maxwellian velocity distribution:

$$f(𝐯)=\frac{1}{\pi ^{3/2}v_{\mathrm{c},\mathrm{h}}^3}\mathrm{exp}(𝐯^2/v_{\mathrm{c},\mathrm{h}}^2),$$
(12)
where $`v_{\mathrm{c},\mathrm{h}}`$ is the contribution of the halo to the local circular velocity, which corresponds to a spherically symmetric density distribution with $`\rho r^2`$ i.e. an isothermal sphere. Observations and numerical simulations indicate that galaxy halos are in fact triaxial and anisotropic however (see Ref.  for a review). As the size  and phase  of the annual modulation signal depend quite sensitively on the halo model assumed, the realistic modeling of the WIMP velocity distribution is crucial when extracting parameters or exclusion limits from data or comparing results from different experiments.

The steady state phase-space distribution function of a collection of collisionless particles is given by the collisionless Boltzmann equation, and the velocity dispersions, $`<v_\mathrm{i}^2>`$, of the system are calculated via the Jeans equations, which are found by taking moments of the collisionless Boltzmann equation (see e.g. ). Solutions to the Jean’s equations have the property that the tensor $`\sigma _{\mathrm{ij}}^2\overline{(v_\mathrm{i}\overline{v}_\mathrm{i})(v_\mathrm{j}\overline{v}_\mathrm{j})}`$ is symmetric, so that at any point a set of orthogonal axes can by chosen such that $`\sigma `$ is diagonal (i.e. $`\sigma _{\mathrm{ij}}^2=\sigma _{\mathrm{ii}}^2\delta _{\mathrm{ij}}`$. However, as discussed in Refs. , there is no equation of state relating the components of the velocity dispersion to the density, so solving the Jeans equations requires assumptions, which may or may not be physically reasonable, about the shape and/or orientation of $`\sigma `$. Evans et. al.  presented the logarithmic ellipsoidal model (which has potential $`\mathrm{\Phi }(x,y,z)=(v_\mathrm{c}^2/2)\mathrm{ln}(x^2+y^2p^2+z^2q^2)`$) which is the simplest triaxial, scale free generalization of the isothermal sphere. They argued that principle axes aligned with conical co-ordinates correspond to physical distribution functions and calculated the corresponding velocity dispersions. On the axes of the halo conical co-ordinates are locally equivalent to cylindrical polar co-ordinates and the local velocity distribution can be approximated by a multi-variate Gaussian with principle axes aligned with the (X,Y,Z) Galactic axes. This form has been widely used in calculations of the WIMP direction detection rate , however a multi-variate Gaussian with arbitrary velocity dispersions need not correspond to a physically sensible halo model. For instance Evans et. al. do not consider parameter values for which the ratio of any two of the velocity dispersions is more than 3:1, so as to avoid models which are afflicted by instabilities .

We will now briefly examine what sets of values of the velocity dispersions are reasonable for the logarithmic ellipsoidal model, if the Sun is located on the intermediate axis of the halo (the deviations from an isotropic velocity distribution are smaller on the major axis ). On the intermediate axis the velocity distribution can be written as

$$f(𝐯)=\frac{1}{(2\pi )^{3/2}\sigma _\mathrm{r}\sigma _\varphi \sigma _\mathrm{z}}\mathrm{exp}\left(\frac{v_\mathrm{r}^2}{2\sigma _\mathrm{r}^2}\frac{v_\varphi ^2}{2\sigma _\varphi ^2}\frac{v_\mathrm{z}^2}{2\sigma _\mathrm{z}^2}\right),$$
(13)
with the velocity dispersions given by

$`\sigma _\mathrm{r}^2`$ $`=`$ $`{\displaystyle \frac{v_{\mathrm{c},\mathrm{h}}^2p^4}{(2+\gamma )(1p^2+q^2)}},`$ (14)
$`\sigma _\varphi ^2`$ $`=`$ $`{\displaystyle \frac{v_{\mathrm{c},\mathrm{h}}^2(2q^2p^2)}{2(1p^2+q^2)}},`$ (15)
$`\sigma _\mathrm{z}^2`$ $`=`$ $`{\displaystyle \frac{v_{\mathrm{c},\mathrm{h}}^2(2p^2)}{2(1p^2+q^2)}},`$ (16)
where $`p`$ and $`q`$ are constants which satisfy $`0qp1`$ and are related to the axial ratios of the density distribution, $`I_{1,2}`$, by

$`I_1^2`$ $`=`$ $`{\displaystyle \frac{p^2(p^2q^2+p^2q^2)}{q^2+p^2p^2q^2}},`$ (17)
$`I_2^2`$ $`=`$ $`{\displaystyle \frac{q^2(p^2q^2p^2+q^2)}{q^2+p^2p^2q^2}},`$ (18)
and $`\gamma `$ is a constant isotropy parameter, which in the spherical limit $`p=q=1`$ is related to the anisotropy parameter

$$\beta =1\frac{<v_\theta ^2>+<v_\varphi ^2>}{2<v_\mathrm{r}^2>},$$
(19)
by $`\gamma =2\beta `$. Note that $`\sigma _\varphi ^2+\sigma _\mathrm{z}^2=v_{\mathrm{c},\mathrm{h}}^2`$ and $`\sigma _\varphi >\sigma _\mathrm{z}`$.

If we require $`\sigma _\mathrm{j}/3<\sigma _\mathrm{i}<3\sigma _\mathrm{j}`$, so as to avoid models which might be affected by instabilities , then the velocity dispersions lie in the ranges $`0.32<\sigma _\mathrm{z}/v_{\mathrm{c},\mathrm{h}}<0.71`$, $`0.71<\sigma _\varphi /v_{\mathrm{c},\mathrm{h}}<0.95`$ and $`0.24<\sigma _\mathrm{r}/v_{\mathrm{c},\mathrm{h}}<2.0`$. If we also require that $`0<\beta <0.4`$, as found at the solar radius in simulated halos , then $`\sigma _\mathrm{r}`$ is further restricted to $`0.71<\sigma _\mathrm{r}/v_{\mathrm{c},\mathrm{h}}<0.92`$. In Fig. 3 we plot the values of $`\sigma _\mathrm{r}`$ and $`\sigma _\varphi `$ ($`\sigma _\mathrm{z}^2=v_{\mathrm{c},\mathrm{h}}^2\sigma _\varphi ^2`$) which satisfy these two requirements. We still need to check that the axial ratios of the density distribution are in reasonable agreement with observed and simulated galaxy halos. There is no simple algebraic relation between the $`I_\mathrm{i}`$ s and the $`\sigma _\mathrm{i}`$ s. If we require $`0.6<I_1<1.0`$ and $`0.3<I_2<1.0`$ (see e.g. Ref.  for a review of constraints on the shape of dark matter halos) then some of the previously allowed sets of velocity dispersions are now excluded.

The logarithmic ellipsoidal model may not be a unique solution to the Jeans equations (and there is no reason to expect the Sun to be located on one of the axes of the halo), but otherwise there is no reason for the principle axes of the velocity distribution to correspond to the axes of the halo. In other words a multivariate Gaussian velocity distribution with axes corresponding to the axes of the halo and arbitrary velocity dispersions may not correspond to a physically reasonable halo model. While the constraints on the axial ratios and anisotropy of the halo we have imposed are not cast in stone, they illustrate that only restricted sets of values of the velocity dispersions correspond to observationally and physically reasonable halo models, and that the ratios of the velocity dispersions can not be too large. Simulations provide further support for this argument. Helmi, White and Springel  examined the simulation particles within a 4 $`\mathrm{kpc}`$ box located 8 $`\mathrm{kpc}`$ from the center of a Milky Way like halo and found that, apart from a clump of fast moving particles from a late accreting subhalo, the velocity distribution was well approximated by a multi-variate Gaussian, with principal axis velocity dispersions in the ratio $`1:1.08:1.27`$. If in fact the local dark matter distribution is not smooth (see below for discussion of this possibility), then the local velocity distribution could not be approximated by a multivariate Gaussian with any set of velocity dispersions. We conclude that the large velocity dispersion ratios found by Fornengo and Scopel to produce extreme distortions of the annual modulation ($`(1:5:4)`$ and $`(10:1:3)`$ are unlikely to correspond to a realistic halo model.

It is possible that the local WIMP distribution may not be completely smooth, or in other words that the phase-space distribution has not reached a steady state. The density and velocity distributions of the “particles” in numerically simulated halos are relatively smooth at the solar radius , however WIMP direct detection experiments probe the dark matter distribution within a local sub-mpc region and even the highest resolution simulations have less than 100 “particles” per $`\mathrm{kpc}^3`$. Furthermore the first neutralino clumps to form have mass more than ten orders of magnitude smaller than the smallest subhalos resolved in numerical simulations , and it is possible that the cores of some of these first small, high density, clumps could survive to the present day . Stiff and Widrow  have recently developed a technique which uses test particles to probe the velocity distribution at a single point within a numerical simulation. They find a local velocity distribution consisting of a series of discrete peaks, and argue that the smooth background found previously in simulations may be, at least partly, a consequence of the use of a finite volume to measure the velocity distribution. Even if the local phase-space distribution bears no trace of the first WIMP clumps to form, more massive subhalos which have been accreted onto the MW relatively recently, and have sufficiently eccentric orbits, could produce coherent streams of high velocity WIMPs at the solar radius .

An additional complication is that the contribution of the visible components of the MW to the circular velocity at the Solar radius is non-negligible, and may even dominate that of the halo (see e.g. Refs. ). Self consistent halo models can in principle be built (taking into account various observations) using Eddington’s formula , however building mass models of the Milky Way is a complex and ill-constrained process ; a small change in the observational constraints can lead to a large change in the properties of the dark halo. In Dehnen and Binney’s mass models, where the halo is described by a spheroidal density distribution, the contribution of the MW halo to the circular velocity at the solar radius lies in the range $`v_{\mathrm{c},\mathrm{h}}(R_0)=110180\mathrm{km}\mathrm{s}^1`$ , while Moore et. al. find, for halos with a central density cusp $`\rho r^{1.5}`$, taking $`190\mathrm{km}\mathrm{s}^1<v_\mathrm{c}(R_0)<230\mathrm{km}\mathrm{s}^1`$ and using a smaller set of constraints, $`100\mathrm{km}\mathrm{s}^1<v_{\mathrm{c},\mathrm{h}}(R_0)<120\mathrm{km}\mathrm{s}^1`$ . We therefore conclude that a relatively conservative range of values is: $`100\mathrm{km}\mathrm{s}^1<v_{\mathrm{c},\mathrm{h}}(R_0)<200\mathrm{km}\mathrm{s}^1`$, and that the “standard” value, $`v_{\mathrm{c},\mathrm{h}}(R_0)=220\mathrm{km}\mathrm{s}^1`$, is probably too large.

## IV Results

To compare theoretical predictions with experimental data requires the detector atomic mass, form factor and resolution, and the relationship between the nuclear recoil energy and the energy detected, to be taken into account (see e.g. Ref. ). As the DAMA collaboration’s raw data is not publicly available instead of working with the differential event rate we look at the detector independent quantity $`T(v_{\mathrm{min}},t)`$, as defined in eq. (2). The relation between $`v_{\mathrm{min}}`$ and the detected energy is given, for various target nuclei and WIMP masses, in Fig. 4 of Ref. . For reference, the DAMA experiment has a measured energy threshold of 2 keV which corresponds (for Iodine which dominates) to a recoil energy threshold of 22 keV, which for $`m_\chi =50(200)`$ GeV requires a minimum WIMP velocity $`v_{\mathrm{min}}=309(146)\mathrm{km}\mathrm{s}^1`$.

### A Motion of the Earth

We will start by examining the effects of the modeling of the motion of the Earth on the phase and shape of the annual modulation signal. We use two fiducial halo models: the standard halo model (with an isotropic Maxwellian velocity distribution) and the logarithmic ellipsoidal model , with the Earth located on the intermediate axis and parameter values $`p=0.9,q=0.8`$ and $`\gamma =0.62`$ corresponding to axial ratios $`1:0.78:0.48`$ and $`\beta =0.4`$. The later model is intended as a specific, somewhat arbitrary but not unreasonable, example of an anisotropic halo model. In both cases we take $`v_{\mathrm{c},\mathrm{h}}=150\mathrm{km}\mathrm{s}^1`$. In Fig. 4 we plot $`T(v_{\mathrm{min}},t)`$ Our plots look superficially different to those of Fornengo and Scopel as they normalize to the event rate on January 1st whereas we subtract the mean event rate. produced by each of the expressions for the Earth’s motion discussed in Sec. II, for five specific values of $`v_{\mathrm{min}}`$. As has long been known (see e.g. Ref. ), for small $`v_{\mathrm{min}}`$ the phase of the modulation is reversed. Furthermore as $`v_{\mathrm{min}}`$ is increased the amplitude of the modulation reaches a maximum before decreasing again. For the standard halo model this maximum occurs at roughly $`220\mathrm{km}\mathrm{s}^1`$, whereas for the fiducial triaxial model the maximum is at around $`250\mathrm{km}\mathrm{s}^1`$ and $`T(\mathrm{\hspace{0.17em}200}\mathrm{km}\mathrm{s}^1,t)T(\mathrm{\hspace{0.17em}300}\mathrm{km}\mathrm{s}^1,t)`$. In Fig. 5 we plot $`T(v_{\mathrm{min}},t)`$ produced by each of the sets of values for the Sun’s motion with respect to the LSR discussed in Sec. II.

To provide a quantitative comparison tables I and II contain the day at which $`T(v_{\mathrm{min}},t)`$ (and hence the differential event rate) is largest, $`t_\mathrm{p}`$, the maximum deviation from $`T(v_{\mathrm{min}},t)`$ calculated using the full expression for the Earth’s motion (eq. (8)), and the maximum deviation from a sinusoidally varying $`T(v_{\mathrm{min}},t)`$ with the same mean, amplitude and $`t_\mathrm{p}`$, for the standard halo model and the fiducial triaxial model respectively. Tables III and IV contain the same data for $`T(v_{\mathrm{min}},t)`$ found using each of the sets of values for the Sun’s velocity with respect to the LSR. For reference the DAMA collaboration have carried out a fit to the modulation of their data using the function $`A\mathrm{cos}\left[2\pi (tt_\mathrm{p})/T\right]`$, and when $`T`$ is fixed at one year they found $`t_\mathrm{p}=144\pm 13`$ days .

The time at which the event rate is maximum \**\**\**For simplicity we neglect the “flip” in the phase which occurs at small $`v_{\mathrm{min}}`$ in this discussion. varies by up to ten days depending on the expression used for the Earth’s velocity, and depends on $`v_{\mathrm{min}}`$ if the velocity distribution is anisotropic . If only the component of the Earth’s velocity in the Y direction is used, then the change in the phase which occurs for anisotropic halos is missed. The three sophisticated expressions produce results which are in reasonably good agreement however, with maximum deviations of around a few per-cent for $`v_{\mathrm{min}}400\mathrm{km}\mathrm{s}^1`$. For larger $`v_{\mathrm{min}}`$ the exact form of the high energy tail of the speed distribution becomes important, so that using a different expression for the Earth’s velocity produces a large fractional change in $`T(v_{\mathrm{min}},t)`$. As the mean event rate is small for large $`v_{\mathrm{min}}`$, the absolute errors are small though and therefore not so important from an experimental point of view.

The maximum deviations from a sinusoidal $`T(v_{\mathrm{min}},t)`$ are no more than a few per-cent for $`v_{\mathrm{min}}400\mathrm{km}\mathrm{s}^1`$, as is expected since the error in neglecting the 2nd and higher order terms in the Taylor expansion is $`𝒪((v_{\mathrm{e},\mathrm{Y}}/v_{})^20.01)`$. The amplitude of the annual modulation is largest for large $`v_{\mathrm{min}}`$ , as in the tail of the speed distribution the fraction of particles with speed greater than some fixed value changes significantly with $`v_\mathrm{e}(t)`$, the Taylor expansion is then inappropriate and the deviation from sinusoidal variation becomes large (up to around $`10\%`$). For the fiducial anisotropic model we have chosen the deviations happen to be smaller than for the standard isotropic halo model, but this is not the case in general.

The time at which the event rate is maximum is relatively weakly dependent on the Sun’s velocity with respect to the LSR (provided that it isn’t neglected entirely), however the maximum deviation from $`T(v_{\mathrm{min}},t)`$ found using the Hipparcos values for the Sun’s velocity depends strongly on the component in the Y direction and can be large. Finally we note that combining a poor approximation for the Earth’s orbit with an erroneous value for the Sun’s velocity with respect to the LSR would produce even larger errors.

### B Local WIMP velocity distribution

The region of WIMP mass-cross-section parameter space consistent with the DAMA annual modulation signal depends strongly on the contribution of the halo to the circular velocity at the solar radius, $`v_{\mathrm{c},\mathrm{h}}(R_0)`$ . To explicitly illustrate the reason for this, in Fig. 6 we plot $`T(v_{\mathrm{min}},t)`$ for $`v_{\mathrm{min}}=100,200`$ and $`300\mathrm{km}\mathrm{s}^1`$ for values of $`v_{\mathrm{c},\mathrm{h}}`$ in the range $`100200\mathrm{km}\mathrm{s}^1`$, for the standard halo model. As $`v_{\mathrm{c},\mathrm{h}}`$ (and hence the typical WIMP velocity) is decreased the value of $`v_{\mathrm{min}}`$ at which the modulation flips phase and also that at which the amplitude of the modulation reaches a local maximum, are both smaller. For instance for $`v_{\mathrm{c},\mathrm{h}}=100\mathrm{km}\mathrm{s}^1`$ the local maximum occurs at roughly $`v_{\mathrm{min}}=155\mathrm{km}\mathrm{s}^1`$ compared with $`v_{\mathrm{min}}=220\mathrm{km}\mathrm{s}^1`$ for $`v_{\mathrm{c},\mathrm{h}}=150\mathrm{km}\mathrm{s}^1`$.

As our aim is to examine how physical and observational constraints restrict changes in the phase and shape of the annual modulation signal (rather than to carry out a detailed comparison with experimental data) we focus on the logarithmic ellipsoidal model, which reproduces some of the important features (triaxiality and velocity anisotropy) of real galaxy halos. We examine the form of $`T(v_{\mathrm{min}},t)`$, with the Sun located on the intermediate axis, for the sets of values of the velocity dispersions that are unlikely to lead to instabilities and those which also correspond to halos with axis ratios and velocity anisotropy consistent with observed and simulated halos, as discussed in Sec. III. In Fig. 7 we plot the shift in the time at which $`T(v_{\mathrm{min}},t)`$ is maximum, $`\mathrm{\Delta }t_\mathrm{p}`$, and the fractional change in the amplitude of the variation of $`T(v_{\mathrm{min}},t)`$ relative to that for the standard isotropic halo model, for $`v_{\mathrm{min}}=100,300`$ and $`500\mathrm{km}\mathrm{s}^1`$, fixing $`v_{\mathrm{c},\mathrm{h}}=150\mathrm{km}\mathrm{s}^1`$. In Fig. 8 we plot the fractional changes in the mean and amplitude of the variation of $`T(v_{\mathrm{min}},t)`$, relative to the standard isotropic halo model.

If we consider all sets of velocity dispersions for which $`\sigma _\mathrm{j}/3<\sigma _\mathrm{i}<3\sigma _\mathrm{j}`$, the shift in $`t_\mathrm{p}`$ has magnitude of up to 40 days and is quite strongly dependent on $`v_{\mathrm{min}}`$. The change in the amplitude is of order $`10`$s of per-cent for $`v_{\mathrm{min}}=100`$ and $`300\mathrm{km}\mathrm{s}^1`$ and up to a factor of 5 for $`v_{\mathrm{min}}=500\mathrm{km}\mathrm{s}^1`$, while the change in the mean is even larger increasing from $`10`$s of per-cent for $`v_{\mathrm{min}}=100\mathrm{km}\mathrm{s}^1`$ to more than an order of magnitude at $`v_{\mathrm{min}}=500\mathrm{km}\mathrm{s}^1`$. The large change in the mean and amplitude for large $`v_{\mathrm{min}}`$ occur because the logarithmic ellipsoidal model has a more extended tail of high velocity particles than the standard halo model  however, as noted in Sec. III, the mean event rate is very small for large $`v_{\mathrm{min}}`$. When only sets of velocity dispersions which correspond to halos with realistic axis ratios and velocity anisotropy are considered the maximum shift in $`t_\mathrm{p}`$ is roughly 50 $`\%`$ smaller and the the maximum change in the mean signal is roughly an order of magnitude smaller, which illustrates the importance of only considering velocity distributions that are physically and/or observationally reasonable. Note that the standard isotropic halo model (which has $`\beta =0`$ and $`I_1=I_2=1`$) lies at the edge of the observational constraints we impose.

In Figs. 9 and 10 we plot the shifts in the mean, amplitude and phase relative to the standard isotropic halo model with $`v_{\mathrm{c},\mathrm{h}}=150\mathrm{km}\mathrm{s}^1`$, for the same sets of velocity dispersion ratios but now with $`v_{\mathrm{c},\mathrm{h}}=200\mathrm{km}\mathrm{s}^1`$. The changes in the signal due to the change in $`v_{\mathrm{c},\mathrm{h}}`$ are of the same order of magnitude as the changes from varying the other parameters of the model.

We do not find shifts in $`t_\mathrm{p}`$ as large as those found by Copi and Krauss  and Fornengo and Scopel . For many of the sets of parameter values we have considered the shift in $`t_\mathrm{p}`$ is large enough to be incompatible with the phase of the DAMA annual modulation signal, however. This suggests that if all three components of the Earth’s velocity are included when calculating the event rate, then the consideration of anisotropic halo models may not lead to as large an increase in the region of WIMP mass and cross-section parameter space corresponding to the DAMA annual modulation signal as found in Ref. . As the DAMA data is not publicly available, however, it is not possible to confirm this.

Finally we overview the effect of streams of particles on the annual modulation signal. A stream of particles with negligible velocity dispersion would produce a contribution to $`T(v_{\mathrm{min}},t)`$ which is a step ††††††The velocity dispersion would in fact be finite, but small, and would lead to a broadening in the step ., the position and amplitude of which varies with time :

$$T_\mathrm{s}(v_{\mathrm{min}},t)=\{\begin{array}{cc}\frac{\sqrt{\pi }v_{\mathrm{c},\mathrm{h}}}{2}\frac{\rho _\mathrm{s}}{\zeta \rho _{0.3}}\frac{1}{v_\mathrm{s}(t)},\hfill & \mathrm{for}v_{\mathrm{min}}<v_\mathrm{s}(t),\hfill \\ 0,\hfill & \mathrm{for}v_{\mathrm{min}}>v_\mathrm{s}(t).\hfill \end{array}$$
(20)
where $`\rho _\mathrm{s}`$ is the density of the stream and $`v_\mathrm{s}(t)`$ is the speed of the stream in the rest frame of the detector, which is time dependent due to the time dependence of the Earth’s velocity. The variation of the amplitude and position of the step depend on the relative orientation of the Earth’s orbit and the orbit of the clump; for a clump of fixed velocity the closer the alignment of the clump’s path with the Earth’s orbit the larger the variation in $`v_\mathrm{s}(t)`$ (in the, extremely improbably, case that the clump’s path was perpendicular to the Earth’s orbit there would be no variation). As outlined in Ref.  if such a step were detected in the differential energy spectrum then the density, speed and direction of the stream responsible could, in principle (with a large number of events and a detector with good energy resolution), be recovered. Multiple streams would produce multiple steps at different energies with different temporal variations, which might be difficult to disentangle. Even if most of the dark matter is smoothly distributed a high velocity ($`v_\mathrm{s}500\mathrm{kms}^1`$) stream of particles from a late accreting clump could produce a step in the differential event rate at large energy, the position of which would be modulated as discussed above . However as the event rate decreases roughly exponentially with increasing energy this would only be detectable with very large statistics (large target mass and exposure). A similar signal would be produced if there is a population of extragalactic WIMPs .

## V Discussion

We have investigated the effect of uncertainties in astrophysical inputs (the motion of the Earth with respect to the Galactic rest frame and the local WIMP velocity distribution) on the calculation of the WIMP annual modulation signal. Accurate calculation of the shape and phase of the annual modulation signal requires all three components of the Earth’s velocity with respect to the Sun to be taken into account properly. If the components perpendicular to the Sun’s motion are neglected, the (energy dependent) shift in $`t_\mathrm{p}`$ which occurs if the local velocity distribution is anisotropic is missed. Neglecting the motion of the Sun with respect to the Local Standard of rest, $`v_{,\mathrm{pec}}`$, leads to an error of around 10 days in $`t_\mathrm{p}`$, and uncertainties in $`v_{,\mathrm{pec}}`$ lead to an error of order a few days in $`t_\mathrm{p}`$ (as stated in Ref. ) and errors in the shape of the signal which grow from a few percent at low energies to more than ten per-cent at high energy. Approximating the modulation with a sinusoid with the same mean, amplitude and phase produces errors of up to ten per-cent. The errors in $`t_\mathrm{p}`$ are crucial when comparing theoretical expectations with the modulation observed by the DAMA collaboration  (which has phase $`t_\mathrm{p}=144\pm 13`$ days). The errors in the shape of the signal induced by uncertainties in the Earth’s motion are currently less important, but may become important for annual modulation searches with tonne size detectors (such as the planned GENIUS detector ), especially if we want to extract information about the local velocity distribution from an observed signal.

If the local WIMP velocity is anisotropic then the phase , amplitude  and even shape  of the annual modulation signal can change. We have investigated, focusing on the logarithmic ellipsoidal halo model, how the form of the annual modulation changes for parameter choices which correspond to physically and observationally reasonable halo models. We found that for reasonable sets of values for the velocity dispersions the shift in $`t_\mathrm{p}`$, which is energy dependent, can be up to 20 days and the mean and amplitude of the signal change by tens of per-cent, at experimentally accessible energies. It is possible that other halo models could produce larger changes in the annual modulation signal, however we have shown, in the context of this model, that restricting the choice of parameter values to those that are physically and observationally reasonable seriously restricts the changes in the signal. We have also argued that a multivariate Gaussian velocity distribution, with axes corresponding to the axes of the halo, and velocity dispersions with large ratios is unlikely to correspond to a physically reasonable halo model.

It is often assumed that the contribution of the luminous components of the MW to the circular velocity at the solar radius is negligible so that $`v_{\mathrm{c},\mathrm{h}}(R_{})v_\mathrm{c}(R_{})220\mathrm{kms}^1`$. This is not the case and in fact there is a large uncertainty in the value of $`v_{\mathrm{c},\mathrm{h}}(R_{})`$ , which is important since (as we saw in Sec. IIIB) the mean and amplitude of the event rate (and hence the region of WIMP mass- cross-section parameters space corresponding to an observed modulation ) depend quite sensitively on $`v_{\mathrm{c},\mathrm{h}}(R_{})`$.

Finally we discussed the possibility that the local WIMP distribution is not completely smooth . If the local dark matter distribution consists of a number of streams of particles with small velocity dispersion, then the event rate would contain a number of steps, the amplitude and position of which would vary with time . Even if there is a smooth background WIMP distribution, high velocity streams from late accreting subhalos  may be detectable with good energy resolution and a large number of events.

In summary, it is important to properly take into account astrophysical uncertainties (not just in the WIMP velocity distribution, but also in the motion of the detector with respect to the Galactic rest frame) when calculating the WIMP annual modulation signal. Analyzing data assuming a sinusoidal modulation with fixed phase could lead to erroneous constraints on, or best fit values, for the WIMP mass or cross-section, even worse a WIMP signal could be overlooked. On the other hand using unrealistic halo models or parameter values could lead to overly restrictive exclusion limits or a misleadingly large range of allowed values of the WIMP mass and cross-section.

## Acknowledgments

A.M.G. was supported by the Swedish Research Council.

## A Erratum

While the halo only contributes some fraction of the circular velocity at the solar radius, it is the total circular velocity $`v_\mathrm{c}(R_0)=220\pm 20\mathrm{k}\mathrm{m}\mathrm{s}^1`$ (and not just the halo’s contribution) which determines the velocity dispersions of the dark matter particles; the dark matter particles feel the gravity of all the components of the MW. Some of the values of $`v_\mathrm{c}(R_0)`$ used for quantitative calculations in this paper are therefore too low. The main conclusions of this paper, regarding the effects of uncertainties in the modelling of the Earths’s motion and the local WIMP velocity distribution on the phase and shape of the annual modulation signal, are unchanged however.

I am grateful to Piero Ullio and Joakim Edsjö for bringing this issue to my attention.


format: text

text:
# Effect of symmetry of the electron states of HTSC on the current-voltage characteristics of SIS junctions.

## Abstract

The current-voltage $`IV`$ characteristics of $`SIS`$ junctions are calculated in the framework of a multiband model with an anisotropic effective order parameter of HTSC. The results of calculations show that the shape of $`IV`$ characteristic and the density of electron states changes significantly depending on the parameters of the model. A theoretical explanation is proposed for the experimentally observed $`s`$-like behavior of the $`IV`$ characteristics of $`SIN`$ and $`SIS`$ junctions with BSCCO-type superconductors. The dependence of the superconducting peaks asymmetry on the mutual arrangement of the bands is analyzed. The difference between the obtained results and the results of single-band models with the $`s`$ and $`d`$ symmetries of the order parameter is discussed.

The lack of a generally accepted explanation of the pairing mechanism in high-temperature superconductors (HTSC) makes the interpretation of numerous experimental data (including the results of tunnel experiments) difficult. The complexity of the crystal structure of HTSC compounds also complicates the situation. In most cases, it is difficult to explain peculiarities in the tunnel characteristics within the framework of the standard Bardeen-Cooper-Schrieffer (BCS) model. Such peculiarities include the variety of the subgap structure ranging from a linear $`d`$-type structure to a nearly planar structure of the $`s`$ type , the asymmetry of superconducting peaks at the edge of the gap , and nontrivial behavior of the tunnel density of states of $`SIN`$ junctions and of the current-voltage characteristics of $`SIS`$ junctions outside the gap region. In some experiments , strongly suppressed $`{\displaystyle \frac{dI}{dV}}`$ characteristics with very narrow peaks typical of $`s`$-type superconductors were obtained in the subgap region, which contradict at first glance the available data on strong anisotropy of the order parameter in the $`CuO_2`$ plane .

A large number of theoretical models have been proposed to explain the variety of experimental data. These include the models based on various features of the band structure (taking into account, for example, the closeness of the Van Hove singularities to the chemical potential ), on the choice of the order parameter symmetry , and on the properties of the tunnel matrix element  in the single-band BCS scheme. Even the slave-boson approach  was used for explanation the asymmetry of the superconducting peaks and the inhomogeneity of the density of electron states in the normal state. A series of publications  is devoted to the role of inelastic scattering, including the inelastic scattering by antiferromagnetic spin fluctuations. The results of theoretical investigations show that the inclusion of Van Hove singularities and the use of the $`(s+d)`$ symmetry of the order parameter are obviously essential for a correct explanation of the most experimentally observed properties of the tunnel conductivity of HTSC. For this reason, the experimental investigation of the properties of the electron spectrum of high-temperature superconductors in the normal state is very important.

In our opinion, it is the analysis of the crystal structure that makes it possible to explain a number of results of tunnel measurements in high-temperature superconductors. This paper aims at explanation the features of tunnel characteristics using the model with an anisotropic effective order parameter, which is based only on the properties of the electron spectrum of HTSC . The spectrum of the considered model approximately corresponds to the band structure of compounds of the BiSrCaCuO (BSCCO) type. An important feature of this approach is that the initial electron-electron interaction leading to pairing is treated as isotropic in the $`CuO_2`$ plane. The anisotropy of the order parameter and of the excitation spectrum is determined by the symmetry properties of the crystal lattice. According to our calculations, many features of the tunnel characteristics typical of HTSC compounds can be explained by using a rather universal approach which will be described below.

Let us consider the effect of symmetry of the initial bands of the superconductor on the current-voltage characteristics of $`SIS`$-junctions. The electron system of the $`CuO_2`$ plane is described by model Hamiltonian of the form:

$`H`$ $`=`$ $`\epsilon _{z^2}{\displaystyle \underset{i,\sigma }{}}c_{i,\sigma }^+c_{i,\sigma }+\epsilon _d{\displaystyle \underset{i,\sigma }{}}d_{i,\sigma }^+d_{j,\sigma }+\epsilon _p{\displaystyle \underset{i,\sigma }{}}(p_x^+(i,\sigma )p_x(i,\sigma )+p_y^+(i,\sigma )p_y(i,\sigma ))+`$ (1)
$`+`$ $`{\displaystyle \underset{i,j,\sigma }{}}(t_{z^2p_x}^{i,j}c_{i,\sigma }^+p_x(j,\sigma )+h.c.)+{\displaystyle \underset{i,j,\sigma }{}}(t_{z^2p_y}^{i,j}c_{i,\sigma }^+p_y(j,\sigma )+h.c.)+`$
$`+`$ $`{\displaystyle \underset{i,j,\sigma }{}}(t_{dp_x}^{i,j}d_{i,\sigma }^+p_x(j,\sigma )+h.c.)+{\displaystyle \underset{i,j,\sigma }{}}(t_{dp_y}^{i,j}d_{i,\sigma }^+p_y(j,\sigma )+h.c.)+`$
$`+`$ $`U_{z^2}{\displaystyle \underset{i}{}}c_i^+c_ic_i^+c_i+U_d{\displaystyle \underset{i}{}}d_i^+d_id_i^+d_i,`$
Here $`c_{i,\sigma }^+`$ and $`d_{i,\sigma }^+`$ are creation operators for electrons with spin $`\sigma `$ on the $`d_{z^2}`$\- and $`d_{x^2y^2}`$\- orbitals of the $`i`$-th site of copper; $`p_x^+(i,\sigma )`$ and $`p_y^+(i,\sigma )`$ are the creation operators for electrons with spin $`\sigma `$ on the $`p_x`$ and $`p_y`$ orbitals of the $`i`$-th site of oxygen (Fig.1). $`\epsilon _p`$, $`\epsilon _d`$ and $`\epsilon _{z^2}`$ are the energies of the $`p`$ levels of oxygen and of the $`d_{x^2y^2}`$ and $`d_{z^2}`$ levels of copper, respectively, measured from the chemical potential ($`\mu =0`$). $`t_{z^2p_x}^{i,j}=t_{z^2p_y}^{i,j}t_{z^2p}`$ and $`t_{dp_x}^{i,j}=t_{dp_y}^{i,j}t_{dp}`$ are the matrix elements of one-particle transitions between the $`d_{z^2}`$ and $`d_{x^2y^2}`$ orbitals of copper and the $`p`$ orbitals of oxygen. For the sake of simplicity, we assume that superconductivity is caused by the isotropic attraction $`U_{z^2}<0`$ of electrons on the $`d_{z^2}`$ orbital. It should be noted that a more complicated consideration of the superconducting interaction would not lead to a qualitative difference from the case investigated here. We also assume the presence of the isotropic effective electron-electron interaction $`U_d`$ on the $`d_{x^2y^2}`$ orbital of copper.

Introducing new operators $`a_{𝐤,\sigma }^{(\beta )}`$ ($`\beta `$=1,2,3) in the $`𝐤`$ representation according to the formulas

$`d_{𝐤,\sigma }=C_d^{(1)}a_{𝐤,\sigma }^{(1)}+C_d^{(2)}a_{𝐤,\sigma }^{(2)}+C_d^{(3)}a_{𝐤,\sigma }^{(3)}`$
$`p_x(𝐤,\sigma )=C_{p_x}^{(1)}a_{𝐤,\sigma }^{(1)}+C_{p_x}^{(2)}a_{𝐤,\sigma }^{(2)}+C_{p_x}^{(3)}a_{𝐤,\sigma }^{(3)}`$ (2)
$`p_y(𝐤,\sigma )=C_{p_y}^{(1)}a_{𝐤,\sigma }^{(1)}+C_{p_y}^{(2)}a_{𝐤,\sigma }^{(2)}+C_{p_y}^{(3)}a_{𝐤,\sigma }^{(3)},`$
we diagonalize the Hamiltonian part describing one-particle transition between the $`d_{x^2y^2}`$ orbitals of copper and the $`p`$ orbitals of oxygen. The coefficients $`C_d`$ and $`C_{p_x(p_y)}`$ in formulas (2) have the form:

$`C_d^{(1)}=0,`$ $`C_d^{(\alpha )}={\displaystyle \frac{\epsilon _p\epsilon _\alpha }{\sqrt{(\epsilon _p\epsilon _\alpha )^2+(t_x)^2+(t_y)^2}}},`$
$`C_{p_x(p_y)}^{(1)}={\displaystyle \frac{it_{y(x)}}{\sqrt{(t_x)^2+(t_y)^2}}},`$ $`C_{p_x(p_y)}^{(\alpha )}=\pm {\displaystyle \frac{t_{x(y)}}{\sqrt{(\epsilon _p\epsilon _\alpha )^2+(t_x)^2+(t_y)^2}}},`$
where $`\alpha =2,3`$ and $`t_{x(y)}=t_{dp}\mathrm{sin}\frac{k_{x(y)}}{2}.`$ Similarly, the matrix element of a transition between the $`d_{z^2}`$ orbitals of copper and the $`p`$ orbitals of oxygen in the $`𝐤`$ representation can be written as $`\stackrel{~}{t}_{x(y)}=t_{z^2p}\mathrm{sin}\frac{k_{x(y)}}{2}`$. Since the matrix elements $`t_{dp}^{i,j}`$ and $`t_{z^2p}^{i,j}`$ depend on the indices of different sites, the quantities $`t_{x(y)}`$ and $`\stackrel{~}{t}_{x(y)}`$ must be functions of $`𝐤`$. Operators $`a_{𝐤,\sigma }^{(1)}{}_{}{}^{+}`$, $`a_{𝐤,\sigma }^{(2)}{}_{}{}^{+}`$ and $`a_{𝐤,\sigma }^{(3)}{}_{}{}^{+}`$ are the creation operators for electrons with spin $`\sigma `$ and quasimomentum $`𝐤`$ in the bands formed by hybridized $`p`$ orbitals of oxygen and $`d_{x^2y^2}`$ orbitals of copper with the dispersion relation

$`\epsilon _1=\epsilon _p,\epsilon _{2(3)}(𝐤)={\displaystyle \frac{\epsilon _p+\epsilon _d}{2}}{\displaystyle \frac{1}{2}}\sqrt{(\epsilon _p\epsilon _d)^2+4(t_x^2(𝐤)+t_y^2(𝐤))}.`$ (3)
As a result of transformation (2), the system Hamiltonian in the mean-field approximation takes the form:

$`H`$ $`=`$ $`{\displaystyle \underset{𝐤,\sigma }{}}\epsilon _{z^2}c_{𝐤,\sigma }^+c_{𝐤,\sigma }+{\displaystyle \underset{𝐤,\sigma ,\beta }{}}\epsilon _\beta a_{𝐤,\sigma }^{(\beta )}{}_{}{}^{+}a_{𝐤,\sigma }^{(\beta )}`$ (4)
$`+`$ $`{\displaystyle \underset{𝐤,\sigma ,\beta }{}}(W_\beta (𝐤)c_{𝐤,\sigma }^+a_{𝐤,\sigma }^{(\beta )}+h.c.)`$
$``$ $`{\displaystyle \underset{𝐤}{}}(\mathrm{\Delta }_{z^2}c_{𝐤,}^+c_{𝐤,}^++h.c){\displaystyle \underset{𝐤,\beta ,\gamma }{}}(\mathrm{\Delta }_{\beta \gamma }a_{𝐤,}^{(\beta )}{}_{}{}^{+}a_{𝐤,}^{(\gamma )}{}_{}{}^{+}+h.c),`$
where

$`W_1(𝐤)={\displaystyle \frac{\stackrel{~}{t}_xt_y+\stackrel{~}{t}_yt_x}{\sqrt{t_x^2+t_y^2}}},W_\alpha (𝐤)={\displaystyle \frac{t_y\stackrel{~}{t}_yt_x\stackrel{~}{t}_x}{\sqrt{(\epsilon _p\epsilon _\alpha )^2+(t_x)^2+(t_y)^2}}},\alpha =2,3`$ (5)
and new order parameters

$`\mathrm{\Delta }_{\beta \gamma }(𝐤)=C_d^{(\beta )}(𝐤)C_d^{(\gamma )}(𝐤)\mathrm{\Delta }_d,`$ $`\beta ,\gamma =1,2,3,`$ (6)
have been introduced, in which

$`\mathrm{\Delta }_d={\displaystyle \frac{1}{N}}U_d{\displaystyle \underset{𝐤}{}}d_{𝐤,}d_{𝐤,},`$ $`\mathrm{\Delta }_{z^2}={\displaystyle \frac{1}{N}}U_{z^2}{\displaystyle \underset{𝐤}{}}c_{𝐤,}c_{𝐤,}.`$ (7)
The arrangement of the energy bands corresponding to Hamiltonian (4) is shown in the inset to Fig.1. In the framework of a realistic description of BSCCO-type compounds , it is assumed that the chemical potential $`\mu `$ and the $`z^2`$ level of copper locate closely to the center of the upper band $`\epsilon _3(𝐤)`$. The remaining parameters of the model are assumed to satisfy the relations $`\epsilon _d<\epsilon _p`$, $`\epsilon _p\epsilon _dt_{dp}`$, $`t_{z^2p}|\epsilon _d|,|\epsilon _p|`$, $`|\mathrm{\Delta }_d|<\mathrm{\Delta }_{z^2}|\epsilon _d|,|\epsilon _p|`$.

Since only the energy range near the chemical potential (on the order of several $`\mathrm{\Delta }_{z^2}`$) plays an important role here, we can disregard the two lower bands $`\epsilon _1`$ and $`\epsilon _2(𝐤)`$ taking into account the relation between parameters and restrict our analysis to the upper band $`\epsilon _3(𝐤)`$ only. In addition to one-particle transitions included in Hamiltonian (1), we can take into account the additional hybridization of the $`d_{z^2}`$ orbital with the atoms surrounding the $`CuO_2`$ plane, including $`BiO`$ and $`SrO`$ complexes of the BSSCO-type compounds as well as direct transition between the $`d_{z^2}`$ orbitals of different copper atoms. In this connection, we assume that the initial $`z^2`$ band is characterized by the dispersion relation $`\epsilon _{z^2}(𝐤)=\epsilon _{z^2}+t_{z^2}(cos(k_x)+cos(k_y))`$ with a width much smaller than for $`\epsilon _2(𝐤)`$ and $`\epsilon _3(𝐤)`$ bands. Thus, in the chosen approximation, the initial problem can be reduced to a two-band model in which the initial $`\epsilon _{z^2}(𝐤)`$ and the $`\epsilon _3(𝐤)`$ bands with the one-particle hybridization $`W_3(𝐤)`$ and the order parameters $`\mathrm{\Delta }_{z^2}`$ and $`\mathrm{\Delta }_{33}(𝐤)`$ are considered.

We introduce the following time Green’s functions:

$`G_{z^2}(𝐤;t,t^{})=iTc_{𝐤,\sigma }(t)c_{𝐤,\sigma }^+(t^{}),`$ $`g_{\sigma \sigma ^{}}F_{z^2}^+(𝐤;t,t^{})=iTc_{𝐤,\sigma }^+(t)c_{𝐤,\sigma ^{}}^+(t^{}),`$
$`G_{z^2,3}(𝐤;t,t^{})=iTc_{𝐤,\sigma }(t)a_{𝐤,\sigma }^{(3)+}(t^{}),`$ $`g_{\sigma \sigma ^{}}F_{z^2,3}^+(𝐤;t,t^{})=iTa_{𝐤,\sigma }^{(3)+}(t)c_{𝐤,\sigma ^{}}^+(t^{}),`$
where $`\widehat{g}=i\widehat{\sigma }^y`$. Using the equations of motion for operators $`a_{𝐤,\sigma }^{(3)}`$ and $`c_{𝐤,\sigma }`$ and transforming to the frequency representation, we obtain the following set of equations for the $`z^2`$ Green’s functions:

$$\left(\omega \epsilon _{z^2}(𝐤)\right)G_{z^2}(𝐤,\omega )W_3(𝐤)G_{z^2,3}^+(𝐤,\omega )+\mathrm{\Delta }_{z^2}F_{z^2}^+(𝐤,\omega )=1,$$
(8)
$$\left(\omega \epsilon _3(𝐤)\right)G_{z^2,3}^+(𝐤,\omega )W_3(𝐤)G_{z^2}(𝐤,\omega )+\mathrm{\Delta }_3(𝐤)F_{z^2,3}^+(𝐤,\omega )=0,$$
(9)
$$\left(\omega +\epsilon _{z^2}(𝐤)\right)F_{z^2}^+(𝐤,\omega )+W_3(𝐤)F_{z^2,3}^+(𝐤,\omega )+\mathrm{\Delta }_{z^2}^+G_{z^2}(𝐤,\omega )=0,$$
(10)
$$\left(\omega +\epsilon _3(𝐤)\right)F_{z^2,3}^+(𝐤,\omega )+W_3(𝐤)F_{z^2}^+(𝐤,\omega )+\mathrm{\Delta }_3^+(𝐤)G_{z^2,3}^+(𝐤,\omega )=0,$$
(11)
where

$$\mathrm{\Delta }_3(𝐤)\mathrm{\Delta }_{33}(𝐤)=\mathrm{\Delta }_d\frac{(\epsilon _p\epsilon _3)^2}{(\epsilon _p\epsilon _3)^2+(W_x)^2+(W_y)^2}.$$
(12)
In the following, we will make use of quasiparticle density of states of the $`z^2`$ band,

$`N(\omega )={\displaystyle \frac{1}{\pi }}{\displaystyle ImG_{z^2}^R(\omega ,𝐤)\frac{d^3k}{(2\pi )^3}}.`$ (13)
The retarded Green’s function is obtained by the solution of the system of equations (8-11) and has the form:

$$G_{z^2}^R(\omega ,𝐤)=\frac{u_{,𝐤}^2}{\omega E_{}(𝐤)+i0}+\frac{u_{+,𝐤}^2}{\omega E_+(𝐤)+i0}+\frac{v_{,𝐤}^2}{\omega +E_{}(𝐤)+i0}+\frac{v_{+,𝐤}^2}{\omega +E_+(𝐤)+i0},$$

where the coherence factors are defined as

$$u_{\pm ,𝐤}^2=\frac{(E_{}+\epsilon _{z^2})(E_{}^2+\epsilon _3^2+\mathrm{\Delta }_3^2)+W_3^2(E_{}\epsilon _3)}{2E_{}(E_+^2E_{}^2)}$$

$$v_{\pm ,𝐤}^2=\frac{(E_{}\epsilon _{z^2})(E_{}^2+\epsilon _3^2+\mathrm{\Delta }_3^2)+W_3^2(E_{}+\epsilon _3)}{2E_{}(E_+^2E_{}^2)}.$$

The dispersion relations for two branches of the excitation spectrum have the form

$`E_\pm ^2(𝐤)`$ $`=`$ $`{\displaystyle \frac{\epsilon _{z^2}^2(𝐤)+\mathrm{\Delta }_{z^2}^2+\epsilon _3^2(𝐤)+\mathrm{\Delta }_3^2(𝐤)+2W_3^2(𝐤)}{2}}`$
$`\pm `$ $`{\displaystyle \frac{\sqrt{(\epsilon _{z^2}^2(𝐤)+\mathrm{\Delta }_{z^2}^2\epsilon _3^2(𝐤)\mathrm{\Delta }_3^2(𝐤))^2+4W_3^2(𝐤)\{(\mathrm{\Delta }_{z^2}\mathrm{\Delta }_3(𝐤))^2+(\epsilon _{z^2}(𝐤)+\epsilon _3(𝐤))^2\}}}{2}}.`$
The set of equations (8-11) can be reduced to the system of two equations for Green’s functions $`G_{z^2}`$ and $`F_{z^2}`$ with the effective order parameter defined by the formula

$$\mathrm{\Delta }_{z^2}(𝐤)=\frac{\mathrm{\Delta }_3(𝐤)W_3^2(𝐤)}{\epsilon _3^2(𝐤)+\mathrm{\Delta }_3^2(𝐤)}+\mathrm{\Delta }_{z^2}.$$
(15)
The sign of the order parameter $`\mathrm{\Delta }_3(𝐤)`$(12) is determined by the sign of the parameter $`\mathrm{\Delta }_d`$(7), which can be positive or negative depending on the type of interaction $`U_d`$ (repulsion or attraction). The parameter $`\mathrm{\Delta }_d=0`$ if $`U_d=0`$ or $`U_d>U_d^{crit}>0`$, where $`U_d^{crit}`$ is a certain critical value of repulsion on the $`d_{x^2y^2}`$ orbital of copper which suppresses superconductivity in the system. As seen from the formula (15), in the case of a nonzero interaction $`U_d`$ on the $`d_{x^2y^2}`$ orbital ($`\mathrm{\Delta }_d0`$), the order parameter depends on quasimomentum, and in the case of repulsion, $`\frac{\mathrm{\Delta }_3(𝐤)}{\mathrm{\Delta }_{z^2}}<0`$ ($`U_d>0`$), the parameter $`\mathrm{\Delta }_{z^2}(𝐤)`$ changes its sign.

The approach considered explains the strong anisotropy of the order parameter as a consequence of the symmetry properties of the one-particle matrix element of the interband hybridization $`W_3(𝐤)`$ , which has nodes along the diagonals of the Brillouin zone due to the difference between the types of symmetry of the initial bands. In this case, the branch $`E_{}(𝐤)`$ of the excitation spectrum vanishes at points located on lines in the $`𝐤`$ space  at which the effective order parameter $`\mathrm{\Delta }_{z^2}(𝐤)`$ is equal to zero. It should be noted that the initial interaction is regarded as isotropic; i.e., the effect considered does not depend on the origin of the pairing mechanism. In addition, the approach does not require a strong anisotropy of the spectrum and, hence, can be used for various types of high-temperature superconductors. The model proposed in  also leads to formulas of the type (14), (15), but it is based on the exotic condition of the interaction sign reversal in various regions of the Fermi surface.

The diagonalization of the one-particle part of Hamiltonian (1) leads to the problem with anisotropic attraction in energy bands. The effective order parameters in the bands may have nodes, but have no pure $`d`$ or $`(s+d)`$ symmetry since their anisotropy is determined by the band representations of the space symmetry group of the lattice .

The formulated model of the order parameter anisotropy makes it possible to obtain the $`{\displaystyle \frac{dI}{dV}}`$ characteristics of $`SIS`$ junctions for superconductors of the BSCCO type. Taking into account the crystal structure of these compounds, we will assume that tunneling along the $`c`$ axis occurs mainly through the $`d_{z^2}`$ orbitals of copper in the $`CuO_2`$ plane (and through the apical oxygen which is not included explicitly in the model). The matrix element of tunneling is assumed to be independent of momentum ($`T_{\mathrm{𝐤𝐩}}=T=const`$) due to the random formation of bonds between the $`d_{z^2}`$ orbitals on both sides of the break junction. Thus, tunneling between two superconductors ($`CuO_2`$ layers) is treated as occuring through a number of point junctions.

The expression for the dependence of the quasiparticle tunnel current on the voltage applied to the junction in this case assumes the standard form:

$`I(V)=4e|T|^2{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}[n(\omega )n(\omega eV)]N(\omega )N(\omega eV)𝑑\omega ,`$ (16)
where $`n(\omega )=[e^{{\scriptscriptstyle \frac{\omega }{T}}}+1]^1`$ is the Fermi distribution function and $`N(\omega )`$ is defined by formula (13). All calculations were performed for temperature $`T=0`$. An important aspect of this approach is that $`N(\omega )`$ is not the average density of states in the conduction band (as in the single-band BCS model), but is a partial density on the $`d_{z^2}`$ orbital.

In order to find the density of states $`N(\omega )`$ and the tunnel conductance $`{\displaystyle \frac{dI}{dV}}`$ in accordance with the chosen approximation, we used the values of parameters from the tight-binding model, which correspond to the calculations of the band structure of HTSC . If all parameters are expressed in terms of $`\mathrm{\Delta }_{z^2}`$ and the position of energy levels is measured from the chemical potential, we have $`t_{dp}=75`$, $`t_{z^2p}=15`$, $`\epsilon _d=75`$ and $`\epsilon _p=50`$. In all calculations, we assumed the existence of a finite relaxation constant $`\gamma =0.05`$. The parameter determining the width of the initial $`\epsilon _{z^2}(𝐤)`$ band is $`t_{z^2}=1.5`$.

For the chosen values of the parameters, we analyzed the behavior of the characteristics $`N(\omega )`$ and $`{\displaystyle \frac{dI}{dV}}`$ depending on the interaction value between electrons on the $`d_{x^2y^2}`$ orbital of copper and on the position of the center $`\epsilon _{z^2}`$ of the initial $`z^2`$ band relative to the chemical potential (Figs. 2-4). The obtained dependencies were compared with these of single-band model with the $`s`$\- and $`d`$-types of symmetry of the order parameter ($`\mathrm{\Delta }=const`$ and $`\mathrm{\Delta }(cos(k_x)cos(k_y))`$, respectively).

Figure 2A shows the $`N(\omega )`$ curves calculated for the normal state ($`\mathrm{\Delta }_{z^2}=0`$). These curves display split peaks associated with the Van Hove singularities of the initial $`\epsilon _{z^2}(𝐤)`$ band for various positions of its center relative to the chemical potential, as well as the following peaks formed by singularities of the $`\epsilon _3(𝐤)`$ band. Figure 2b shows the same curves in the case when the $`\epsilon _{z^2}(𝐤)`$ spectrum is replaced by a dispersionless level ($`t_{z^2}=0`$).

Fig.3.A presents the $`N(\omega )`$ dependencies in the superconducting state for the case when the $`z^2`$ band is far from the chemical potential ($`\epsilon _{z^2}=10`$). It can be seen from the figure that if there is electron-electron repulsion on the $`d_{x^2y^2}`$ orbital $`\mathrm{\Delta }_d=0.5`$), the density of quasiparticles is similar to that calculated in a model with the $`d`$ symmetry of the order parameter (in particular, it is a linear function of $`\omega `$ in the vicinity of $`\omega =0`$). If there is no interaction on the $`d_{x^2y^2}`$ orbital, the distance between the superconducting peaks and their height increase. According to analytical estimations, $`N(\omega )\omega ^{3/2}`$ for very small frequencies. Similar results can also be obtained in the single-band model with a nontrivial spectrum and an effective order parameter of certain symmetry. As noted above, taking into account the above approximations, the diagonalization of the single-particle part of Hamiltonian (1) leads to the single-band model with anisotropic pairing. In this case, the symmetry of the corresponding order parameter can be approximated by an $`(s+d)`$-type symmetry.

If the center of the initial $`z^2`$ band is close to the chemical potential (see Figs. 3B and 3C), the difference between the $`a`$ and $`b`$ curves (with and without taking into account of repulsion) vanishes due to the dominant role of the Van Hove singularities on energy scales of the order of $`\mathrm{\Delta }_{z^2}`$. In both cases, the behavior of the density of quasiparticles becomes of the $`s`$ type. The calculated density of states demonstrates the asymmetry of peaks, which is associated with the position of Van Hove singularities in the normal density of states (see Fig. 2). The Van Hove singularity nearest to the chemical potential increases the height of the corresponding peak in the density of states in the superconducting state. A comparison of Figs. 3B and 3C shows that, for a certain intermediate value of $`\epsilon _{z^2}`$, there is the mirror switching of the asymmetry of peaks associated with the model spectrum structure.

The curves presented in Fig. 3A show that even without the repulsion on the $`d_{x^2y^2}`$\- orbital if the $`\epsilon _{z^2}`$ level is far from the chemical potential ($`\epsilon _{z^2}=10`$), the distance between the superconducting peaks is smaller than $`2\mathrm{\Delta }_{z^2}`$. As the $`\epsilon _{z^2}`$ level approaches the chemical potential, the distance between the superconducting peaks increases, and the position of the peaks on the $`N(\omega )`$ curve for $`\epsilon _{z^2}=0`$ corresponds to the parameter $`\mathrm{\Delta }_{z^2}`$. In this case, the height of the peaks increases significantly. The distance between the peaks changes as a result of the displacement of the $`\epsilon _{z^2}`$ level relative to the chemical potential for the fixed values of the remaining parameters. However, a more detailed analysis of this phenomenon requires the solution of the self-consistent equations for order parameters as functions of the parameters of the model as well as the determination of the interaction constants $`U_{z^2}`$ and $`U_d`$ taking into account the specific mechanism of pairing, which is beyond the scope of the present paper. The dependence of the distance between the superconducting peaks on the position of the chemical potential in the model with the $`(s+d)`$-type of order parameter symmetry was also found theoretically in .

A comparison of the $`N(\omega )`$ dependencies found in the model with single-band models with the $`s`$\- and $`d`$-types of the order parameter symmetry was carried out taking into account the location of the Van Hove singularities associated with the initial $`\epsilon _{z^2}(𝐤)`$ band. This is pointed out by similar asymptotic behavior of densities of states calculated for both the normal and the superconducting state in $`s`$-, $`d`$\- and considered models (see Fig. 3C) on the scales of energy higher than $`\mathrm{\Delta }_{z^2}`$.

A comparison of the results presented in Figs. 3B and 3D shows that, as the width of the initial $`z^2`$ band determined by the parameter $`t_{z^2}`$ decreases ($`t_{z^2}=1.5`$ in Fig.3B and $`t_{z^2}=0`$ in Fig.3D), the heights of the superconducting peaks increases, since the singularity in the normal density of states in this case is displaced towards the chemical potential (see Fig. 2).

Figure 4 shows the $`{\displaystyle \frac{dI}{dV}}`$ characteristics of $`SIS`$ junction calculated by formula (16) using the determined $`N(\omega )`$ dependencies for three positions of the $`\epsilon _{z^2}`$ level relative to the chemical potential $`\epsilon _{z^2}=10`$, $`\epsilon _{z^2}=1.5`$, $`\epsilon _{z^2}=0`$, which can correspond to different HTSC types. The common feature of these characteristics is that in the case of repulsion the superconducting peak appears at lower voltages and is lowered, but the curve at low voltages lies higher than in the absence of interaction on the $`d_{x^2y^2}`$ orbital of copper. If the center of the initial $`z^2`$ band is far from the chemical potential ($`\epsilon _{z^2}=10`$), the $`{\displaystyle \frac{dI}{dV}}`$ characteristics are similar to those calculated in a model with the $`d`$ symmetry of the order parameter (Fig. 4a). As the $`\epsilon _{z^2}`$ level approaches the chemical potential, the curves tend to the dependencies obtained in the model with the $`s`$ symmetry of the order parameter (Figs. 4b and 4c). For low voltages $`V<\mathrm{\Delta }_{z^2}`$, the $`IV`$ characteristic of an $`SIS`$ junction is strongly suppressed. In analogy with the $`N(\omega )`$ dependencies, the $`{\displaystyle \frac{dI}{dV}}`$ characteristics exhibit extremely narrow superconducting peaks. This corresponds to a number of experimentally observed results for BSCCO-type compounds, e.g., optimally doped , overdoped , and underdoped  $`Bi_2Sr_2CaCu_2O_{8+\delta }`$.

Thus, the electron density of states and current-voltage characteristics of $`SIS`$ junctions are calculated in the framework of the multiband HTSC model with an anisotropic effective order parameter. We took into account the hybridization between the $`d_x`$ and $`d_y`$ orbitals of oxygen and the $`d_{x^2y^2}`$ orbitals of copper in the $`CuO_2`$ plane and assumed the presence of electron-electron attraction on the $`d_{z^2}`$ orbital of copper. The performed calculations demonstrate the possibility of explanation the results of experimental measurements of the electron density of states and the tunnel properties of $`SIN`$ and $`SIS`$ junctions (e.g., $`s`$-type characteristics of BSCCO-type compounds) on the basis of the model taking into account the structure of the electron spectrum of HTSC. It should be stressed that anisotropy of the order parameter and of the excitation spectrum is explained only by the symmetry of the crystal lattice and atomic orbitals and does not depend on the nature of pairing.

The main result of this work is that for a strongly anisotropic order parameter $`\mathrm{\Delta }(𝐤)`$ (15) (e.g. of $`d`$-type) and in the presence of nodes in the excitation spectrum of a superconductor, the simple inclusion of the real band structure of the HTSC in a wide range of the model parameter values results in the experimentally observed  $`s`$-type behavior of the current-voltage characteristics $`SIN`$ and $`SIS`$ junctions. A comparison of the curves calculated for various values of the parameters of the model and in the simplest cases of models with the $`s`$ and $`d`$ symmetries of order parameters indicates a strong dependence of both the density of states and the $`{\displaystyle \frac{dI}{dV}}`$ characteristics on the electron band structure of high-temperature superconductors. These characteristics also depend on the crystal lattice symmetry and on the presence of an additional interaction between electrons, which was introduced in the model on the $`d_{x^2y^2}`$ orbital of copper. Above analysis points to the necessity of considering of the band structure of HTSC for interpretation of the results of tunnel experiments .

The authors are grateful to B.A. Volkov, E.G. Maksimov, and Ya.G. Ponomarev for numerous discussions in the course of the research. This work was supported by the Russian Foundation for Basic Research (N 02-02-16925, N 00-15-15998) and Landau Scholarship Grant.