text:
# Selected aspects of the classical molecular dynamics of 𝑁 bodies in strong interaction.

## I Description of the code.

Let us start by describing the Classical N-Body Dynamics code (labeled CNBD) used in this article. The basic ingredients of such a code are very simple. The dynamical evolution of each body of the system is driven by the classical Newtonian equations of motion. The two-body potential used in the present work is a third degree polynomial whose derivatives are null at the range $`r_1`$ and at the distance of maximum depth $`r_0`$. The depth value is $`V_{min}`$ and the value at $`r=0`$ is finite and equal to $`V_0`$. The shape of the two-body potential is shown on figure 1. This potential has the basic properties of the Lennard-Jones potential used in other works PlagnolBormio ; Dorso99 : a finite range attractive part and a repulsive short range part. To follow the dynamical evolution of the system an adaptative stepsize fourth-order Runge-Kutta algorithm is used NumRec . The main difference with other works is that the time step $`\mathrm{\Delta }t`$ can vary: if the potential varies strongly, $`\mathrm{\Delta }t`$ is small and when the potential varies gently, $`\mathrm{\Delta }t`$ becomes larger. This allows a very high accuracy with shorter CPU time than for fixed time step algorithms. It requires an additional simulation parameter $`ϵ`$ which is adjusted to ensure the verification of conservation laws (energy, momentum, angular momentum) with a reasonable simulation time. For an $`ϵ`$ value of $`10^5`$, typical CPU times for a collision of two clusters with 50 particles each with an ending time equal to 200 time simulations units on a Compaq DS20 computer under the UNIX True64 operating system ranges from $`30`$ to $`400`$ seconds, depending mainly on the impact parameter. The energy gap between the beginning and the the ending simulation time is lower than 0.001%. This simulation has five free parameters: four linked to the physics (the interaction) and one linked to the numerical algorithm ($`ϵ`$).

Since one wants to study the simplest case, neither long range repulsive interaction nor quantum corrections like a Pauli potential have been introduced Dorso87 . Additionally, no statistical decay code is applied on the excited fragments formed during the collision. The final products have to be regarded as “primary” products which will decay afterwards.

In order to avoid any confusion with nuclear physics, the units used here are called Simulation Units and noted $`S.U.`$. The distance will then be in Distance Simulation Units ($`D.S.U.`$), the energies in Energy Simulation Units ($`E.S.U.`$), the velocities in Velocity Simulation Units ($`V.S.U.`$) and the reaction time in Time Simulation Unit ($`T.S.U.`$). We will be only interested in the relative evolutions of the observables and in their link to the properties of the stable systems. The main goal of the present work is not to reproduce the experimental data of nucleus-nucleus collisions, but rather to see to which extent this simple simulation is qualitatively similar, or not, to experimental data.

## II Static properties of “ground states”.

Once the basic ingredients are defined, one can build stable systems. Since the two-body potential only depends on the distance between the two bodies, such systems are small crystals. The “ground states” of such systems are defined as the configuration in position space which minimizes their total energy. This is obtained by using a Metropolis simulated annealing method NumRec . The locations of particles obtained this way are very close to those obtained by using a basin-hopping algorithm for Lennard-Jones clusters WALES97 .

On top-left of figure 2 is displayed the energy per body $`E_{Bind}/N`$ of these “ground states” as a function of the number $`N`$ of bodies in the cluster (solid line). The dashed and dotted line corresponds to a fit using a liquid-drop formula. The dashed line corresponds to the energy of the less bound body $`E_{LessBound}`$ in the cluster. It will be seen that this energy seems to play a particular role in cluster-cluster collisions. The overall dependence of $`E_{Bind}/N`$ with $`N`$ is very similar to what is seen in nuclear physics. The main difference is seen for high values of $`N`$ because of the absence of a Coulomb-like interaction: $`E_{Bind}/N`$ continues to decrease with $`N`$ whereas it increases for nuclei. These values are close to those found by using the basin-hopping algorithm WALES97 .

The dependence with $`N`$ of the root mean radius radius $`r(N)=\sqrt{<r^2>}`$ of N-body clusters is shown on the top-right panel of figure 2 (solid line). The dashed and dotted line corresponds to a fit using a $`N^{\frac{1}{3}}`$ term. Here again, this is similar to what is known for nuclei. The main difference is the $`r_0`$ term which is due to the large size of the “repulsive core” compared to the range of the attractive part of the interaction.

Since a size can be defined, an effective density can be calculated. One can stretch and squeeze the N-body cluster and build an effective “zero temperature equation of state”. This is displayed on the low-left panel on figure 2 for various system sizes. The density $`\rho _{eff}`$ on this abscissa is the density relative to the density for the “ground state”. The curvature of this “equation of state” can be computed and defined as the “compressibility modulus” $`K(N)`$ of the system. Its variations with $`N`$ are displayed on the low-right panel of figure 2. One can see on these two plots that the “equation of state” and $`K(N)`$ are strongly dependent on $`N`$ for values of $`N`$ below 20$``$30 and then less dependent on $`N`$ above. One could be able to define an “equation of state of infinite matter” by computing the limit of these evolutions for huge values of $`N`$. Here again, the behavior of these classical N-body systems is very similar to that of nuclei.

This first study shows that the classical N-body clusters have strong similarities with atomic nuclei. This suggests that all the parametrisation used to describe nuclei (liquid drop parametrisation, radii, equation of state, etc…) could be directly deduced from the parameters of the two-body interaction. Since this interaction remains unchanged, the differences seen for different clusters has to be attributed only to $`N`$ and to the geometrical configuration of the bodies in the clusters.

## III Cluster-cluster collisions.

After studying the static properties of stable N-body systems, collisions between clusters can be investigated to see what kind of reaction can be obtained. Roughly 20,000 collisions have been generated for different system sizes, different entrance channel asymmetries and different impact parameters. For each collision, the orientation of the inertia axis of clusters are randomly chosen. This avoid to run twice the same collision if the impact parameter, the target and the projectile size and the incident energy are the same for two different collisions. The simulations have been performed for the following systems: 13 + 13, 34 + 34, 50 + 50, 100 + 100 and 18 + 50. The incident energies have been chosen in a way that the available energy in the center of mass frame ranges from an energy far below the binding energy of the fused system ($`30E.S.U`$) to an energy well above the binding energy of the fused system ($`120E.S.U.`$). As it will be shown, these collisions can be studied in the same way as the experimental data of nucleus-nucleus collisions are. These collisions can be seen as “numerical experiments”. Only few examples of such analyses will be shown. More detailed analyses will be done in forthcoming articles. Firstly, a small list of reaction mechanisms obtained will be done. Secondly, the excitation energy stored in excited clusters will be studied and a link to the static properties will be done. Finally, the question of the statistical description of this dynamical model will be briefly addressed.

### III.1 Reaction mechanisms.

On figure 3 are displayed the pictures of cluster-cluster collisions obtained at a fixed time. The ending time of these collisions is $`t=200T.S.U.`$. On each picture, the bodies which were originally belonging to the projectile are in dark grey and those which were belonging to the target are in light grey. The most striking observation is that all kinds of reaction mechanisms observed in nucleus-nucleus collisions seem to be present in this very simple simulation. One can see low energy processes like the fusion/evaporation process (first row left panel), pure binary collision (second row left panel), stripping/pick-up mechanism (second row right panel) and deep inelastic collisions (third row right panel). Intermediate energy processes, like neck formation and break-up (first row right panel) and multi-fragmentation (third row left panel), are seen. Finally, high energy processes like the participant/spectator process (lower-most panel) are also seen.

This similarity can also be seen on the so-called Wilczynski plots Wilczynski  shown on figure 4. These plots display the correlation between the flow angle $`\theta _{flow}`$ and the total kinetic energy of the clusters $`E_{kin}=E_k^i`$ where $`E_k^i`$ is the kinetic energy of cluster $`i`$ in the center of mass frame. Two bodies are assumed to belong to the same cluster if they are in interaction, i.e. if their relative distance is below the range $`r_1`$ of the interaction (see figure 1). This algorithm of cluster recognition is the simplest one and is called Minimum Spanning Tree algorithm in other works Dorso99 ; Dorso93 . The most dissipative collisions correspond to the smallest $`E_{kin}`$ values: the available energy is converted in internal energy of clusters. These plots were built for a $`N=34`$ projectile colliding a $`N=34`$ target at four available energies $`E_{cm}/N`$ in the center of mass: $`E_{cm}/N=60E.S.U.`$, $`E_{cm}/N=90E.S.U.`$, $`E_{cm}/N=120E.S.U.`$ and $`E_{cm}/N=30E.S.U.`$ which corresponds respectively to the energy of the less bound body for the fused $`N=68`$ cluster, to the binding energy per body for $`N=68`$, to the energy of the most bound body for $`N=68`$ and to an energy below the energy of the less bound body for $`N=68`$. For each energy, 1,000 collisions have been computed assuming a flat impact parameter distribution ranging from $`b=0D.S.U.`$ to the sum of the two cluster radii plus the range of the interaction ($`b_{max}=r(N_{proj})+r(N_{targ})+r_1`$, where $`r_1`$ is the range of the two-body interaction). In the analyses, each collision is weighted assuming a triangular impact parameter distribution between 0 and $`b_{max}`$ ($`weightb`$).

At $`E_{cm}/N=30E.S.U.`$ (upper left panel), the flow angle is always negative. This means that the projectile and the target like clusters are deflected to the opposite direction relative to their original one. At this energy, the attractive part of the interaction is dominant. For the lowest $`E_{kin}`$ values, all possible values of $`\theta _{flow}`$ are covered: this corresponds to the fusion/evaporation process. At $`E_{cm}/N=60E.S.U.`$ (upper right panel), the picture is changed and the range of $`\theta _{flow}`$ values is smaller than for the previous energy: the repulsive part starts to act. The fusion/evaporation area (small $`E_{kin}`$ values, all $`\theta _{flow}`$ values) is still present. For the two highest energies (lower panels), the picture is roughly the same. For the less dissipative collisions, the attractive part is still dominant (negative $`\theta _{flow}`$ values). But when the dissipation increases, the repulsive part becomes dominant and the projectile and the target bounce on each other (positive $`\theta _{flow}`$ values). For the most dissipative collisions, only positive $`\theta _{flow}`$ values are seen which indicates the disappearance of the fusion/evaporation process. These evolutions are qualitatively similar to what is seen in nucleus-nucleus collisions (see for example Metivier2000 ).

These two studies indicate that the description of the reaction mechanisms in terms of mean-field at low energies and in terms of nucleon-nucleon collision at high energies can be deduced from the properties of the nucleon-nucleon interaction only. The similarities of this classical N-body simulation with nucleus-nucleus collisions suggest that, apart for the quantum mechanics effects, an effective description of nucleus-nucleus collisions could be obtained with a reduced number of parameters. Providing the interaction has a finite range attractive part and a short range repulsive part, the overall behavior of the N-body systems seems to be independent of the values of the parameters of this interaction. For a more quantitative agreement with nucleus-nucleus collisions, additional physical ingredients (quantum mechanics, Coulomb interaction) have to be included.

### III.2 Energy deposit in clusters.

Let us now focus our study of cluster-cluster collision on selected topics. Of particular interest is the excitation energy stored in clusters. How this energy can be linked to the available energy and how it is linked to the properties of the “ground state” characteristics of these clusters?

One can plot for example the evolution of the excitation energy $`E^{}/N`$ of the cluster with its parallel velocity $`V_{//}`$ for $`N=34`$ \+ $`N=34`$ collisions (figure 5) for the whole impact parameter range. The excitation energy of each cluster is simply the difference between the total energy (potential plus kinetic) and the “ground state” energy of the cluster. It reads:

$$E^{}=\underset{i=1}{\overset{N}{}}E_{kin}^i+\underset{i,j>i}{}V(r_{ij})E_{Bind}(N)$$
(1)
where $`E_{kin}^i`$ is the kinetic energy of the body $`i`$ in the cluster center of mass frame, $`r_{ij}`$ the relative distance between the bodies $`i`$ and $`j`$, $`V(r_{ij})`$ the potential energy and $`E_{Bind}(N)`$ the binding energy of the “ground-state” of the cluster with $`N`$ bodies. This excitation energy is determined at the end of the calculation corresponding to $`t=200T.S.U.`$. This energy is very close to the one obtained at the separation time of the clusters (the smallest time at which clusters can be identified), since in this time range the evaporation is very weak and the clusters have no time to cool down significantly Dorso99 . On each panel of figure 5, the full line corresponds to the expected correlation between $`E^{}/N`$ and $`V_{//}`$ for a pure binary scenario (the excitation energy is only due to the velocity damping of each partner), the horizontal dashed line to the energy of the less bound body for the fused system $`N=68`$ and the small circle is centered around the expected values of velocity and excitation energy for the fused system.

At $`E_{cm}/N=30E.S.U.`$, the points are slightly below the full line. This means that the excitation energy is strongly linked to the velocity damping. The small shift is due to mass transfers between the projectile and the target, and to promptly emitted clusters. The area corresponding to the fused system is well populated showing that a complete fusion process occurs. At $`E_{cm}/N=60E.S.U.`$, the distribution of points is roughly compatible with the pure binary process hypothesis, but the complete fusion process area is empty. The horizontal line seems to be an upper limit to the excitation energy which can be stored in clusters. Clusters with rather small excitation energies are found at velocities around the center of mass velocity. For the two highest energies, this trend is enhanced: a set of clusters are around the full line, and when the full line is above the horizontal line, one finds clusters at small excitation energies around the center of mass velocity. The energy of the less bound body seems to be a limit to the excitation energy which can be stored in these clusters.

This can be more clearly seen when the excitation energy $`E^{}/N`$ is plotted as a function of $`N`$, as in figure 6. On each panel, the full line corresponds to the energy of the less bound body $`E_{LessBound}`$ in the cluster and the dashed line to the binding energy per body $`E_{bind}/N`$. As in figure 5, the small circle corresponds to the expected values for the fused system. At $`E_{cm}/N=30E.S.U.`$, the area corresponding to complete fusion is filled and all the available energy can be stored as excitation energy. But for higher energies, one can clearly see that for each fragment size, $`E^{}/N`$ never overcomes $`E_{LessBound}`$. At $`E_{cm}/N=60E.S.U.`$ clusters with sizes higher than the projectile size and the target size can be seen. This area corresponds to an incomplete fusion process. For the two highest energies ($`E_{cm}/N=90E.S.U.`$ and $`E_{cm}/N=120E.S.U.`$), the plots are almost identical: there is no more fusion and the clusters are smaller than the target and than the projectile. One can notice that $`E^{}/N`$ never reaches the binding energy $`E_{bind}/N`$ except for small clusters where $`E_{bind}/N`$ and $`E_{LessBound}`$ are equal.

This limitation of excitation energy can be understood quite easily. The less bound body remains bound to the cluster only if its total energy is negative, i.e. its kinetic energy due to the excitation is below its potential one. If one assumes that the excitation energy is roughly equally shared over all bodies in the cluster, when the kinetic energy balances the potential energy of the less bound body, this body is no more bound to the cluster and can escape. To be observed for a long time, the excited cluster must have an excitation energy per body below the energy of the less bound body.

The mechanism of energy deposit in classical N-body clusters seems to be the following one: the excitation seems to be mainly driven by the velocity damping of the two partners and to a lesser extent by exchanges of bodies between them. Once the energy of the less bound body is reached, unbound bodies and/or clusters escape quickly and keep an excitation energy per body below the energy of the less bound one. As a consequence, the highest energy deposit per body can only be obtained at energies close to $`E_{LessBound}`$. For higher available energies, the system fragments quickly, leaving rather “cold” clusters around the center of mass velocity. This could be an explanation to the quite low excitation energies of clusters found in central collisions for the Xe + Sn system at 50 A.MeV Marie98 . This subject will be more completely covered in a forthcoming article.

### III.3 Statistical description of collisions.

Let us end with fragment size analyses of cluster-cluster collisions. In nucleus-nucleus collisions studies, such analyses are very often used to fix the parameters of statistical models and to verify the compatibility of statistical decay models with experimental data (see for example Marie97 ; DAgostino99a ; Bocage ). The so-called Dalitz plot for central ($`0b/b_{max}<0.1`$) $`N=34`$ \+ $`N=34`$ collisions are shown on figure 7. Each panel corresponds to a fixed available energy in the center of mass ranging from $`E_{cm}/N=30`$ to $`120E.S.U.`$. For this analysis, only the three heaviest clusters are taken into account. Each event is associated with a point in this plot. The distance of one point with respect to each edge of the triangle is proportional to the size of each of the three clusters. The corners correspond to events with a large size cluster and two small ones (fusion/evaporation process), an event in the middle of an edge corresponds to an event with two equal size clusters and a small one (fission process or binary collision) and the center of the triangle corresponds to three equal size clusters (multi-fragmentation process). One can see that when the available energy increases, the reaction mechanisms goes continuously from fusion/evaporation to binary collisions and finally to multi-fragmentation and/or vaporization.

This picture is qualitatively very similar to the one obtained in SMM calculations SMM , describing the decay of a single source, as shown on figure 8. In this case, each panel corresponds to an excitation energy. As in figure 7, the lowest energies correspond to fusion/evaporation like processes, and goes with a continuous transition towards fission and multi-fragmentation when the excitation energy increases.

The strong qualitative similarity between these two pictures suggests that it could be possible to build a statistical decay code which would give the same results as the dynamical simulation. Of particular interest would be the study of the relations between the parameters of the statistical decay code (size of the source, excitation energy, deformation, radial flow, freeze-out volume, external constraints) with the parameters of the dynamical one (projectile and target sizes, impact parameter, parameters of the two-body interaction). This could also allow to check under which conditions the statistical decay code can be applied. This study could finally allow us to have a consistent description of cluster-cluster collisions, where the parameters of the statistical decay model are deduced from the parameters of the dynamical simulation. In this case, the cluster-cluster collisions could be described completely with a reduced number of parameters.

## IV Conclusions

As it has been seen in this brief overview of the characteristics and the reaction mechanisms of classical $`N`$ bodies in strong interaction, there are strong qualitative similarities between these systems and the atomic nuclei. The static properties and the reaction mechanisms observed for the atomic nuclei and for these classical clusters are found very close to each other. This could mean that the experimental observations made for nucleus-nucleus collisions are mainly governed by the N-body character of the system and by the overall shape of the two-body interaction (finite range attractive part and short range repulsive part). More physical ingredients (Coulomb interaction, quantum mechanics) are of course necessary if one wants to have a quantitative agreement with experimental data. But it is surprising to have such a qualitative agreement while essential physical ingredients are missing in the simulation.

Since all kind of reaction mechanisms are observed in these classical simulations, from the low energy fusion/evaporation processes to the high energy participant/spectator processes, they may also allow to connect in a consistent way the mean-field and the nucleon-nucleon approaches. One could for example study the relation between the incompressibility modulus and the size of the repulsive part of the interaction ($`K_{\mathrm{}}`$ and $`\sigma _{nn}`$ in the nuclear case). It is well known in transport calculations that one can “stiffen” the effective equation of state by just increasing the value of the nucleon-nucleon cross section $`\sigma _{nn}`$ DLM92 .

Finally, such simulations may reconcile two approaches which were up to now opposed in nuclear physics: the dynamical description and the statistical description of fragment emissions. The dynamical models are unique tools to establish the link between the parameters of the interaction and those of the statistical models. Additionally, one can determine under which conditions the statistical models can be used, and what are the true meanings of the parameters used in statistical models (temperature, chemical potentials, freeze-out volumes, etc.). This study may allow to have a complete and consistent description of colliding N-body systems with a reduced number of free parameters. This could be not only useful for nuclear physics, but also for other fields of physics like for example the cluster physics.

## Acknowledgments

I would like to thank my colleagues of the INDRA collaboration and from the Laboratoire de Physique Corpusculaire de Caen for a careful reading of this article. I especially thanks Denis Lacroix for suggestions and fruitful discussions. I would like also to thank O.Lopez and A.Botvina for the SMM calculations.


format: text

text:
# 1 Introduction

## 1 Introduction

Coherent X–ray pulsations from 1E 2259$`+`$586 were discovered at the end of seventies. Only few years ago it was recognized however that a number of X–ray pulsators, including 1E 2259$`+`$586, posses peculiar properties which are very much at variance with those of known accreting pulsars in X–ray binaries. These objects, initially suggested as a homogeneous new class of pulsators (Mereghetti & Stella 1995; van Paradijs et al. 1995), have been named in different ways, reflecting our ignorance of their nature: Very Low Mass X–ray Pulsars, Braking Pulsars, 6 s Pulsars, Anomalous X–ray Pulsars. The latter designation (AXPs) has became the most popular and will be used hereafter.

The nature of AXPs is one of the most challenging unsolved problems of Galactic high energy astrophysics. Over the last few years there has been a remarkable observational and theoretical effort to unveil their nature. Although we can be reasonably confident that AXPs are magnetic rotating neutron stars (NSs), their energy production mechanism is still uncertain; it is also unclear whether they are solitary objects or are in binary systems with very low mass companion. As a consequence, different production mechanisms for the observed X–ray emission have been proposed, involving either accretion or the dissipation of magnetic energy.

The properties that distinguish AXPs from known magnetic ($``$ 1012 G) accreting X–ray pulsars found in High and Low Mass X–Ray Binaries (HMXBs and LMXBs) are the following:

* spin periods in a narrow range ($``$6–12 s) compared with the much broader distribution (0.069 – $``$104 s) observed in HMXRB pulsars;
* no conspicuous optical counterparts (see Section 6), with upper limits which rule out the presence of massive companions, like OB (super)giants or Be stars;
* very soft and absorbed X–ray spectra: a black body with characteristic temperature in the 0.4–0.6 keV range and a steep power–law with photon index in the 2.5–4 interval (suggesting that perhaps a large part of the total luminosity is hidden in the EUV band);
* relatively low X–ray luminosity ($`10^{34}`$–1036 erg s-1) compared with that of HMXB pulsars;
* relatively low flux variability on timescales from hours to years;
* relatively stable spin period evolution, with long term spin–down trend (see Section 5);
* a very flat distribution in the Galactic plane and three clear associations with supernova remnants (SNRs; suggesting a young population).

There are currently (February 2001) five ascertained members of the AXP class (see Table 1) plus one likely candidate. This review, after a brief presentation of theoretical models (Section 2), concentrates on the recent results inferred from observations of AXPs at different wavelengths (X–rays, Section 4 and 5; optical and IR, Section 6; radio, Section 7). The implications of these results for the proposed theoretical models are briefly discussed.

## 2 Theoretical Models

Theoretical models for AXPs can be classified into two main classes depending on the mechanism that is supposed to power their X–ray emission: accretion or magnetic field decay. The former class includes both isolated NSs and NSs in binary systems.

Accretion from matter in a overdense region was first proposed for 4U 0142$`+`$614 based on its apparent association with a molecular cloud. However this would imply that AXPs move at low velocities in the ambient medium (Israel et al. 1994). Mereghetti & Stella (1995) proposed that AXPs form an homogeneous subclass of accreting neutron stars, perhaps members of very low mass X–ray binaries (VLMXBs), which are characterized by lower luminosities and higher magnetic fields ($`B`$$``$1011 G) than accreting neutron stars in classical LMXBs. On the other hand, Van Paradijs et al. (1995) proposed that AXPs are the result of a common envelope and spiral–in evolution of a neutron star and its massive companion. This ends up in the complete disruption of the companion star after the so–called Thorne–Zẏtkov stage. Based on the AXPs–SNRs association, Chatterjee et al. (2000) and Perna et al. (2000), proposed instead that AXPs accrete from a fossil disk made of matter falling back onto the neutron star after its birth. Alpar (2000) proposed that AXPs are NSs accreting from the debris of their SNRs, in a phase following the SGR stage, or, alternatively, NSs accreting from a companion star in a rare path of LMXB evolution. In these scenarios the P/Ṗ does not provide a reliable good estimator of the pulsar age; rather it should represent an asymptotic value approached by the pulsar. The equilibrium period should depend on the trapped angular momentum in the residual accreting matter (Alpar 2000). In the accretion scenario the narrow period distribution of AXPs can be accounted for by either limiting the magnetic field of the NS and strength of the propeller wind emission (Marsden et al. 1999) or using ADAF models and an appropriate distribution of magnetic field, initial spin and accretion disk mass (Chatterjee & Hernquist 2000).

Finally, Thompson & Duncan (1993, 1996) proposed that AXPs are “magnetars”, isolated neutron stars with a super-strong magnetic field ($`10^{1416}`$ G; for a review see also Thompson, this book). Magnetars were originally proposed to explain the properties of soft $`\gamma `$–ray repeaters. These were later determined to exhibit pulse periods (8.05 s, 7.5 s and 5.16 s in SGR 0526–66, SGR 1806–20 and SGR 1900+14, respectively) and period derivatives similar to those of AXPs (Kouveliotou et al. 1998; Hurley et al. 1999). If this connection proved correct, AXPs might be some sort of quiescent analogous of soft $`\gamma `$–ray repeaters. In the framework of the magnetar model, Colpi et al. (2000) showed that it is possible to account for the narrow period distribution of AXPs if the initial magnetic field (in the 1015–1016 G range) of magnetars decays significantly on a timescale of the order of 104 years.

## 3 The AXP sample

Table 1 lists the main characteristics of the five X–ray pulsars which form the AXP class. The new AXP candidate, AX J1844$``$0258, is also included.

Recently a possible connection of AXPs with Soft $`\gamma `$–Repeaters (SGRs, for a review see Kouveliotou, this book) has been proposed (Kouveliotou et al. 1998, 1999; Hurley et al. 1999); this builds on several observational similarities, such as the range of spin periods, their derivatives, and the possible association with SNRs. In the case of SGRs large offsets from the SNR centers were measured which might imply that AXPs may eventually evolve into SGRs as they age and move away from the SNR centre (Gaensler 2000). This would require a difference of 10–100 kyr in the age of AXPs and SGRs; therefore it would be difficult to explain, in the light of the inferred period derivatives, the similarity of periods. We also note that a 226 ms radio pulsar has been recently discovered in the SNR previously associated with the Soft $`\gamma `$–Repeater SGR 1900+14 (Lorimer & Xilouris 2000). Although it is unclear which of the two objects (if any) is associated to the SNR, a recent reanalysis of the distance indicators (Case & Bhattacharya 1998) places the SNR at 10$`\pm `$3 kpc, consistent with the position of the newly discovered radio pulsar. The alternative possibility is that SGRs are born with a different velocity distribution than AXPs and, therefore, the two classes cannot be drawn from the same parent population. The AXPs–SGRs connection is still an open issue.

## 4 The BeppoSAX view

Observations by different X–ray satellites have often been used to check for the presence of spectral and flux variations in AXPs. These, however, are affected by the uncertainties introduced by comparing the results of different instruments covering different energy ranges, etc. The BeppoSAX satellite, with a coherent set of instruments covering the 0.1–10 keV range with spectral resolution of 3–10, observed the whole sample of AXPs. Spectral results from BeppoSAX observations have already been reported for three AXPs, 1E 2259$`+`$586 (Parmar et al. 1998), 1E 1048.1$``$5937 (Oosterbroek et al. 1998), and 4U 0142$`+`$614 (Israel et al. 1999). In the following we will report on the spectral and timing properties of another AXPs, 1RXS J170849$``$400910 (Israel et al. 2001a), observed by BeppoSAX in March–April 1999. The BeppoSAX observation of 1E1841–045 is not included in this report due to strong contamination from the SNR Kes 73 (timing results have been reported by Gotthelf et al. 1999). The candidate AXP AX J1844$``$0258 has not yet been observed by BeppoSAX.

The relatively bright source 1RXS J170849$``$400910 was discovered by ROSAT at the beginning of its mission; however only in 1997 this source attracted much attention because of the discovery of $``$11 s pulsations with ASCA (Sugizaki et al. 1997). Based on the pulse period and unusually soft X–ray spectrum the source was tentatively classified as a candidate AXP. This interpretation was confirmed thanks to ROSAT HRI observations which provided the first measurement of the period derivative Ṗ$``$2$`\times `$10-11 s s-1 and a more accurate X–ray position from which it was possible, through optical imaging, to confidently exclude the presence of a massive companion star (Israel et al. 1999b). 1RXS J170849$``$400910 has been monitored with the RossiXTE since 1998; a sudden spin–up event, suggestive of a “glitch” from a highly magnetised NS, was recorded in September 1999 (Kaspi et al. 2000a; see also Section 5).

BeppoSAX observed this source on 31 March – 1 April 1999 with the Narrow field Instruments: the Low–Energy Concentrator Spectrometer (LECS; 0.1–10 keV; Parmar et al. 1997; 26 ks effective exposure time) and the Medium–Energy Concentrator Spectrometer (MECS; 1.3–10 keV; Boella et al. 1997; 52 ks effective exposure time). A simple power–law model (as well as any other single component model) did not fit well the data ($`\chi ^2`$ of 1.42 for 199 degrees of freedom). A better fit (see right panel of Figure 1) was obtained including a soft thermal component in analogy with similar results obtained for other AXPs ($`\chi ^2`$ of 1.12 for 197 degrees of freedom). The best fit was obtained for an absorbed, $`N_H`$=(1.46$`\pm `$0.02)$`\times `$1022 cm-2, power–law with photon index 2.6$`\pm `$0.2 and a blackbody component with a characteristic temperature of 0.45$`\pm `$0.03 keV (90% c.l. reported; Israel et al. 2001a). The unabsorbed 1–10 keV flux was 6$`\times `$10-11 erg s-1. The blackbody component accounts for about 30% of the total observed flux. Figure 1 compares the spectral shape and components of 1RXS J170849$``$400910 and 4U 0142$`+`$614 as seen by BeppoSAX. 1RXS J170849$``$400910 is therefore the fourth AXP for which a two componet spectrum (steep power–law plus soft blackbody) has been detected.

A further interesting result of the 1RXS J170849$``$400910 BeppoSAX observation was the detection of a clear pulse shape change as a function of energy (see left panels of Figure 2). In particular the minimum in the lowest energy interval (0.1–2 keV) corresponds to a maximum in the 6–10 keV band. The pulsed fraction (semiamplitude of modulation divided by the mean source count rate) decreases from 40% to 30% from the lowest to the highest energy band, respectively. Pulse phase spectroscopy was carried out with the MECS data (due to poor statistics at low energies the data from the LECS data were not used). By keeping the parameters of the blackbody component fixed at the phase–averaged value, a $``$3$`\sigma `$ significant variation in the power–law photon index was found. The large uncertainties prevent the detection of possible variations in $`N_H`$ (see right panels of Figure 2). Although small energy–dependent changes in the pulse shape were already suggested in the past for 4U 0142$`+`$614 (Israel et al. 1999a; Paul et al. 2000), the variations of 1RXS J170849$``$400910 are highly significant and likely arising from a changing power–law slope. The lack of any conspicuous change in the pulse profiles of AXPs was used in the past to argue against the possibility that these sources are accreting X–ray objects. More sensitive studies will yield additional information on the spectral changes causing the pulse shape variations, extend the energy range of the source detection above 10 keV, and allow to look for cyclotron signatures.

The spectral properties of 1RXS J170849$``$400910 reported above are plotted in Figure 3 (left panels) together with the BeppoSAX results obtained for the other three AXPs for which good spectral data are available. It is apparent that the spectral properties of AXPs are relatively homogeneous. 1E 1048.1$``$5937 has the highest blackbody temperature and largest (absorbed) $`L_{BB}`$/$`L_{tot_{\mathrm{abs}}}`$ ratio. Note that the $`L_{BB}`$/$`L_{tot_{\mathrm{abs}}}`$ ratio is not strongly dependent on the power–law photon index. The AXPs with the largest (1RXS J170849$``$400910) and smallest (1E 1048.1$``$5937) absorption have nearly the same photon index, suggesting that the detection of a steep power–law in AXPs is not an artifact due to a large $`N_H`$ value.

The picture becomes quite different when extrapolating the total unabsorbed fluxes in the 0.1–10 keV range and considering the bolometric blackbody flux; this gives in all cases a $`L_{BB_{\mathrm{bol}}}`$/$`L_{tot_{\mathrm{unabs}}}`$ ratio smaller than 30%. Despite the $``$103 ratio in total unabsorbed flux between the brightest (4U 0142$`+`$614) and faintest (1RXS J170849$``$400910) AXPs, the blackbody component varies by only a factor of $``$10, strongly suggesting that the energy budget of AXPs is largely dominated by the power–law. Another interesting result is shown in the right panel of Figure 3 where the 0.1–10 keV pulsed fractions are plotted as a function of $`L_{BB_{\mathrm{bol}}}`$/$`L_{tot_{\mathrm{unabs}}}`$. This correlation strongly suggests that the thermal component is mainly responsible for the pulsations component. This is also in agreement with the fact that 4U 0142$`+`$614, the brightest AXP, has also the lowest known pulsed fraction (Israel et al. 2001b).

## 5 Further X-ray observations

Several results from X–ray observations of AXPs and related objects were reported in the last year. In particular the pulsation stability of AXPs was compared to that of conventional accreting X–ray pulsars and isolated radio pulsars. 1RXS J170849$``$400910 and 1E 2259$`+`$586 were the first two AXPs for which a systematic study was carried out and phase–coherent timing solutions obtained by means of RossiXTE data (Kaspi et al. 1999). These two sources were found to be quite stable rotators with phase residuals of only $``$1%, comparable to or smaller than those measured for most radio pulsars. However, in September 1999 the RossiXTE satellite detected a sudden spin–up event in 1RXS J170849$``$400910 which was originally interpreted as a “glitch” similar to that observed in the Vela and other young radio pulsars (Kaspi et al. 2000a). However, we note that in principle glitches could also be detected in accreting spinning–down X–ray sources (with a sufficiently high magnetic field strength) if they are in a low noise level phase, as indeed AXPs are known to be. A way to distinguish, in the near future, whether 1RXS J170849$``$400910 experienced a radio pulsar–like glitch or, perhaps, an accreting X–ray pulsar–like spin–up behaviour would be to accurately monitor the period history after the event. Unfortunately these results have not been reported so far.

Also the pulsations stability of 1E 1048.1$``$5937 was studied in great detail with ASCA (Paul et al. 2000) and RossiXTE (Kaspi et al. 2000b) data. In the latter case the sampling proved insufficient to find a phase–coherent solution due to a high noise level (changing spin–down rate). The observed deviations from simple spin–down were found to be inconsistent with a single glitch event. Such period changes do not seem to be accompanied by pulse shape or spectral changes (in ASCA) or even flux variations as observed in accreting X–ray pulsars (Paul et al. 2000; Kaspi et al. 2000b). In this respect 1E 1048.1$``$5937 was suggested to be a transition object between SGRs and AXPs (Kaspi et al. 2000b). Based on earlier RossiXTE data, Baykal et al. (2000a) found that the level of pulse fluctuations in 1E 1048.1$``$5937 is consistent with the typical noise level of accretion powered pulsars (Baykal & Ögelman 1993; Bildsten et al. 1997). Also the results of further RossiXTE observations of 1E 2259$`+`$586 suggest that the source is perhaps in a fairly stable accretion phase with a constant X–ray luminosity and spin–down rate (Baykal et al. 2000a). An interesting result has been recently reported for the high–mass accreting X–ray pulsar 4U 1907+09. This is the only known X–ray pulsar (besides AXPs) which has always been spinning–down since the discovery of pulsations at 440 s in its X–ray flux in 1983. This result, based again on RossiXTE data, shows that 4U 1907+09 during 1996–1998 has been a rather stable rotator with a spin–down rate only a factor four greater than that of 1E 2259$`+`$586 (Baykal et al. 2000b). Furthermore the long–term noise strength of 4U 1907+09 is one order of magnitude lower than that of 1E 2259$`+`$586; the conclusion is that the existence of an accreting source displaying persistent spin–down shows that quiet spin–down trends do not necessarily imply that the sources are not accreting.

No new results have been reported concerning the search of orbital signatures by means of delays in the pulse arrival times of AXPs. In fact after the tight upper limits inferred on the possible companion star mass of 1E 1048.1$``$5937 and 1E 2259$`+`$586 (Mereghetti et al. 1998), and 4U 0142$`+`$614 (Wilson et al. 1999) based on RossiXTE data, none of the other three AXPs has been studied in this respect. The lack of detectable Doppler effect has been used by many authors to favor the magnetar scenarios. It is worth mentioning that the pulse arrival time delays remain undetected also in 4U 1626–674, an X–ray pulsar with a very low–mass companion star (likely a light He–burning or He white dwarf; see Chakrabarty 1998) in a 42 min orbital period binary. If AXPs were in binary systems with companion stars similar to that of 4U 1626–674, expected delays would be below the current upper limits.

Finally Marsden & White (2001) recently found a correlation between spin–down rates of AXPs and SGRs and the L$`_{\mathrm{BB}_{\mathrm{bol}}}`$/LPL ratio concluding that, regardless the nature of these sources, they are likely objects with similar emission mechanisms.

## 6 Optical/IR follow–up observations

Similar to other Galactic X–ray sources, the identification of AXPs at optical and IR wavelengths suffers from the comparatively poor spatial resolution of X–ray telescopes which is often inadequate to cope with the crowded Galactic plane fields in which AXPs lie, and the strong extinction in their direction. Although the X–ray positional accuracy is not yet sufficient to unambiguously identify the possible optical/IR counterparts of AXPs, it is now possible to sort out a number of good candidates or, at least, put tight constraints on the optical/IR–to–X–ray energy distribution of these sources.

This is the case of the recent proposed optical counterpart of 4U 0142$`+`$614, a faint relatively blue object ($`R`$=$`25`$; $`VR`$=$`0.63`$) which, in the color–magnitude diagram, lies half–way between the main sequence and the track of 0.6 M white dwarfs (Hulleman et al. 2000a; Keck observations). However the optical counterpart candidate position of 4U 0142$`+`$614 (within the Einstein HRI error circle; White et al. 1987) is outside the EXOSAT LEIT error circle (White et al. 1987) and at the edge of the ROSAT PSPC one. A more accurate X–ray position is clearly needed in order to confirm this result. Since none of the accretion–based scenarios considered by Hulleman et al. (2000a) fit the optical and X–ray data these authors conclude that the measurements may be in agreement with a magnetar. However no detailed models have yet been developed for the optical emission of magnetars (Hulleman et al. 2000a). We note that the optical data are also in agreement with at least two other scenarios: (a) a binary system hosting a white dwarf companion, which accounts for the optical emission, and an accreting neutron star producing the X–ray emission, (b) a disk around an isolated NS with a more realistic choice of values for the disk size, illumination and inclination (Israel et al. 2001c).

Hulleman et al. (2000b) carried out also a relatively deep search for the optical counterpart of 1E 2259$`+`$586 based on Keck observations. No object was found down to limiting magnitudes of $`R`$=$`25.7`$ and $`I`$=$`24.3`$. These authors conclude that it is unlikely that 1E 2259$`+`$586 is an isolated NS accreting from a residual disk. Also in this case, however a binary system with a helium–burning 0.3 M or a white dwarf companion star, or a magnetar scenario cannot be excluded yet.

Optical observations of the field of 1RXS J170849$``$400910 were obtained by Israel et al. (1999b). These authors found that the possible counterpart cannot be a massive early type star (a distant and/or absorbed OB star would appear more reddened). However the images were taken from a 1.5 m telescope and are not deep enough to constrain any other proposed theoretical scenario such as a low mass companion, a residual disk or a magnetar.

Recently Mereghetti et al. (2001) reported on the first search for optical/IR counterpart in the field of 1E 1841$``$045 located at the center of the SNR Kes 73 (mainly with the 1.5 m and 3.5 m class ESO telecopes). Similar to the previous case, the results are not very constraining, due to the high extinction in the direction of the source ($`A_V`$$``$11). No detailed reports on optical/IR observations of AX J1844$``$0258 have yet been published. Table 2 summarises the current optical/IR upper limits and measurements of AXPs.

## 7 Radio follow–up

If AXPs are isolated NSs that emit radiation in a fashion similar to rotation powered pulsars then they might shine and pulsate also in the radio band. We observed the fields of four southern AXPs, namely 1E 1048.1$``$5937, 1E 1841$``$045, AX J1844$``$0258 and 1RXS J170849$``$400910, with the Parkes Observatory and typical exposures of 10–20 ks for each source. The sampling time was in the 1–1.2 ms range with a beam aperture of $``$14’ (at 1.4 GHz). Sampled duty cycles were in the $``$0.001–20% range, while dispersion measures (DM) up to 10 times larger than the Galactic values in the direction of observed AXPs were sampled in order to take into account also the possible presence of local matter. Setting the detection threshold at a signal to noise ratio of 10, an average value of $``$70$`\mu `$Jy at 1.4 GHz can be reasonably assumed for the upper limit to the radio flux of the sample of observed AXPs observed. The search for coherent signals at the X–ray period (we used the $`P`$ and Ṗ values reported in Table 1 to infer the period interval for the search) did not detect any significant signal (Burderi et al. 2001).

## 8 The Future

After more than 20 years from the discovery of pulsations from 1E 2259$`+`$586 the nature of AXPs is still uncertain. Important limiting factors in the study of AXPs are: (i) the relatively low X–ray positional accuracy provided so far by X–ray telescopes together with the difficulties in obtaining deep optical/IR images (needed to sample faint objects) from the largest ground–based telescopes, (ii) the lack of spectral information above 10 keV where many cosmic X–ray sources (X–ray pulsars especially) display important features, (iii) the paucity of unambiguous predictions by theoretical models, and (iv) the small numbers in the currently known AXP population. The present generation of X–ray astronomy satellites (Chandra and NewtonXMM) and large ground–based telescopes (VLT, Gemini, Subaru and KeckII) opens new prospectives in the field. The unrivalled spatial resolution and absolute source positioning accuracy ($``$1”) of Chandra will provide the best chance to unambiguously identify the counterparts of AXPs. Moreover assessing that the AXP candidate AX J1844$``$0258 (which has shown a factor of variability larger than 10) is indeed an AXP would imply that the number of AXPs and their formation rate have been largely underestimated. The high throughput of NewtonXMM and the relatively wide field of view of its instruments will presumably allow to find other AXPs even at fainter fluxes, thus providing new important informations to understand the puzzling nature of AXPs.

## Acknowledgements

The results reported in this review have been obtained with the contribution of a number of scientists, who are part of large international collaboration aimed at studying the AXPs. It is a pleasure to thank all the people who have contributed to this project: L. Angelini, L. Burderi, S. Campana, S. Covino, N. D’Amico, F. Haberl, G. Marconi, R. Mignani, I. Negueruela, T. Oosterbroek, A.N. Parmar, A. Possenti, and N.E. White.


text:
# Electrons on the double helix: optical experiments on native DNA

## Abstract

Optical experiments on calf thymus DNA films subjected to different buffer environments are reported. The optical conductivity is that of a disordered or lightly doped semiconductor with a well-defined band-gap for charge excitations and low frequency transport determined by a small number of strongly localized electron states.

The question whether electrons are delocalized along the DNA duplex has attracted substantial recent interest and controversy. Both short and long-range migration was found by studies that were chemical in nature. DC conductivity measurements on single DNA strands led to equally controversial results - with metallic, semiconducting and insulating behavior all observed. Recent experiments  are also suggestive of a proximity-induced superconductivity at low temperatures. DNA networks in contrast were found to be highly resistive. Our earlier experiments conducted on $`\lambda `$-DNA in the microwave spectral range gave evidence for a large resistance associated with the duplex and for a thermally driven transport process.

These conflicting results may, to a large extent be due to different DNA varieties. A native DNA duplex with random (or nearly random) base pair sequences is expected to have electron states that are different than an oligomer with identical base pairs, such as a poly(C)-poly(G) track. In addition, the DNA configuration is sensitive to the buffer environment, which surrounds the duplex. In a dry form the duplex has non-parallel base pairs, and at higher water concentrations the duplex undergoes a structural change to a more ordered system. This then may lead to widely different transport properties and also charge excitations of a different nature.

In order to examine the nature of electron states in native DNA, we have conducted optical measurements spanning a wide spectral range, from microwave frequencies to the UV part of the electromagnetic spectrum, on oriented films, fabricated from DNA extracted from calf thymus. Our findings concerning the vibrational modes will be discussed elsewhere. Here we focus on the electronic excitations at high and at low energies. Our data gives evidence for well-defined charge excitations at high energies involving the p-orbitals of the base pairs, and at low energies charge excitations displaying all the characteristics of conduction due to a small number of localized electron states. Thus our experiments suggest that DNA is a wide bandgap semiconductor, with intrinsic disorder, counterion fluctuations, and possibly other sources leading to a small number of localized electronic states on the base pair sequence.

Wet-spun, free-standing oriented samples were prepared from calf thymus Li-DNA (Pharmacia) with a molecular weight of $`10^7`$ (corresponding to contour length of 5 mm or some 100 persistence lengths) by a method described by Rupprecht et al. This spinning allows controlled production of sufficient amounts of highly oriented thin films by spooling DNA fibers that are continuously stretched during precipitation into an aqueous alcohol solution. Films of thickness of 0.2 mm and lateral dimensions between 5 and 10 mm were used.

In all applicable cases placing it in a chamber with appropriate relative humidity controlled the hydration of the DNA sample. Several configurations were employed to evaluate the optical conductivity $`\sigma _1\left(\omega \right)`$. In all cases we treat the collection of DNA strands as a collection of thin wires, of diameter 2 nm, and we define the conductivity $`\sigma `$ as $`j/E`$ where $`j`$ is the electric current density induced along the helix axis. For randomly coiled DNA strands the loss due to motion of electric charges $`W`$ is, to a good approximation, given by $`W=\frac{1}{3}V\sigma \left(E_0\right)^2`$, where $`V`$ is the volume of the sample and $`E_0`$ is the time averaged applied ac field at the position of the sample. At 12 GHz the conductivity was evaluated from the measured loss of highly sensitive resonant cavities that were loaded with the material. The technique and the analysis, which leads to evaluation of the conductivity from the measured losses, is well established. In the millimeter spectral range, 100 GHz to 1000 GHz, backward wave oscillators were employed as coherent sources in a transmission configuration. The oriented calf-thymus DNA film was placed on a 1 mm thick sapphire substrate and then held in place by a sheet of 6 $`\mu `$m thick mylar to form a three layer system. Transmission as a function of frequency was recorded in the specified frequency range. In our analysis we utilized the fact that for plane waves incident normally on a slab of material, resonances occur whenever the slab is an integer number of half wavelengths. Thus, using the sapphire as a substrate, resonances occurred approximately every 50 GHz. Having analyzed the transmission through the sapphire alone prior to mounting the sample, the optical properties of the substrate were well characterized. The index of refraction of mylar was taken to be approximately 1.5, and its extinction coefficient was neglected. Thus using a three-layer transmission model, each resonance was analyzed for the optical properties of the DNA film, allowing for a 1.5 $`cm^1`$ resolution of the spectrum. Thin tungsten wire grids with a spacing much smaller than the wavelength of our radiation acted as adjustable polarizers to probe the sample anisotropy. Transmission measurements in the UV spectral range were conducted using a Beckman Coulter DU 640 Spectrophotometer.

The optical conductivity, measured over a broad spectral range is displayed in Fig. 1. We have omitted the spectral range between 500 and 10,000 $`cm^1`$, where various intramolecular vibrational excitations occur, which will be the subject of a separate publication. In the figure we have also displayed the conductivity extracted from the transmission data of Wittlin et al. conducted on dry specimens. Several comments on the experiments we have conducted are in order. First, during experiments with an electric field polarized parallel and perpendicular to the duplex axis we did not observe a substantial anisotropy of the conductivity, and the absorption peaks associated with vibrational modes were also found to be isotropic. While this is surprising, it can be explained by assuming that the DNA duplex while oriented macroscopically does not assume a straight configuration but displays a substantial local directional variation. This has been observed in certain DNA films. As this issue is unresolved, we have not included a factor of $`\frac{1}{3}`$ in the loss equation; this however does not affect our overall conclusion. Second, it is evident from the figure that data at low frequencies, measured on dry DNA species is reproducible, that is, the conductivity we have evaluated is, within experimental error identical to those found by others . Third, the low frequency optical conductivity depends on the water environment with the conductivity in a wet environment significantly larger than that in a dry environment. We have found similar results earlier for $`\lambda `$-DNA.

The mode with the onset at 30,000 $`cm^1`$ is due to intra-base electronic excitations associated with the $`\pi \pi ^{}`$ molecular orbital transitions. This we have confirmed by conducting optical absorption measurements on the individual A, T, C and G bases; the absorption as shown in Fig. 1 is virtually identical to the sum of the optical transitions associated with the four base pair species. The spectral weight

$$\sigma \left(\omega \right)𝑑\omega =\frac{\pi Ne^2}{2m}$$
(1)
of the mode, where $`N`$ is the concentration of charge carriers, and $`m`$ is the (electronic) mass, has been evaluated by using the measured extinction coefficient. By a numerical integration of the spectra for samples of known dimensions and concentration of DNA, we obtain a value for the concentration of charge carriers of order $`N10^{21}cm^3`$, which compares favorably with the concentration evaluated by assuming that there is one electron per base associated with this transition, $`N=2\times 10^{21}cm^3`$.

On the basis of this analysis we conclude that the absorption feature in the UV spectral range represents nearly all the spectral weight associated with the electronic excitations of the base pairs of the DNA duplex. We note that this excitation cannot be associated with the bandgap in the usual sense, as this mode does not represent a transition between the highest occupied and lowest unoccupied states associated with the entire electronic structure of the DNA duplex. The reason for this is the following: for a duplex, the optical transition we observe corresponds to a transition between energy levels of the various single bases, i.e. intra-base excitations, while the transition matrix element involving energy levels of different bases (such as A to T or C to G optical transitions) is vanishingly small. The bandgap, on the other hand, corresponds to the energy difference between the top of the HOMO band and the bottom of the LUMO band, with these bands in general corresponding to different bases.

The optical conductivity below about 500 $`cm^1`$, represented by the full line in Fig. 1 is due to low frequency vibrations involving the double helix structure. These vibrations have been described earlier and are not discussed here. In Fig. 2 the frequency dependent conductivity, as measured below these vibrational modes in the micro and millimeter wave spectral range is displayed for both dry and fully hydrated DNA. The dashed lines represent a conductivity given by the expression

$$\sigma _1\left(\omega \right)=A\left(T\right)\omega ^\alpha $$
(2)
with the exponent $`\alpha `$ given in the figure. Such a power law dependence has been observed in a variety of disordered solids, ranging from ionic glasses to materials where the electron states are localized . We believe that contributions to the conductivity due to the counter ions and water molecules, which surround the DNA duplex, can be neglected, and that the conductivity is due to localized electrons or holes, for several reasons. First, both the counter ions and the water in the hydration layers are strongly bound to the (negatively charged) DNA duplex and consequently excitations due to these should occur only at higher frequencies in the infrared spectral range and above. Second, the number of counter ions does not vary during the hydration, and thus they cannot be responsible for the strong increase of the conductivity with the increasing water content. The increased conductivity measured for the hydrated DNA may, in principle, originate from the response of the water molecules surrounding the duplex in the hydration layer. We believe this unlikely as we have measured the dependence of the optical conductivity as a function of humidity, and have found a strongly nonlinear dependence of said conductivity on the number of water molecules surrounding the duplex. Experiments on frozen samples of the same composition, where water molecules and counter-ions are immobilized, show similar behavior indicating that the contribution of the counter-ions to the optical conductivity is negligible. A more likely origin of the difference between the conductivity of the dry and water saturated DNA is as follows. In a water rich or high humidity environment, the DNA duplex takes on a more ordered structure (bases stacked parallel to each other) referred to as its B-form as compared to when the DNA duplex is in a dry or low humidity environment, where the structure is less ordered (bases are stacked with different angles with respect to the main axes of the molecules) and referred to as its A-form. We believe that the increased conductivity is due to the more ordered arrangement of the DNA system in its B-form as compared to the disordered A-form allowing for increased electronic charge transport in spite of the fact that distances between neighboring bases are shorter in the A-form. In Fig. 3 we display the temperature dependence of the conductivity measured at 12, 150 and 500 GHz as a function of temperature. Several features are of importance. First we find a strongly temperature dependent conductivity, and the temperature dependence is typical of a conductivity $`\sigma _1`$ determined by temperature driven transport processes. Second, the temperature dependence itself depends only weakly on the frequency over a substantial frequency range.

The temperature and frequency dependence of the conductivity as observed at low frequencies has all the hallmarks of a transport process that is determined by transitions between localized electronic states. The frequency dependence of the conductivity under such circumstances is well described by Eq. (2) with the exponent $`a`$ depending on the conduction process. For instance, in the case of variable range hopping (VRH) the strength of the electron-electron interactions can change the exponent $`a`$. For VRH, in the case of non-interacting electrons one finds an approximately linear dependence of the conductivity on frequency whereas in the case of interacting electrons, one finds a quadratic dependence of the conductivity on frequency (plus logarithmic corrections). We also find a conductivity that is well described by Eq. (2) with exponents different for dry and water saturated DNA. The reason for this difference is not clear; nevertheless the frequency dependence observed is clearly similar to that of a random assembly of charged entities. In the parameter region $`k_BT>\mathrm{}\omega `$, the prefactor $`A\left(T\right)`$ in Eq. (2), should display a non-exponential temperature dependence. The function $`A\left(T\right)`$ depends on the overall energy scale $`E_0`$ associated with the localized states. Variable range hopping is the dominant dc conduction mechanism if $`E_0>k_BT`$, for which

$$\sigma _{DC}=B\mathrm{exp}\left[\left(\frac{T_0}{T}\right)^\beta \right]$$
(3)
in one dimension, with $`\beta =\frac{1}{2}`$, no matter the strength of electron-electron interactions, and the characteristic temperature $`T_0`$ depending on the localization length and also on the density of states of the carriers. The temperature dependence is non-exponential, and in addition becomes progressively weaker with increasing frequency. The measured temperature dependence displayed in Fig. 3 is in qualitative accordance with such behavior. This then strongly suggests states localized by disorder. Both static disorder associated with the random base pair sequences, and also fluctuations involving the DNA duplex may be responsible for determining the overall temperature and frequency dependence.

Our experimental results, obtained over a broad spectral range strongly suggest that native DNA is a wide bandgap semiconductor, with disorder and counterion fluctuations (among other possible sources) leading to a small number of localized carriers on the base pairs - a situation not unlike what occurs in lightly doped crystalline or amorphous semiconductors or doped polymers. Under such circumstances two important contributions to the optical conductivity emerge: a well defined transition associated with the electron states of the ”pure” system (this transition occurring at finite energy) together with a low frequency contribution to the conductivity with prominent frequency and temperature dependencies. We observe both features in native DNA with some modifications from what is found for a typical lightly doped (crystalline or amorphous) semiconductor or polymer. In the latter case one observes a well-defined optical transition between the valence and conduction band, as discussed earlier, for DNA, the transition is between orbitals associated with individual bases. Thus the fact that band structure calculations give a bandgap smaller than that corresponding to the energy for the onset of absorption in Fig. 1 is not surprising.

The number of carriers involved in the low frequency transport process is difficult to estimate as the parameters that enter into the equations for the temperature and frequency dependent conductivity are determined by both the number of carriers, and by the localization length. A comparison with other disordered one-dimensional conductors however allows one to draw some conclusions and to make some order of magnitude estimates. The temperature dependence of the conductivity that we find for DNA is similar to a strongly disordered organic conductor $`Qn(TCNQ)_2`$ with disorder induced by irradiation, for which the localization length is $`\xi 10`$Å, and the number of carriers is $`N10^{22}cm^3`$. The magnitude of the conductivity we find here is approximately three orders of magnitude smaller than that observed in $`Qn(TCNQ)_2`$. For a conduction process determined by temperature driven transitions between localized electron states, the localization length, $`\xi `$, mainly determines the overall temperature dependence. The magnitude of $`\sigma _1`$ reflects both $`\xi `$ and the number of carriers, $`N`$. Thus the comparison between calf-thymus DNA and $`Qn(TCNQ)_2`$ suggests that in native DNA electron states are characterized by a short localization length not exceeding one lattice constant (the distance between base pairs), while the number of carriers participating in the low frequency transport process is of the order of $`N10^{19}cm^3`$. The higher conductivity of DNA in a water rich environment may reflect a different carrier number, but most likely it is due to the more regular B-form of DNA which occurs for high water content, thus also leading to weaker localization effects, and thus to a larger localization length.

The source of the charge entities (electrons or holes) which are associatedwith the low frequency conductivity is not obvious. The random base pair sequences, together with disorder associated with the finite persistence length, and counter ion fluctuations all may lead to a small number of localized charges on the base pair stack. The fact that the states, which contribute to the low frequency optical response, are localized has different ingredients. The random base pair sequences which occur in native DNA lead to a random potential along the duplex, and thus to charge localization. The potential energy fluctuations associated with the base pair sequences are of the order of $`0.5eV`$, about ten times larger than the overlap integral between the electron states on neighboring base pairs along the double helix. Under such circumstances the localization length is of the order of one lattice constant, comparable to the value we inferred before. One should also note that equally important is the dynamic disorder associated with base pair fluctuations, of which the influence on charge transport along the DNA duplex has been conjectured earlier.

We wish to thank Phu Tran for assisting with the cavity measurements. This research was supported by the National Science Foundation grant DMR-0077251.


text:
# Multiband model for tunneling in MgB2 junctions

## Abstract

A theoretical model for quasiparticle and Josephson tunneling in multiband superconductors is developed and applied to MgB2-based junctions. The gap functions in different bands in MgB2 are obtained from an extended Eliashberg formalism, using the results of band structure calculations. The temperature and angle dependencies of MgB2 tunneling spectra and the Josephson critical current are calculated. The conditions for observing one or two gaps are given. We argue that the model may help to settle the current debate concerning two-band superconductivity in MgB2.

Soon after the discovery of superconductivity in MgB2 akimitsu , first principle calculations were performed to determine the electronic structure of this material. It was found that 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 Kortus ; An ; Kong ; Bohnen . The multiband picture has given rise to the concept that two superconducting energy gaps can coexist Liu ; Shulga  in MgB2.

Two-band superconductivity is a phenomenon that has been observed in Nb doped SrTiO3 Bednorz . Recent experimental STM and point-contact spectroscopy Klein ; Bugoslavsky ; Szabo , high-resolution photo-emission spectroscopy Tsuda , Raman spectroscopy Chen , specific heat measurements sph  and muon-spin-relaxation studies of the magnetic penetration depth Nieder  support the concept of a double gap in MgB2 (see Ref. Buzea,  for a review of experiments). However, there is an ambiguity in the interpretation of point-contact data concerning the existence of two gaps Klein ; Bugoslavsky ; Szabo . Moreover, some tunneling measurements Gonelli  show only one gap with a magnitude smaller than the BCS value of $`\mathrm{\Delta }=1.76`$ $`k_BT_c`$.

In order to resolve this discrepancy, we address the question, how multiband superconductivity will manifest itself in tunneling. We present the theoretical model for quasiparticle and Josephson tunneling in MgB2-based junctions. Using the results of band-structure calculations, we apply an extended Eliashberg formalism to obtain the gap functions in different bands, taking strong coupling effects into account. Tunneling from a normal metal (N) into MgB2 is considered in an extended Blonder-Tinkham-Klapwijk (BTK) model BTK . The temperature dependencies and absolute values of the $`I_cR_N`$ product ($`I_c`$ is the critical current and $`R_N`$ is the normal state resistance) are calculated in MgB2-based SIS tunnel junctions, where S denotes a superconductor and I an insulator. Tunneling in the direction of the $`a`$-$`b`$ plane, in the $`c`$-axis direction and under arbitrary angle is considered. Furthermore, the Josephson supercurrent between a single-gap superconductor and MgB2 is calculated.

According to the labeling of Liu et al. Liu , the four Fermi surface sheets in MgB2 are grouped into quasi-two-dimensional $`\sigma `$ bands and three-dimensional $`\pi `$ bands. Hence, normal and superconducting properties of MgB2 can be described by an effective two-band model. Within this model, Liu et al. Liu  estimated the coupling constants and energy gap ratio in the weak coupling regime. More recently, the band decomposition of the superconducting and transport Eliashberg functions $`\alpha _{ij}^2F_{ij}\left(\omega \right)`$ (where $`i`$ and $`j`$ denote $`\sigma `$ or $`\pi `$ bands), which describe the electron-phonon coupling in MgB2 as function of the frequency $`\omega `$, was provided in Ref. Dolgov, . This allows to perform a strong coupling calculation of the superconducting energy gap functions $`\mathrm{\Delta }_i\left(\omega _n\right)`$ in different bands. The functions $`\mathrm{\Delta }_i\left(\omega _n\right)`$ in turn determine the Josephson critical current in a tunnel junction between multiband superconductors, which is given by a straightforward generalization of the well-known result AB  to the case of several conducting bands MGZ  as well as strong coupling. The critical current component for tunneling from band $`i`$ into $`j`$ is given by

$$I_{ij}=\frac{\pi T}{eR_{ij}}\underset{\omega _n}{}\frac{\mathrm{\Delta }_i\left(\omega _n\right)\mathrm{\Delta }_j\left(\omega _n\right)}{\sqrt{\omega _n^2+\mathrm{\Delta }_i^2\left(\omega _n\right)}\sqrt{\omega _n^2+\mathrm{\Delta }_j^2\left(\omega _n\right)}},$$
(1)
where $``$ and $``$ denote left and right superconductors respectively, $`R_{ij}^1=\mathrm{min}\{R_{ij}^1,R_{ij}^1\}`$ is the normal-state conductance of a junction for the bands $`(i,j)`$ which is given by the integral over the Fermi surface $`S_{i(j)}`$

$$\left(R_{()ij}𝒜\right)^1=\frac{2e^2}{\mathrm{}}_{v_x>0}\frac{D_{ij}v_{n,i(j)}d^2S_{i(j)}}{\left(2\pi \right)^3v_{F,i(j)}},$$
(2)
where $`𝒜`$ is the junction area, $`v_n`$ is the projection of the Fermi velocity $`v_F`$ on the direction normal to the junction plane, and $`D_{ij}`$ is the probability for a quasiparticle to tunnel from band $`i`$ in $``$ into band $`j`$ in $``$. The total critical current is the sum of the components $`I_c=_{ij}I_{ij}`$.

The gap functions $`\mathrm{\Delta }_i\left(\omega _n\right)`$ can be calculated with an extension of the Eliashberg formalism Carbotte  to two bands

$$\mathrm{\Delta }_i\left(\omega _n\right)Z_i\left(\omega _n\right)=\pi T\underset{j}{}\underset{\omega _m}{}\frac{\left(\lambda _{ij}\stackrel{~}{\mu }_{ij}^{}\right)\mathrm{\Delta }_j\left(\omega _m\right)}{\sqrt{\omega _m^2+\mathrm{\Delta }_j^2\left(\omega _m\right)}},$$
(3)
$`Z_i\left(\omega _n\right)=1+{\displaystyle \frac{\pi T}{\omega _n}}{\displaystyle \underset{j}{}}{\displaystyle \underset{\omega _m}{}}\lambda _{ij}{\displaystyle \frac{\omega _m}{\sqrt{\omega _m^2+\mathrm{\Delta }_j^2\left(\omega _m\right)}}},`$ (4)
where $`\lambda _{ij}=2_0^{\mathrm{}}\omega \alpha _{ij}^2F_{ij}\left(\omega \right)𝑑\omega /[\omega ^2+\left(\omega _m\omega _n\right)^2]`$, $`Z_i\left(\omega _n\right)`$ are the Migdal renormalization functions and $`\omega _n=\pi T(2n+1)`$. These equations are solved numerically with the electron-phonon $`\alpha _{ij}^2F_{ij}\left(\omega \right)`$ functions from Ref. Dolgov, . The cutoff frequency $`\omega _c`$ is taken equal to 10 times the maximum phonon frequency. The functions $`\stackrel{~}{\mu }_{ij}^{}`$ represent a matrix of the Coulomb pseudopotentials defined at $`\omega _c`$, calculated up to a common prefactor that is used as an adjustable parameter to get $`T_c`$ = 39.4 K. The matrix $`\mu ^{}`$ at the frequency $`\omega _{\mathrm{ln}}`$ (relevant for the McMillan expression for $`T_c`$ in the isotropic case) is given by $`\mu ^{}=[1+\stackrel{~}{\mu }^{}\mathrm{ln}(\omega _c/\omega _{\mathrm{ln}})]^1\stackrel{~}{\mu }^{}`$, where $`\omega _{\mathrm{ln}}`$ follows from $`0=_0^{\mathrm{}}\mathrm{ln}(\omega /\omega _{\mathrm{ln}})\omega ^1\alpha _{ij}^2F_{ij}(\omega )𝑑\omega `$. The corresponding matrix elements are $`\mu _{\sigma \sigma }^{}=0.13`$, $`\mu _{\sigma \pi }^{}=0.042`$, $`\mu _{\pi \sigma }^{}=0.03`$, $`\mu _{\pi \pi }^{}=0.11`$ and $`\lambda _{ij}(\omega _m=\omega _n)`$ from Ref. Dolgov,  are $`\lambda _{\sigma \sigma }=1.017`$, $`\lambda _{\sigma \pi }=0.213`$, $`\lambda _{\pi \sigma }=0.155`$, $`\lambda _{\pi \pi }=0.448`$. Due to the interband coupling terms in Eqs. (3, 4) both gaps close at the same $`T_c`$. The resulting temperature dependencies of the energy gaps, $`\mathrm{\Delta }_i\left(T\right),`$ are plotted in the inset of Fig. 1 and it is found that $`\mathrm{\Delta }_\sigma \left(T=0\right)`$ = 7.09 meV and $`\mathrm{\Delta }_\pi \left(T=0\right)`$ = 2.70 meV, with the $`2\mathrm{\Delta }/T_c`$ ratios being equal to 4.18 and 1.59, respectively. For comparison, also the BCS curve is shown for $`T_c`$ = 39.4 K. The BCS value for the gap that corresponds to $`T_c`$ = 39.4 K is 6.0 meV at 0 K. It can be seen that the temperature dependencies are qualitatively different from the BCS temperature dependence. The ratio of the gaps $`\mathrm{\Delta }_\sigma /\mathrm{\Delta }_\pi `$ increases for increasing temperatures, as was experimentally observed for example in Ref. Klein, .

The influence of impurities can be incorporated into the model. Intraband scattering does not change the two gaps (Anderson’s theorem), while the interband scattering can be included by terms $`\gamma _{ij}\mathrm{\Delta }_j/\sqrt{\omega _n^2+\mathrm{\Delta }_j^2}`$, $`\gamma _{ij}\omega _n/\sqrt{\omega _n^2+\mathrm{\Delta }_j^2}`$ in the Eliashberg equations (3, 4) respectively. The smallness of $`\gamma _{ij}`$ compared to $`\pi T_c`$ indicates that the double-gap feature should experimentally be observable, also in thin films, even for a certain amount of impurity scattering. A large amount of impurity scattering ($`\gamma _{ij}`$ exceeding the maximum phonon frequency) will cause the gaps to converge to the same value. From Eqs. (3, 4) and including the scattering terms an asymptotic value of $`\mathrm{\Delta }_\sigma =\mathrm{\Delta }_\pi =`$ 4.1 meV and $`T_c=`$25.4 K is found, giving a $`2\mathrm{\Delta }/k_BT_c`$ ratio of 3.7.

In order to obtain the normal state resistance, we have to evaluate the effective junction transparency components $`D_{ij}`$. In the case of a specular barrier, $`U(x)=U_0\delta (xx_0)`$, $`D_{ij}`$ is given by

$$D_{ij}=\frac{v_{n,i}v_{n,j}}{\frac{1}{4}\left(v_{n,i}+v_{n,j}\right)^2+U_0^2/\mathrm{}^2}.$$
(5)
It follows from Eqs. (2, 5) and as first pointed out in Ref. Mazin, , that the normal state conductance $`R_{ij}^1`$ in the large $`U_0`$ limit is proportional to the Fermi-surface average $`Nv^2_i`$. The latter is proportional to the contribution of the electrons in band $`i`$ to the squared plasma frequency $`\left(\omega _p^i\right)^2`$. This essentially simplifies the task of summing up the interband currents since the partial plasma frequencies are available from the band structure calculations Kortus ; An ; Kong ; Bohnen . The normal state junction conductance is thus proportional to $`\left(\omega _p^i\right)_{}^2v_n_j`$, where $`v_n_j`$ is the average Fermi velocity projection in the corresponding band (see Table 1). In order to sum up the contributions of different bands, we restrict ourselves to the weighing factors $`\left(\omega _p^i\right)^2`$, neglecting the difference in $`v_n_j`$. This is a reasonable approximation since the difference between $`v_F`$ in the $`\sigma `$ and $`\pi `$ bands in the $`a`$-$`b`$ plane is rather small, while for $`c`$-axis tunneling only the $`\pi `$ band contributes, as will be shown later, so that the problem of summation does not appear in this case.

SIN tunneling. The conductance in a MgB2-I-N tunnel junction is the sum of the contributions of two bands. Each of the conductances is given by the BTK model BTK , where the corresponding normal state conductances $`R_{Ni}^1`$ are proportional to the minimum of the square of the plasma frequencies at the N and MgB2 sides. Since the plasma frequency in a typical normal metal (e.g. Au, Ag) is larger than the plasma frequencies in MgB2, the conductances are limited by the electrons on the MgB2 side

$$R_{N\sigma }^1/R_{N\pi }^1=\left(\omega _p^\sigma \right)^2/\left(\omega _p^\pi \right)^2.$$
(6)
Finally, the normalized conductance of an N-I-MgB2 contact is given by

$$\sigma \left(V\right)\frac{\left(\frac{dI}{dV}\right)_{NIS}}{\left(\frac{dI}{dV}\right)_{NIN}}=\frac{(\omega _p^\pi )^2\sigma _\pi \left(V\right)+(\omega _p^\sigma )^2\sigma _\sigma \left(V\right)}{(\omega _p^\sigma )^2+(\omega _p^\pi )^2}.$$
(7)
Here, the dimensionless conductances $`\sigma _{\sigma ,\pi }\left(V\right)`$ are provided by the BTK model, with the calculated values for the gaps and plasma frequencies, as shown in Table 1.

In the conductance versus voltage plot (Fig. 1) for tunneling in the $`a`$-$`b`$ direction, two peaks are clearly visible, in qualitative agreement with the experimental data Klein ; Bugoslavsky ; Szabo . The ratio of peak magnitudes is not only determined by the ratio of the plasma frequencies, but also by thermal rounding and by the barrier strength $`Z_{\text{BTK}}=U_0/\mathrm{}v_F`$ (where $`v_F`$ is taken constant for the different bands for the same reason as was given in the determination of $`R_{ij}`$). In particular, the peak at the smaller gap dominates in the small $`Z_{\text{BTK}}`$ regime (point-contact), while the second peak dominates at large values of $`Z_{\text{BTK}}`$ (tunneling), as may be seen in Fig. 1 at 4.2 K.

Due to the smallness of $`\omega _p^\sigma `$ in the $`c`$-direction, it can be seen from Eq. (7) that the conductance in the $`c`$-direction is only determined by the $`\pi `$ band. In this case, no double-peak structure is expected in the conductance spectrum. This explains, together with the dependence on $`Z_{\text{BTK}}`$, why in some experiments only one peak was observed Gonelli  or why the second peak was weak Klein ; Bugoslavsky ; Szabo .

Note, that the assumption of the ratio of the normal state conductivities being equal to the ratio of the square of the plasma frequencies holds when the interface is a $`\delta `$-function shaped tunnel barrier, with large $`Z_{\text{BTK}}`$. This means that for small $`Z_{\text{BTK}}`$, the results should be considered as a qualitative indication only. In the latter case, as well as for other types of barriers, a numerical integration of Eqs. (2, 5) must be performed.

SIS Josephson tunneling. We consider Josephson tunneling between two MgB2 superconductors. With the values for the plasma frequencies, $`\omega _p^\sigma <\omega _p^\pi `$, this gives $`R_{\sigma \pi }=R_{\pi \sigma }=\mathrm{max}(R_{\sigma \sigma },R_{\pi \pi })=R_{\sigma \sigma }`$ and $`R_{\sigma \pi }/R_{\pi \pi }=R_{\sigma \sigma }/R_{\pi \pi }=(\omega _p^\pi /\omega _p^\sigma )^2>1`$. The total conductance is given by $`R_N^1=_{ij}R_{ij}^1`$.

For tunneling in the $`a`$-$`b`$ plane (as can be realized for example in an edge configuration), with $`R_{\sigma \pi }=R_{\pi \sigma }`$ and $`I_{\sigma \pi }=I_{\pi \sigma }`$, the total $`I_cR_N`$ product becomes

$`I_cR_N={\displaystyle \frac{I_{\sigma \sigma }R_{\sigma \sigma }+2I_{\sigma \pi }R_{\sigma \pi }+I_{\pi \pi }R_{\pi \pi }\left(\omega _p^\pi /\omega _p^\sigma \right)^2}{3+\left(\omega _p^\pi /\omega _p^\sigma \right)^2}}.`$ (8)
The results of numerical calculations are presented in Fig. 2. Due to strong-coupling and interband coupling effects, the temperature dependencies of $`I_{ij}R_{ij}`$ differ from the well-known Ambegaokar-Baratoff result for an SIS junction between isotropic superconductors, most clearly demonstrated by the positive curvature of the $`I_{\pi \pi }R_{\pi \pi }`$ contribution. The $`I_cR_N`$ value at $`T`$ = 4.2 K is 5.9 mV.

For tunneling along the c-axis, the only contribution to the $`I_cR_N`$ product comes from $`I_{\pi \pi }R_{\pi \pi }`$, because of the negligible value for $`\omega _p^\sigma `$ in the c-axis direction. This gives $`I_cR_N`$ = 4.0 mV at $`T`$ = 4.2 K.

The plasma frequency in a certain direction, under an angle $`\phi `$ with the $`a`$-$`b`$ plane, can be determined from the ellipsoid equation $`(\omega _p^i)^2=(\omega _{p,x}^i)^2+(\omega _{p,z}^i)^2,`$ and $`(\omega _{p,z}^i/\omega _{p,c}^i)^2+(\omega _{p,x}^i/\omega _{p,ab}^i)^2=1`$, where $`\omega _{p,x}^i`$ and $`\omega _{p,z}^i`$ form the decomposition of $`\omega _p^i`$. Because of the negligible value of $`\omega _{p,c}^\sigma `$, it is evident that $`\omega _p^\sigma `$ is negligible for nonzero values of $`\phi =\mathrm{arctan}(\omega _{p,z}^i/\omega _{p,x}^i)`$. This implies that tunneling under a nonzero angle with the $`a`$-$`b`$ plane gives the same result as tunneling in the $`c`$-axis direction, namely $`I_cR_N`$ = 4.0 mV at $`T`$ = 4.2 K. For angles approaching zero (of the order of 0.6), $`I_cR_N`$ rapidly increases towards the maximal value for tunneling from $`a`$-$`b`$ plane to $`a`$-$`b`$ plane, namely $`I_cR_N`$ = 5.9 mV at $`T`$ = 4.2 K. For a large amount of impurity scattering the $`I_{ij}R_{ij}`$ values converge to the same value. It follows in that case from Eq. (8), with the plasma frequencies from Table 1, that $`I_cR_N`$ becomes almost isotropic.

Finally, tunneling from MgB2 into a superconductor $`S^{}`$ with a single gap will be considered (we take Nb as an example). The resulting $`I_{iS^{}}R_{iS^{}}`$ temperature dependencies are calculated numerically, using 1.4 mV for the energy gap in Nb. The ratio of resistances is determined from Eq. (2). Since typical values of plasma frequencies in other superconductors are bigger than in MgB2 (e.g. 9.47 eV for Nb, 12.29 eV for Al and 14.93 eV for Pb, see Ref. Maksimov, ), the following expression is obtained

$$I_cR_N=\frac{I_{\sigma S^{}}R_{\sigma S^{}}+I_{\pi S^{}}R_{\pi S^{}}\left(\omega _p^\pi /\omega _p^\sigma \right)^2}{1+\left(\omega _p^\pi /\omega _p^\sigma \right)^2},$$
(9)
when tunneling occurs into the $`a`$-$`b`$ plane of the MgB2. In the case of $`c`$-axis tunneling, only the $`I_{\pi S^{}}R_{\pi S^{}}`$ contribution remains. The results for tunneling from Nb to MgB2 are also indicated in Fig. 2. Other superconductors give qualitatively similar results. The only scaling parameter is the critical temperature of the superconducting counter-electrode.

Our results for Josephson tunneling provide an upper bound for $`I_cR_N`$ products, being 5.9 mV and 4.0 mV for tunneling into the $`a`$-$`b`$ plane and $`c`$ direction respectively. There have already been several observations of Josephson currents in MgB2 junctions Josephson , with $`I_cR_N`$ values that are much lower than our predictions. This can be due to extrinsic reasons such as a degradation of the $`T_c`$ of surface layers in the vicinity of the barrier, the barrier nature and barrier quality. From our model, however, it follows that polycrystallinity does not reduce the Josephson coupling very much, as indicated by the calculated value of $`I_cR_N`$ of 4.0 mV for c-axis transport, neither does strong impurity scattering because of the relatively large average gap of 4.1 meV in this case.

In conclusion, Josephson tunneling in MgB2-based junctions is discussed theoretically in the framework of a two-band model. The gap functions in different electronic bands are calculated using the Eliashberg formalism together with band structure information. This provides a basis to interpret electronic transport in MgB2. We have shown the possibility to observe either one or two gaps in point-contact spectra of MgB2, depending on the tunneling direction, barrier type and amount of impurities. The results are also relevant for the electronic application of MgB2 since they provide the limit for the Josephson coupling strength in MgB2 based junctions. For MgB2 in the clean limit we have shown that $`I_cR_N`$ values as high as 5.9 mV can be expected for MgB2 tunnel junctions if tunneling occurs in the direction of the $`a`$-$`b`$ plane. In other cases the limiting $`I_cR_N`$ values will not exceed 4.0 mV. Our predictions for the gap and $`I_cR_N`$ anisotropy and for the $`I_c`$ vs $`T`$ dependence in MgB2-based junctions can be verified experimentally and thus may help to settle the current debate on two-band superconductivity in MgB2.

The authors thank D.H.A. Blank, H. Hilgenkamp and I.I. Mazin for useful discussions. This work was supported by the Dutch Foundation for Research on Matter (FOM).