text:
# Electronic properties of metal induced gap states at insulator/metal interfaces — dependence on the alkali halide and the possibility of excitonic mechanism of superconductivity

## I Introduction

Heterointerfaces, namely solid-solid interfaces between very dissimilar materials, have been an issue of great interest in condensed matter physics, where various fascinating possibilities may be considered. Especially, the electronic structure of metal/insulator or metal/semiconductor interfaces has a long history of intensive studies dating back to the 1930’s.

One of the most crucial factors which govern the electronic structure of heterointerfaces is the metal-induced gap states (MIGS), whose eigenenergies lie in the band gap of the semiconductor (or insulator). While MIGS spread well into the metal side of the interface, they decay exponentially into the semiconductor. The decay length is an order of a few Å for typical semiconductors, where MIGS have been believed to determine the Schottky barrier height.

If we have an insulator/metal interface rather than a semiconductor/metal interface, on the other hand, it has generally been believed that the decay length of the MIGS into the insulator is negligible. This view has been based on quite a plausible model calculation, where the insulator/metal is represented by a simple tight-binding model/jellium. For a typical ionic insulator such as LiCl, the band gap ($`O(10)`$ eV) is an order of magnitude greater than the transfer energy ($`O(1)`$ eV) between the cations and anions, and the decay length is estimated to be smaller than 1 Å, which would suggest that MIGS should be irrelevant at insulator/metal interfaces.

Experimentally, it is difficult to detect MIGS at the heterointerface between wide gap insulators and metal, because clean interfaces are difficult to be fabricated. While MIGS have been reported in MgO/Cu(222) or MgO/Ag(001), the question of whether the observed gap states are formed by local chemical bonding remained, since MgO(222) is a polar surface so that a strong hybridization between the O 2p band and the Ag 5sp band is expected for MgO/Ag.

Recently, Kiguchi et al have succeeded in observing MIGS at LiCl/Cu and LiCl/Ag interfaces. They have grown a LiCl film on a metal substrate in a layer-by-layer fashion, and exploiting the near edge x-ray absorption fine structure (NEXAFS) which is an element-selective method, they have observed an extra peak below the bulk edge onset for LiCl film, which was interpreted as an evidence for MIGS. The reason why the detection eluded earlier studies is that, in conventional methods such as ultraviolet photoemission spectroscopy(UPS), inverse photoemission spectroscopy, or electron energy loss spectroscopy(EELS), the signal from the interface is obscured by significant contributions from the substrate. Kiguchi et al. have estimated the decay length of MIGS, which are $`2.6\pm 0.3`$ Å for LiCl/Cu(001) and $`2.9\pm 0.7`$ Å for LiCl/Ag(001), values much greater than those expected in simple models. An ab initio electronic structure calculation for LiCl/metal system has supported the existence of the MIGS that is strongly localized at the interface.

The purpose of the present paper is to study the electronic structure of interfaces between various alkali halides and a metal by means of first principles calculation, and to search for novel possibilities specific to metal/insulator interfaces. We shall show that (i) MIGS have a tail on halogen sites with a $`p_z`$-like character, whose penetration depth ($`\lambda `$) is as large as half the lattice constant of the bulk alkali halide, which are totally consistent with the experimental result, (ii) while $`\lambda `$ is insensitive to the carrier density in the metal, the penetration depth is dominated by the energy gap of the alkali halide. Specifically, the band gaps are $`E_{\mathrm{LiF}}>E_{\mathrm{LiCl}}>E_{\mathrm{LiI}}`$, and the the penetration depth has the expected inequality, $`\lambda _{\mathrm{LiF}}<\lambda _{\mathrm{LiCl}}<\lambda _{\mathrm{LiI}}`$.

We can consider various possibilities for many-body effects unique to heterointerfaces. There is in fact a long history for those, where the notable ones include metal-insulator transition and superconductivity. In the context of the latter, here we propose that large density of states arising from the MIGS inside the insulator is a good news for the superconductivity due to the exciton mechanism. The history of the exciton mechanism of superconductivity dates back to the proposals by Little, Ginzburg, and Bardeen et al. For this mechanism, the coexistence, in real space, of dispersionless excitons and metallic carriers is required. While this is difficult to achieve in general, if the density of states due to MIGS inside the insulator is large enough, we may envisage that the insulator/metal interface such as those discussed here may provide a possible ground for superconductivity. We shall discuss this toward the end of the paper.

## II method

The electronic structure of MIGS is studied with a first-principles calculation based on the density functional theory. We adopt the exchange-correlation functional introduced by Perdew, Burke and Wang and ultra-soft pseudo-potentials in a separable form. We have first calculated the bulk system and determined the lattices. The obtained values (against experimental ones in parentheses) are 4.09 (4.02) Å for LiF, 5.13 (5.13) Å for LiCl, 5.95 (6.00) Å for LiI and 5.63 (5.63) Å for NaCl.

We have then introduced a slab model, where we put a jellium in the middle of the unit cell sandwiched from top and bottom by an alkali halide layer (Fig.1). Starting from this structure as an initial state, we have performed structure optimization. To confirm that the simplified model with the jellium capture main features in heterostructures with real metals, we have also studied the system where five layers of Cu is considered in place of the jellium. The lattice mismatch between LiCl and Cu happens to be small, so that we can take a small unit cell. The cut-off energy of the plane-wave expansion for the wave functions are 20.25 Ry for LiF, LiCl and LiI/jellium, 25.00 Ry for NaCl/jellium, and 56.25 Ry for LiCl/Cu. The atomic configurations and the corresponding electronic states in the ground states are obtained with the conjugate gradient scheme.

The local density of states is calculated as $`_i|\varphi _i(x,y,z)|^2`$, where $`\varphi _i`$’s are eigenfunctions and the summation is taken over a window $`E_F<E_i<E_F+0.5\mathrm{eV}`$ with $`E_F`$ being the Fermi energy. We have set the window above $`E_F`$ to compare with the experiment, where the NEXAFS detects the unoccupied states.

## III Results

### A LiCl

Let us first discuss the case of LiCl/jellium. Since MIGS have been observed both in LiCl/Cu(commensurate) and LiCl/Ag(incommensurate), the commensurability in the lattice constants across the heterointerface should not be essential for the formation of MIGS, as we shall confirm below from the results for LiCl/Cu and LiCl/jellium. So we start with the jellium model.

We show in Fig.2(a) the band structure of LiCl/jellium where LiCl layers, each one monolayer thick, sandwiches the jellium that has $`r_s=2.5`$. We have also displayed the band structure for an isolated LiCl monolayer. From the comparison we can see that the band structure of LiCl is not severely affected by the presence of the jellium.

It is well-known that GGA or LDA generally underestimate the band gap, while GW approximation improves this. For instance, the band gap for bulk LiCl is 9.2 eV in GW, which is very close to the experimental value of 9.4 eV. However, the band-gap underestimation may be amended via the self-energy correction (i.e., without any corrections for LDA wave functions), the shape of the wave functions should be reliable even in LDA, hence the present calculation for LDOS around the $`E_F`$ is expected to be a good approximation.

In Fig.2(b,c) we have plotted the local density of states (LDOS) integrated along $`x`$ or $`y`$-directions. We can see that LDOS has a significant ($``$ lattice constant) tail with a $`p_z`$ character on the halogen site in addition to a usual evanescent-wave like part that exponentially decays away from the interface. In the experiment, while there is no information on Li sites since the element-selective NEXAFS was tuned for Cl in that experiment, an unoccupied state on Cl sites is observed near $`E_F`$. They have a $`p_z`$ character, and spread with a decay length of a few ÅṪhese properties are totally consistent with the present results.

To confirm that the penetration depth of the tail of MIGS is as large as a few Å, we have also calculated a larger system, where LiCl double layers sandwich the jellium. From the band structure and the LDOS in Fig.3 we can see that the MIGS do not penetrate into the second layer of LiCl.

If we now replace the jellium with (five layers of) Cu in Fig.4, we can see that the LDOS for the LiCl/Cu has a long tail just as in the result for the jellium model, so that the characteristic feature of MIGS should not be an artifact of the jellium model.

Why MIGS have amplitudes only on Li sites can be explained as follows. In general, MIGS can be divided into two groups: one having a conduction-band character and the other having a valence-band character in the insulator side of the interface. Noguera et al have introduced the energy $`E_{\mathrm{ZCP}}`$, above which MIGS have a valence-band character. While the charge neutrality in the insulator side is satisfied when $`E_{\mathrm{ZCP}}`$ is equal to $`E_F`$, charge transfer from the insulator to the metal occurs if $`E_{\mathrm{ZCP}}>E_F`$. In Fig.5, we plot the difference in the valence-charge density $`\delta \rho `$ between the LiCl+jellium and the isolated LiCl. We can see that, while $`\delta \rho `$ is positive for most of the insulator region ($`0<z<2`$), it is negative around the Cl site. Thus $`E_{\mathrm{ZCP}}`$ should be higher than $`E_F`$ in this system, and MIGS should have a tail on the anion sites.

Let us move on to the $`r_s`$ dependence. In Fig.6, we show the result when $`r_s`$ is increased from $`2.5`$ to $`6.0`$. We can see that, although LDOS at the interface is smaller than that for $`r_s=2.5`$, the decay length for $`r_s=6`$ is similar to that for $`r_s=2.5`$. Thus we may expect that penetration depth does not sensitively depend on the nature of the metal.

### B Other alkali halides

Let us move on to the cases of other alkali halides. We have performed calculations for LiF, LiI and NaCl. The results are shown in Fig.7 for LiF, Fig.8 for LiI, and Fig.9 for NaCl. The results are similar to that for LiCl(Fig.2). Namely, MIGS have a long tail on anion sites with a penetration depth $`\lambda `$ as large as the layer spacing $`a_z=a/2`$. The lattice constants are $`a=4.09`$ Å for LiF, 5.13 Å (LiCl), 5.95 Å (LiI) and 5.63 Å (NaCl), and $`\lambda `$ for LiI is indeed the largest (1.5 times as large as that for LiF). Note that anions with higher electron affinity have smaller ionic radii in general, so that alkali halides with largeer band gap have smaller lattice constants.

## IV Possibility of exciton mediated pairing

One of the most fascinating possibilities which should be explored for heterointerfaces is, we believe, the exciton-mediated superconductivity. The exciton mechanism has been originally proposed by Little for metallic, quasi-one dimensional spines of conducting electrons to which dielectric organic molecules are conceived to be attached chemically. Subsequently, two-dimensional version of this mechanism has been proposed, where a metallic film covered by a dielectric coating was considered.

To realize the exciton mechanism in general, carriers have to coexist and strongly interact with the excitons. Thus the penetration length must be sufficiently large, which is why it has been believed that the band gap of the coating material should not be too large. On the other hand, Inkson and Anderson have pointed out that, if we adopt the model dielectric function for typical semiconductors, exciton-electron coupling becomes small, partially because the exciton has a large energy dispersion so that the phase space over which the mechanism works for the formation of paring is restricted. In fact, in Ref. the energy dispersion of the exciton is ignored. In other words, they have considered a Frenkel exciton, which is usually unrealistic for (narrow gap) semiconductors.

Now, the present result for MIGS at metal/insulator interfaces can reconcile the very dilemma. Namely, while a conventional wisdom tells that the penetration depth should be extremely small for metal/wide-gap insulator interfaces, the present result suggests that the depth is of the order of the lattice constant for the MIGS, so that a substantial portion of the carriers can interact with the exciton in the insulator side. Thus it becomes an intriguing problem to estimate the electron-exciton coupling for the alkali halide/metal system, so we do this following the argument of ref..

According to ref., the exciton-electron coupling $`\lambda _{\mathrm{ex}}`$ is estimated as $`s\gamma ^2(D/L)\mu (\omega _p/\omega _g)^2`$, where $`s`$ is the screening factor for the Coulomb interaction, $`\gamma `$ the ratio between the LDOS at the interface and that for the bulk, $`D`$ the average depth of penetration, $`L`$ the thickness of the metal film, $`\mu `$ the screened Coulomb interaction averaged over the Fermi surface multiplied by the density of states at $`E_F`$, $`\omega _p`$ the electronic plasma frequency, and $`\omega _g`$ the band gap of the insulator. If we assume that the dielectric function can be approximated as

$`\epsilon (\omega )=1{\displaystyle \frac{(\epsilon _{\mathrm{optical}}1)\omega _0^2}{\omega ^2\omega _0^2}},`$

where $`\omega _0\omega _g`$ is the typical exciton energy, and $`\epsilon _{\mathrm{optical}}`$ the optical dielectric constant, then $`\lambda _{\mathrm{ex}}`$ becomes $`s\gamma ^2(D/L)\mu (\epsilon _{\mathrm{optical}}1)`$.

Let us first look into the factor $`\gamma D/L`$. This factor is, in the simplified model due to Bardeen, the fraction of the time over which the metal electrons spend in the insulator. In the more realistic model we have at hand, we can more precisely integrate the LDOS in the insulator region ($`\rho _i`$, shaded region in Fig.2(c) etc) since MIGS do not possess a monotonic (exponential) LDOS, and we can then divide the quantity by that for the metal region ($`\rho _m`$) to obtain the factor corresponding to $`\gamma D/L`$. Namely, the factor should read $`b2\rho _i/\rho _m`$, in place of $`\gamma D/L`$, where $`a_z=a/2`$ (see Fig.1). The present model is a sandwich structure where the metallic layer is covered by the dielectric coatings from both sides, so that the factor 2 appears in r.h.s. Here we set that the thickness of metallic layer to be 10 Å. The resulting values are $`b=0.08`$ (for LiF), 0.18(LiCl), 0.16(LiI), and 0.12(NaCl).

On the other hand, $`(\epsilon _{\mathrm{optical}}1)`$ are 0.92, 1.68, 2.4 and 1.33 for LiF, LiCl, LiI and NaCl, respectively. Therefore, if we adopt the same values for $`s(=1/2)`$, $`\mu (=1/3)`$, $`\gamma (=1/2)`$ as in Ref., $`\lambda _{\mathrm{ex}}`$ can become as large as 0.06 for LiI, while we have $`\lambda _{\mathrm{ex}}=0.01`$(LiF), 0.05( LiCl), and 0.03(NaCl).

If we adopt weak-coupling conventional superconductors, e.g. Al, as the metal, the superconducting transition temperature can be roughly estimated as

$`T_c=\theta \mathrm{exp}\left({\displaystyle \frac{1}{(\lambda _{\mathrm{ph}}+\lambda _{\mathrm{ex}}\mu ^{})}}\right),`$

where $`\lambda _{\mathrm{ph}}`$ is the electron-phonon coupling, $`\lambda _{\mathrm{ex}}`$ the electron-exciton coupling, and $`\mu ^{}`$ the renormalized $`\mu `$. So $`\theta `$ is roughly expressed as

$`\mathrm{log}\theta {\displaystyle (\mathrm{log}\omega )\left[\rho _{\mathrm{ex}}(\omega )+\rho _{\mathrm{ph}}(\omega )\right]\frac{d\omega }{\omega }},`$

where $`\rho _{\mathrm{ph}}`$ ($`\rho _{\mathrm{ex}}`$) is the density of phonon (exciton) states. In general, $`\rho _{\mathrm{ph}}`$ has a peak in a low energy region, so that, when phonons and excitons coexist, the contribution of the phonon part to $`\theta `$ should be dominant than that of the exciton part. Thus $`\theta `$ basically scales with that of phonon, but the exciton contribution on top of that can enhance the $`T_c`$ significantly. For example, when $`\lambda _{\mathrm{ph}}0.4`$, $`\mu ^{}0.05`$, and $`\lambda _{\mathrm{ex}}0.1`$, $`T_c`$ becomes 1.6 times greater than that in the bulk. Thus we conclude that heterointerface between metal and LiI or LiCl can be a promising ground for detecting the existence of exciton-enhanced superconductivity.

In addition, the proximity of different systems enables us to expect entirely different mechanism. Namely, two of the present authors have earlier proposed that a class of repulsively interacting systems that consist of a carrier band and an insulating band can become superconducting, where the system can effectively mapped to an attractive Hubbard model. The point is, while the conventional boson-exchange pairing arises due to a retarded attraction, superconductivity in the carrier-insulator model occurs due to a nonretarded attraction. In their model, the metallic band is assumed to interact only with the anion sites (or to interact only with the cation sites). At the interface between an alkali halide and a metal, MIGS have their amplitudes on anion sites, so that the MIGS should interact primarily with anion sites, which just corresponds to the situation considered in Ref.. The detailed study on the combined effect of the nonretarded attraction (which decreases $`\mu ^{}`$) and the exciton effect (which increases $`\lambda `$ j is an interesting future problem.

## V Summary

We have studied the electronic properties of metal-induced gap states(MIGS) at interfaces between various alkali halides and a metal. We have performed first-principles density functional calculation to show that (i) in addition to a usual evanescent-wave component, MIGS have a tail on halogen sites with $`p_z`$-like character, whose penetration depth is as large as the lattice constant of the alkali halide, (ii) while $`\lambda `$ does not significantly depend on the carrier density of jellium, it depends on the energy gap of alkali halides, with $`\lambda _{\mathrm{LiF}}<\lambda _{\mathrm{LiCl}}<\lambda _{\mathrm{LiI}}`$. We have estimated the electron-exciton coupling, and have found that we may envisage that insulator/metal interfaces such as those discussed in the present study may provide a possible ground for exciton-mediated superconductivity.

## VI Acknowledgment

We would like to thank K. Saiki, M. Kiguchi and D. R. Jennison for fruitful discussions. The GGA calculation was performed with TAPP (Tokyo Ab-inito Program Package), where RA and YT would like to thank Y. Suwa for technical advices. Numerical calculations were performed on SR8000 in ISSP, University of Tokyo. This research was partially supported by the Ministry of Education, Science, Sports and Culture, Grant-in-Aid for Creative Scientific Research.


format: text

text:
# Prisoners’ dilemma in real-world acquaintance networks: Spikes and quasi-equilibria induced by the interplay between structure and dynamics

## I Introduction

The prisoners’ dilemma (PD) game is a powerful metaphor for the situation where mutual trust and cooperation is beneficial in a long perspective but egoism and guile can produce big short-term profit. Applications range from evolution of RNA viruses rna , through westernization of central Africa africa , to the social situation of soldiers in trench warfare axe:evo  and NASCAR drivers nascar . One of the major achievements has been to establish the condition for cooperative behavior when two players meet repeatedly axe:evo . Another direction has been to investigate the criteria for cooperation to be stable in social space. In this approach the spatial PD game by Nowak and May nowakmay  has been the basic model; a model where there are a number $`N`$ of players interacting only with players in its immediate surrounding. Traditionally, spatial games have been studied on regular lattices; but real social networks are very complexly organized wf —being random to some extent, but also having structure reflecting the social forces. As a step in this direction people have studied the spatial PD on more realistic model networks abramson ; our:pd . In this Communication we pursue this idea to the end and let the players’ encounters follow the ties of empirical social networks. The outcome is a very complex behavior arising from the interplay between the PD dynamics and the underlying network structure.

We represent the underlying networks as graphs $`G=(V,E)`$ where $`V`$ is the set of $`N`$ players (or vertices) and $`E`$ is the set of $`M`$ ties (or edges—unordered pairs of vertices) between them. The networks we use are mostly acquired from online interaction—through contacts within two Internet communites (pussokram.com and nioki.com) and through email exchange. (A summary of the used networks can be seen in Table 1.) Even if the structure of online interaction network differs from regular acquaintance networks, we believe that our results will hold for a quite large class of social networks. As a small test network we also use an acquaintance data constructed from a field survey (the “karate club” network).

In Nowak and May’s spatial PD game nowakmay  each player, at each time step, adopts one of two strategies: cooperate $`C`$ or defect $`D`$. To catch the dilemma one lets an encounter between two cooperators result in unity gain for both players, whereas a $`D`$-$`D`$ encounter gives zero gain for both. However, if a cooperator meets a defector the cooperator scores zero and the defector scores $`b(1,2)`$ ($`b`$ is called ‘temptation’). The gain $`g(v)`$ of a player $`v`$ is the sum of the gains from the encounters with its neighbors. If a neighbor scores higher than a player, the player follows the high scoring neighbor. But, to drive the system (and model occasional irrational moves) a player chooses the opposite strategy with a probability $`p_m`$. We use synchronous updating of the players, i.e. one time step of the algorithm consists of one sweep (over all the players) to calculate the individual gains, and one sweep to update the strategies.

## II Evolution of the cooperator level

Our key quantity is the cooperator density $`\rho `$—the fraction of players adopting the strategy $`C`$. In regular networks the cooperator level is characterized by random fluctuations and high-frequency oscillations (mostly of period two, so called ‘blinkers’ nowakmay ). But, as seen in Fig. 1, when a regular underlying network is replaced by empirical social networks we get a complex behavior of both upward and downward spikes and steps between a number of quasi-equilibrium levels. Other observations of quasi-stable states in games have been obtained for more complex update rules and larger strategy space nowak:wsls ; brauchli ; egui:evo , or for some special types of directed networks our:pd .

The steps between quasi-equilibria in Fig. 1 suggest a scenario of relatively few stable states, where mutations on important vertices may cause a shift from one quasi-stable state to another. To get a better picture we plot a histogram over $`\rho `$ (Fig. 2) that indeed has distinct peaks. For a given $`b`$ the position and relative height of the peaks are $`p_m`$-independent (if $`p_m`$ increases, the noise level of the histograms gets higher) and thus forms a kind of fingerprint of the network. We note that the histograms obtained from one long run or many shorter runs (with different random seeds) are the same, this indicates that the system is self-averaging.

## III Characteristics of the quasi-stable states

Now we turn to the structure of the quasi-stable states. To start with a concrete example, in Fig. 3 we study the karate-club network at the intermediate $`\rho `$-level of Fig. 1(d). In the Fig. 3(a)-configuration 111We speak of a *configuration* as the set of strategies of the players at a given time, and a *state* as the set of configurations comprising a period (if the system would evolve without mutations). $`v_A`$ and $`v_B`$ scores highest among the border players (players with a neighbor of the opposite strategy). This makes many players, including the highly connected $`v_C`$, change strategy to $`C`$. In the Fig. 3(b)-configuration, the defector $`v_B`$ has enough cooperators in their neighborhood to score higher than $`v_A`$, and thus complete the cycle. We note that a crucial point in sustaining the stability is that $`v_A`$ and $`v_B`$ is not connected. We also note that if $`v_D`$ is mutated in the Fig. 3(a) configuration the system would loose its periodicity and stay in a constant configuration (until the next mutation). This suggests that the quasi-equilibria in Fig. 1 does not correspond only to one state, but a set of states that are all close in Hamming distance.

How does the high-$`\rho `$ states differ from the low-$`\rho `$ states? From the update rule to follow the player with locally highest gain, one expects that defectors gain comparatively much in low-$`\rho `$ and vice versa. This is indeed true as shown in Fig. 4 where we plot the average defector gain $`g_D`$ and cooperator gain $`g_C`$ along with $`\rho `$. One may therefore talk about a high-$`\rho `$ “cooperator controlled” and a low-$`\rho `$ “defector controlled” state. The average defector gain is directly related to the fraction of boundary ($`C`$-$`D`$) edges $`\varphi `$—the total defector gain is the sum of all boundary edges times $`b`$, which gives the following expression for $`\varphi `$:

$$\varphi =\frac{g_D(1\rho )N}{bM}.$$
(1)
The values of $`g_D`$ (seen in Fig. 4(c)) means that $`\mathrm{\hspace{0.17em}40}\%`$ of the edges are border between cooperators and defectors in the “defector controlled” state. This rather many considering that the biggest cooperator cluster consists of $`95\%`$ of all cooperators (the same figure for defectors is $`63\%`$). These observations suggests numerous situations where many defectors exploit high-degree cooperators, like $`v_C`$ is being exploited in Fig. 3(b). (This picture is confirmed under general circumstances our:forthcoming .)

## IV The steps in the cooperator level

To investigate the mechanism behind the steps of $`\rho (t)`$ we start by again considering the example of Fig. 3. If the most connected vertex $`v_B`$ (with degree 17) is mutated from $`D`$ to $`C`$ in the configuration of Fig. 3(a) then $`g(v_B)=4`$ while $`g(v_C)=5.6`$. In the next time step $`v_B`$ would shift to $`v_C`$’s strategy $`D`$ while $`v_C`$ follows the cooperator $`v_A`$, so in fact all vertices except one degree-2 vertex will have the same strategies as without a mutation on $`v_B`$. If on the other hand $`v_B`$ was mutated in the Fig. 3(b) configuration, all players would cooperate in the succeeding time step. By this example we conclude that a vertex’ importance with respect to mutation is not only dependent on its degree but also on the current configuration.

For a little more quantitative approach to this problem we study the result of a one-player mutation as follows: Starting at time $`t`$ we run the system without mutations until one period is completed—this is to start from a specific configuration $`c`$ of the local attractor state. Then we mutate vertex $`v`$ and let the system evolve (without mutations) until one period is completed. Finally we let the system complete yet a period and look for the configuration $`\widehat{c}`$ with the minimal Hamming distance to $`c`$ over the period (this is to make the difference zero if the system would relax back to the same state as prior to the mutation), and set $`\widehat{\rho }_v=\rho (\widehat{c})`$. In Fig. 5(a)-(d) we plot the ten $`\widehat{\rho }_v`$ that differs most from $`\rho (c)`$ at four different time steps in four different quasi-stable states. We observe that the system can change quasi-stable state due to mutations on single vertices, but at the same time not go to *any* other quasi-stable state. Furthermore, there is only a small number of mutations that actually causes a transition. Just as in the karate club example above the configurations comprising a state may vary much from one time-step to another. This fact explains the apparent lack of transitions to the $`\rho 0.35`$ state—in fact, the transition to this state (around $`t=\mathrm{8\hspace{0.17em}250}`$) occurs from a single-player mutation from a configuration in the $`\rho 0.65`$ peak. As expected, all vertices contributing to the lines in Fig. 5(a)-(d) have high degree (the lowest degree is $`111`$). But, as mentioned above, it is not true that a high degree implies a high importance—for example, of the ten most connected vertices, on average $`5.5`$ contributes to the ten most different $`\widehat{\rho }_v`$ values plotted in Fig. 5(a)-(d). It is thus not only the structure but the interplay between the structure and dynamics that gives rise to the quasi-stable equilibria. As a last point we note that the quasi-equilibria are stable in the absence of mutations—if the system at the time steps of Fig. 5(a)-(d) is evolved without mutation, the largest shift in $`\rho `$ is $`0.008`$.

## V Minimal structure causing steps and spikes

What features of the network structure are causing the nonstationary evolution of the cooperator density? As we have seen, the vertices have to have a highly skewed degree distribution—so that mutation on one vertex can cause many vertices, directly or indirectly, change strategy. Furthermore it should not be the case that the vertices of highest degree all are interconnected—then these would all follow the same strategy, and most of the other vertices too, so nothing more than small fluctuations of $`\rho `$, or steps involving almost all vertices, could occur. To test that there is a bias among the most connected vertices not to attach to each other, we sample different networks with the same degree sequence by randomly rewiring the original network with the restriction that the degree of each vertex is preserved roberts:mcmc ; maslov:pro . The average degree of the maximal subgraph of the ten vertices of highest degree, $`k_{10}`$, and the expectation value of the same quantity for random graphs with the same degree sequence as $`G`$, $`\overline{k}_{10}`$, for the four large networks are shown in Table 1. (The value ten is chosen arbitrarily—any number of the same order would do.) We see that $`k_{10}`$ is indeed significantly lower than $`\overline{k}_{10}`$ for all four networks. Online interaction represent a somewhat special type of acquaintance network, even if such networks are interesting in their own right well:comp , the question if this structure exists in ‘offline’ acquaintance networks is open (and will be discussed further in Ref. our:forthcoming ). We also note that if one runs the PD dynamics on a rewired network one sees spikes but no quasi-stable states our:forthcoming .

Networks with a broad degree distribution can be generated by, for example, Barabási and Albert’s (BA) scale-free network model ba:model : Start with $`m_0`$ disconnected vertices, and grow the network by adding a vertex $`v`$ with degree $`m`$ per time step. An edge from the new vertex are added to an old vertex $`w`$ with a probability proportional to $`w`$’s degree.

We try to model the behavior of the cooperator density by starting from BA model networks with $`m=m_0=3`$ and $`\mathrm{20\hspace{0.17em}000}`$ vertices. To tune how much mutually connected the vertices of highest degree are, we make the subset of the $`n_+`$ (or $`n_{}`$) vertices of highest degree complete by adding the missing edges (or completely disconnected by removing existing edges). As seen in Fig. 6 the steps and spikes of the real-world networks can be qualitatively reproduced. The resemblance is closest for removed edges among a small number of hubs Fig. 6(c). If the edges are removed within a larger set of hubs, the steps are increasing in frequency and number until the the structure is completely blurred out. If all missing edges are added within the top $`n_+`$ connected vertices, the steps vanishes but the spikes increase in size and width. An explanation is that when the vertices of highest degree are fully connected, a mutation of one of these would momentarily increase the number of cooperators by a large amount and also decrease the speed of restoration of cooperation. At the same time it is well-known that cooperation is promoted in highly-connected regions axe:socstr  which explains the high $`\rho `$-value of the quasi-stable state.

## VI Summary and conclusions

In summary, we have studied Nowak and May’s spatial PD game on empirical social networks. We find a complex time-evolution of the cooperator level characterized by large spikes and transitions between a number of quasi-stable states. The quasi-stable states of low $`\rho `$ are characterized by comparatively large gain for the defectors and large fraction of boundary edges, in quasi-stable states with high $`\rho `$ this situation is reversed. We find that mutations on high-degree players are likely to be responsible for transitions between the quasi-stable states; but the outcome of a mutation of a high-degree vertex is also much dependent on the current configuration. These findings are presented for some specific sets of parameter values, but they are qualitatively the same for a broad range of values. We argue that the structure causing this behavior is a inhomogeneous degree distribution, and that the number of edges within the subset highest-degree vertices is relatively low. Based on this observation we also construct model networks that reproduces the steps and spikes of the real-world networks. To epitomize, we believe that we have illustrated how underlying acquaintance patterns can give rise to a complex time-development of social instability—a picture that could be empirically testable by carefully arranged social observations.

## Acknowledgements

The authors are grateful for comments from C. Edling, F. Liljeros and K. Sneppen; and acknowledges support from Swedish Research Council through Contract No. 2002-4135 and (B.J.K.) the Korea Science and Engineering Foundation through Grant No. R14-2002-062-01000-0.


text:
# Full counting statistics of multiple Andreev reflection

## Abstract

We derive the full counting statistics of charge transfer through a voltage biased superconducting junction. We find that for measurement times much longer than the inverse Josephson frequency, the counting statistics describes a correlated transfer of quanta of multiple electron charges, each quantum associated with the transfer of a single quasiparticle. An expression for the the counting statistics in terms of the quasiparticle scattering amplitudes is derived.

Due to the discrete nature of electric charge, the current in mesoscopic conductors generally fluctuates. Over the last decade, there has been an increasing interest, theoretical as well as experimental Buttrev , in the physics of current fluctuations. Most studies have been focused on noise, the second moment of the fluctuations, but recently a considerable interest has been shown for the full distribution of charge fluctuations, the full counting statistics (FCS) Levitov1 . A variety of theoretical approaches to the FCS, ranging from quantum mechanical Levitov2 ; Beenakker , via quasiclassical Nazarov  to classical Pilgram , have been developed. The third moment of current fluctuations was very recently measured Reulet , opening up the road, as well, to experimental investigation of the higher moments of the fluctuations.

Noninteracting electrons in purely normal conductors are transferred one by one Levitov1 . In normal-superconducting junctions, the charge transfer mechanism across the normal-superconducting interface, at energies below the superconducting gap, is Andreev reflection. As a consequence, the FCS include terms describing correlated transfer of pairs of electrons MuzKhmel . Recently, Belzig and Nazarov Belzig2  studied the FCS in superconducting junctions with a fixed phase difference between the superconducting electrodes. They found that the classical interpretation of the FCS, the probability to transfer a given number of electrons across the junction during the measurement, could imply negative probabilities. Coupling the junction to a detector Nazarov1 , they showed that this resulted from an attempt to interpret the phenomena of supercurrent with classical means.

In voltage biased superconducting junctions, the physical situation is quite different. Due to the applied voltage bias, the superconducting phase difference oscillates with the Josephson frequency $`2eV/\mathrm{}`$, giving rise to both dc and ac-components of the current. For measurement times much longer than the inverse Josephson frequency, only the dc-current, which is dissipative, contributes to the net charge-transport. Microscopically, the charge is transported between the two superconductors via coherent multiple Andreev reflections (MAR) OldMAR . The current has been studied in various junctions, both theoretically TheoryMAR  and experimentally ExpMAR . Recently, also the noise was studied Noisetheory ; Naveh ; Noiseexp .

In this paper we present the FCS of charge transfer through a voltage biased junction, in terms of the amplitudes for quasiparticle scattering. Each quasiparticle scattering process results in an integer number of electron charges being transferred across the junction. As a consequence the FCS can be interpreted in classical terms. At temperatures much lower than the superconducting energy gap $`\mathrm{\Delta }`$, many-quasiparticle scattering processes are exponentially suppressed, resulting in a simple probability distribution, containing only the probabilities for single quasiparticle scattering. This distribution reproduces known results for dc-current and zero-frequency noise. We discuss in detail the third cumulant for single channel junctions and diffusive junctions shorter than the superconducting coherence length.

We consider a superconducting junction consisting of two superconducting reservoirs connected via a normal, mesoscopic conductor (see Fig. 1). For simplicity of notation, we consider a junction with a single transport mode, the multi-mode generalization is discussed below. A voltage $`V`$ is applied between the two reservoirs.

The single-particle wavefunctions in the junction, solutions to the time dependent Bogoliubov-de Gennes equation, are scattering states labelled by the incoming quasiparticle type, injection energy and reservoir. The scattering states are superpositions of amplitudes for quasiparticles at energies $`\pm neV`$ from the injection energy, counted from the local chemical potential in each contact. The amplitude for an incoming quasiparticle of type $`\alpha `$ at energy $`E_m`$ to exit the junction as a quasiparticle of type $`\beta `$ at energy $`E_n=E_m+(nm)eV`$ is denoted $`s_{nm}^{\beta \alpha }`$ (see Fig. 1). This amplitude is a function of the scattering matrix of the normal conductor and the Andreev reflection amplitudes Johansson .

To access the FCS, we make the first important observation that quasiparticle scattering in a voltage biased superconducting junction is formally identical to scattering in a normal voltage biased junction, with an applied harmonic ac-field. The FCS in such a system was investigated in detail by Ivanov and Levitov. Following Ref.  we note that for measurement times much longer than the quasiparticle scattering time, the inelastic single mode scattering problem can be mapped onto an elastic scattering problem with many modes. The scattering between all energies and quasiparticle types of a “ladder” (see Fig. 1) is correlated, while different ladders contribute incoherently. We denote the ladder by its energy $`E_0`$ in the interval $`\mathrm{\Delta }2eVE_0<\mathrm{\Delta }`$ counting its “leg” in the left lead (e.g. $`E_0`$ in Fig 1). Thus we may concentrate on the scattering matrix $`S`$, with elements $`s_{nm}^{\alpha \beta }`$, of a single ladder, and then integrate over the ladder energy $`E_0`$. In general $`S`$ has infinite dimensions, but in our case the vanishing probability of Andreev reflection far outside the gap naturally cuts the number of relevant modes to the order of $`\mathrm{\Delta }/eV`$. Quasiparticle current (but not charge current) is conserved in the scattering processes at the normal-superconductor interfaces. As a consequence $`S`$ is unitary comment .

We are then, in line with Levitov2 ; LevitovTD , able to directly write down the characteristic function in terms of all different many-particle scattering probabilities

$$\chi _{E_0}^{qp}(\mathrm{\Lambda })=\underset{𝐢,𝐨}{}e^{i\left(_{n\beta 𝐨}\lambda _{n\beta }_{m\alpha 𝐢}\lambda _{m\alpha }\right)}P_{𝐢|𝐨},$$
(1)
where $`\mathrm{\Lambda }`$ is the set of counting fields $`\lambda _{n\alpha }`$, one for each mode $`n\alpha `$. The outer sum runs over all possible sets of incoming modes $`𝐢=\{m_1\alpha _1,m_2\alpha _2,\mathrm{}\}`$ and outgoing modes $`𝐨=\{n_1\beta _1,n_2\beta _2,\mathrm{}\}`$.

The many-particle scattering probabilities are given by

$$P_{𝐢|𝐨}=|s_𝐢^𝐨|^2\underset{m𝐢}{}f(E_m)\underset{m^{}𝐢}{}[1f(E_m^{})],$$
(2)
where $`f(E)=(1+\text{exp}[E/kT])`$ is the Fermi distribution function and $`|s_𝐢^𝐨|`$ is the determinant of the matrix formed by taking the columns $`𝐢`$ and rows $`𝐨`$ of $`S`$.

Eqs. (1) and (2) gives us the FCS for quasiparticle transfer in a voltage biased superconducting junction. The object of main interest is however the FCS of the charge transfer. We then make the second important observation that for measuring times much longer than the inverse Josephson frequency, the net charge transfer is directly related to the quasiparticle transfer. A quasiparticle of type $`\alpha `$ incident at energy $`E_m`$, which is scattered into an outgoing quasiparticle of type $`\beta `$ at energy $`E_n`$, transports exactly $`mn`$ electrical charges across the junction (for $`m<n`$, the transported charge is thus negative for the bias in Fig. 1). This is independent of the type of quasiparticles $`\alpha `$ and $`\beta `$, and also of the charge of the quasiparticles. This can be shown by studying any quasiparticle scattering path (see Fig. 1), keeping in mind that the process of Andreev reflection transfers exactly two electrons across the normal-superconductor interface, while a normal transmission transfers exactly one electron.

This can also be seen from energy conservation: An electron traversing the normal part of the junction from left to right will absorb the energy quantum $`eV`$ from the electric field, while an electron moving in the opposite direction will emit the quantum $`eV`$. The effective number of quanta $`eV`$ absorbed in a quasiparticle scattering process scattered thus equals the number of electrons transferred from left to right. We emphasize that this approach correctly counts all the electrons transferred, including the electrons entering the superconductor as Cooper pairs at energies within the gap.

We can thus count the transferred electrons by counting the transferred quasiparticles weighted by the number of electrons transferred in each quasiparticle scattering event. Thus by chosing the counting fields $`\lambda _{n\alpha }=n\lambda `$ in Eq. (1) we can directly write down the characteristic function for charge transfer

$$\chi _{E_0}(\lambda )=\underset{𝐢,𝐨}{}e^{i\lambda \left(_{n\beta 𝐨}n_{m\alpha 𝐢}m\right)}P_{𝐢|𝐨},$$
(3)
Following Ref Levitov2 , the characteristic function can be written on a determinant form as

$$\chi _{E_0}(\lambda )=\text{det}\left[1\overline{f}+\overline{f}S_\lambda ^1S_\lambda \right],$$
(4)
where the elements of the matrix $`S_\lambda `$ are given by $`(S_\lambda )_{mn}^{\alpha \beta }=s_{mn}^{\alpha \beta }e^{i\lambda (nm)/2}`$ and the diagonal matrix $`\overline{f}`$ has elements $`\overline{f}_{nn}^{\alpha \alpha }=f(E_n)`$. This results provides a general solution to the temperature and voltage dependence of the full counting statistics of voltage biased superconducting junctions, for measurement times much longer than inverse Josephson frequency, and the quasiparticle scattering time, as defined by the energy dependence of the scattering matrix $`\mathrm{}_E\mathrm{ln}|s_{mn}^{\alpha \beta }|`$. Eqs. (3) and (4) shows that the charge is transferred in correlated quanta of multiple electron charges.

At low temperatures $`kT\mathrm{\Delta }`$ all quasiparticle states below the gap are filled, while all states above the gap are empty. This simplifies Eq. (2) significantly since it fixes the set of incoming modes $`𝐢`$ to all modes below the gap, but still all different sets of outgoing modes should be considered.

The MAR ladder forms a single mode for transport in energy space Johansson . The scattering amplitudes $`s_{nm}^{\beta \alpha }`$ can be decomposed into amplitudes for entering the normal region, $`t_{m\alpha }^{m\pm }`$, propagation along the MAR ladder, $`t_{mn}^/`$, and leaving the normal region, $`t_{n\pm }^{n\alpha }`$, as

$$s_{nm}^{\beta \alpha }=t_{m\alpha }^{m+}t_{mn}^{}t_n^{n\beta },t_{mn}^{}=t_{mk}^{}t_{kn}^{}(m<k<n),$$
(5)
and analogously for $`E_n<E_m`$. Using this decomposition we find that choosing two or more outgoing modes above the gap in the set $`𝐨`$ gives $`P_{𝐢|𝐨}=0`$, since the matrix $`s_𝐢^𝐨`$ then contains two or more parallell rows. In other words, the many-particle scattering process of two or more quasiparticles passing the gap through the single mode in energy space is prohibited by the Pauli exclusion principle. For the outgoing sets with all outgoing modes below the gap $`|s_𝐢^𝐨|`$ is evaluated through making $`s_𝐢^𝐨`$ unitary by inserting a row and a column with the amplitudes $`t_{m\alpha }^{m+}t_{mg}^{}`$, where $`E_g`$ is the first energy above the gap in the ladder. Finally, the sets of one outgoing mode above the gap, using similar manipulations, give the single particle scattering probabilities. The characteristic function can then be expressed in single-particle scattering probabilities

$$\chi _{E_0}(\lambda )=\underset{m0,n>0}{}e^{i(nm)\lambda }\underset{\alpha ,\beta }{}|s_{nm}^{\beta \alpha }|^2.$$
(6)
This simple expression is the main result of this paper.

The cumulants of the charge transfer distribution function are obtained by taking derivatives of $`\mathrm{ln}\chi _{E_0}(\lambda )`$, and summing over all ladders (integrating over the energy $`E_0`$),

$$n(\tau )^m=\frac{\tau }{h}_\mathrm{\Delta }^{\mathrm{\Delta }2eV}𝑑E_0(i_\lambda )^m\mathrm{ln}\chi _{E_0}(\lambda )|_{\lambda =0},$$
(7)
where $`\tau `$ is the measurement time.

The integrands for the various moments are conveniently expressed in the physically relevant $`n`$-electron scattering probabilities Johansson ,

$$P_n(E_0)=\underset{m0}{}\mathrm{\Theta }(E_{m+n}\mathrm{\Delta })\underset{\alpha ,\beta }{}|s_{(m+n)m}^{\alpha \beta }(E_0)|^2(n1),$$
(8)
with $`P_0(E_0)`$, the probability of no quasiparticle passing the gap, giving the correct normalization $`_{n=0}^{\mathrm{}}P_n(E_0)=1`$. In Fig. 2 the probability of different $`n`$-electron scattering processes, as a function of voltage, are shown.

The integrands for the first three cumulants in Eq. (7) are the spectral current density $`I(E)=_{n=0}^{\mathrm{}}nP_n(E)`$, the noise spectral density $`S_I(E)=_{n=0}^{\mathrm{}}n^2P_n(E)I(E)^2`$, and the spectral density of the third cumulant $`C_3(E)=_{n=0}^{\mathrm{}}n^3P_n(E)I(E)(3S_I(E)+I(E)^2)`$. The expression for the current is identical to the one presented in Ref. Johansson , and the expression for the zero-frequency current noise reproduces the known results Noisetheory  (see the upper panel in Fig. 3).

The noise and the third cumulant, calculated numerically, are plotted in Fig. 3 for different transparencies of the normal contact. We see that the subgap structure is generally more pronounced in the third cumulant, compared to the current and noise. Furthermore $`C_3`$ changes sign both as a function of voltage and transparency.

In the tunnel limit $`D1`$ where all scattering probabilities in Eq. (6) are small, the corresponding probability distribution is Poissonian, i.e. describes an uncorrelated transfer of quantas of multiple charge. Within the $`n`$-th subgap region, i.e. $`n1<2\mathrm{\Delta }/eV<n`$, $`n`$-electron transfer dominate, giving $`S=2enI`$ and $`C_3=(2en)^2I`$, see the upper panel in Fig. 4. As seen in the right panel of Fig. 2 the picture is more complicated on the borders between the regions. Due to resonances up to three adjacent $`n`$-particle processes have similar strength Johansson , and the FCS is no longer Poissonian.

So far we considered only single mode junctions. In the limit of a short junction, when the scattering matrix of the normal conductor is independent on energy on the scale of $`\mathrm{\Delta }`$, it is possible, just as for the current and noise Naveh , to write the generating function in terms of the transmission eigenvalues $`D_n`$ of the normal conductor only. For all mesoscopic conductors where the transmission eigenvalue distribution is known, the FCS can thus be obtained via averaging the single mode result in Eq. (6). As an example, we present in the lower panel of Fig. 4 the first three moments for a short, diffusive junction. We note that the third moment is positive for all voltages. Just as the noise, the third moment has a divergency for $`eV0`$. However, in an experiment, we expect the fluctuations at small voltages to be strongly suppressed by inelastic scattering Gunsenheimer , effectively removing the divergency Naveh .

In conclusion, we have derived an expression for the full counting statistics of a voltage biased superconducting junction, describing a correlated transfer of quantas of multiple electron charge. It is found that the counting statistics can be expressed in terms of probabilities for quasiparticles to scatter between the two superconductors.

We acknowledge discussions with Yuriy Makhlin and Dimitri Bagrets. This work was supported by the Humboldt Foundation, the BMBF and the ZIP programme of the German government (GJ), by the Swiss network MANEP (PS), and by the Swedish Foundation STINT (ÅI).


text:
# BeppoSAX serendipitous discovery of the X-ray pulsar SAX J1802.7–2017

## 1 Introduction

High mass X-ray binaries (hereafter HMXBs) are young systems, and the neutron stars (NSs) that are often hosted in them usually have a strong magnetic field ($`B10^{12}`$ Gauss). Accretion onto these NSs occurs via capture of stellar wind matter and Roche-lobe overflow. HMXBs comprise two main subgroups: a) the supergiant systems, and b) the Be star systems. In the supergiant systems, the companion is of spectral type earlier than B2 and has evolved off the main sequence. Orbital periods are generally less than 10 days, orbits are circular, and the mass transfer takes place because of the strong stellar wind from the OB star and/or because of “incipient” Roche Lobe overflow. The Be star systems are characterized by emission lines (mainly the Balmer series) which originate in the equatorial circumstellar envelope of the companion star. The orbital periods in these systems tend to be longer than those of group a) and correlate well with the pulsar spin period (Corbet 1986). The orbits are usually moderately eccentric. Transient activity is common in Be HMXBs; different types of outburst have been observed from different sources, and occasionally also from the same source. Giant outbursts involve high peak luminosities, occur at any orbital phase and show only little (if any) orbital modulation of the X-ray flux. Recurrent outbursts usually involve lower peak luminosities, tend to occur close to periastron and sometimes show a strong modulation of the X-ray flux with the orbital phase (see e.g. Stella, White, & Rosner 1986). To date, about 80 HMXBs are known; $`40`$ of them show periodic X-ray pulsations with spin periods distributed over a wide range from 69 ms to $`24`$ min (see Nagase 1989 for a review).

While studying the atoll source GX 9+1 with BeppoSAX, we have discovered a new X-ray pulsator in a HMXB. We report here the results of the timing analysis of this new X-ray source.

## 2 Observations

The observation of the GX 9+1 field was carried out from 2001 September 16 02:01:30.0 (UTC) to 2001 September 20 03:00:08.5 (UTC), using the co-aligned Narrow Field Instruments (NFIs) on board BeppoSAX. These are: a Low Energy Concentrator Spectrometer (LECS; energy range 0.1–10 keV; Parmar et al. 1997), two Medium Energy Concentrator Spectrometers (MECS; energy range 1–10 keV; Boella et al. 1997), a High Pressure Gas Scintillation Proportional Counter (HPGSPC; energy range 7–60 keV; Manzo et al. 1997), and a Phoswich Detector System (PDS; energy range 13–200 keV; Frontera et al. 1997). The exposure times were $`60`$ ks, $`149`$ ks, $`142`$ ks, $`71`$ ks, for LECS, MECS, HPGSPC and PDS, respectively. The circular field of view (FOV) of the LECS and MECS is $`37^{^{}}`$ and $`56^{^{}}`$ in diameter, respectively, while those of the HPGSPC and PDS are hexagonal with FWHM of $`78^{^{}}`$ and $`66^{^{}}`$, respectively. The LECS and MECS detectors are position sensitive counters with imaging capability. The position reconstruction uncertainty for MECS is $`0.5^{^{}}`$ in the central area of $`9^{^{}}`$ radius, and $`1.5^{^{}}`$ in the outer region of the FOV (Boella et al. 1997). The HPGSPC and PDS systems do not have imaging capabilities, and their data are therefore difficult to interpret and analyse for individual sources when the FOV includes more than one source.

Figure 1 shows the BeppoSAX/MECS image (1-10 keV) centered at the position of the bright atoll source GX 9+1. A fainter source is visible at R.A.(2000.0)= 18h 02m 39.9s and Dec.(2000.0)= $`20^{}`$ $`17^{^{}}`$ $`13.50^{^{\prime \prime }}`$ (position uncertainty $`2^{^{}}`$), at an angular distance from GX 9+1 of $`22^{^{}}`$. The source was outside the FOV of the BeppoSAX/LECS. We have verified its presence in both the MECS2 and MECS3 images, which probably excludes that this was a ghost image of a source outside the MECS FOV. Moreover, we can be sure that the source was within the PDS FOV, because the source X-ray pulsations were also detected in the PDS data (see below). We searched for known X-ray sources in a circular region of $`30^{^{}}`$ centered at GX 9+1 in the SIMBAD data base. We found no known sources with a position compatible with that of the faint source; we therefore designate this serendipitous source as SAX J1802.7–2017.

We extracted the MECS events of SAX J1802.7–2017 from an elliptical region centered at R.A.(2000.0)= 18h 02m 39.9s and Dec.(2000.0)= $`20^{}17^{^{}}13.50^{^{\prime \prime }}`$ and the background events from an elliptical region similar to the one used for the source, centered at a symmetric position with respect to the center of the MECS FOV and not contaminated by GX 9+1 (see Fig. 1). In Figure 2 the lightcurves of SAX J1802.7–2017 plus background (crosses) and of the background (stars) are shown. The source was not detected during the first 50 ks when the count rate within the extraction region ($`0.034\pm 0.001`$ c/s) was compatible with the background count rate. In the following $`300`$ ks the source count rate showed a large variability with an average intensity significantly above the background level. Two flaring events took place $`110`$ ks and $`300`$ ks from the beginning of the observation.

## 3 Temporal Analysis

The arrival time of all events in the MECS and PDS were corrected to the solar system barycenter. We searched for periodicities by computing a Power Spectrum Density (PSD) in the range between $`4\times 10^4`$ Hz and 1 kHz from Fast Fourier Transforms (FFTs) performed on MECS data of SAX J1802.7–2017.

A main peak at a frequency of $`7.20`$ mHz, corresponding to a period of $`139`$ s, is evident (see Fig. 3). The second harmonic is also clearly visible. We performed a folding search for periods centered around $`139`$ s finding a not well defined period, because two $`\chi ^2`$ peaks were present at $`139.49`$ s and $`139.70`$ s, respectively.

To study possible delays in the pulse arrival times, or, equivalently, pulse period variations we divided the whole data set into 35 consecutive intervals each having a length of about 10 ks. We found that only 27 of these intervals had enough statistics to carry out the analysis described below. A folding search was performed in each interval for a range of trial periods centered around 139 s. The corresponding best periods were obtained by fitting the $`\chi ^2`$ versus trial period curve with a Gaussian function. The best periods varied significantly and with continuity between 139.44 s and 139.86 s with an average value of 139.608 s.

Phase delays (Fig. 4) were obtained by cross-correlating the folded lightcurves obtained for each of the 27 intervals with that of the whole observation used as a template (an average pulse period of 139.608 s was adopted in the folding). Initially we fitted the phase delays $`\mathrm{\Delta }\varphi `$ using

$$\mathrm{\Delta }\varphi =a_0+a_1t_n+a_2t_n^2,$$
(1)
where $`t_n`$ is the arrival time of the $`n`$th pulse. In this formula the linear term $`a_1=\mathrm{\Delta }P_{pulse}/P_{pulse}^2`$ and the quadratic term $`a_2=\dot{P}_{pulse}/2P_{pulse}^2`$ are related to a correction to the pulse period and the derivative of the pulse period, respectively. We obtained a $`\chi _{red}^2`$ of 2.6 indicating that the modulation of the phase delays cannot only be explained by the presence of a derivative of the pulse period. Moreover we find a value of $`\dot{P}_{pulse}/P_{pulse}`$ of $`0.40\pm 0.01`$ yr-1, more than an order of magnitude larger than the largest measured value in any known X-ray pulsar (i.e. $`10^2`$ yr-1 in GX 1+4, see Pereira et al. 1999). We checked whether the modulation of the phase delays could be explained by the propagation delays due to the orbital motion of the X-ray pulsar around a companion star by fitting the phases with

$$\mathrm{\Delta }\varphi =a_0+a_1t_n+B\mathrm{cos}\left[\frac{2\pi (t_nT_{\pi /2})}{P_{orb}}\right],$$
(2)
The results of the fit are reported in Table 1. The $`\chi _{red}^2`$ was 1.3. Since $`a_x\mathrm{sin}i=P_{pulse}B`$, we find $`a_x\mathrm{sin}i70`$ lt-s; the corrected pulse period was $`139.612`$ s.

We tried to add a quadratic term to expression (2) to take into account a spin period derivative, but, because of the small number of points and a relatively short observation, we could not constrain all parameters of the fit. In any case, by fixing the quadratic term in the range $`9\times 10^3`$$`9\times 10^3`$ days-2 (corresponding to a $`\dot{P}_{pulse}/P_{pulse}`$ in the range $`10^2`$$`10^2`$ yr-1) we found that the orbital parameters do not change significantly ($`2\%`$ changes in the orbital period value and $`4\%`$ in the $`a_x\mathrm{sin}i`$ value) with respect to the parameters reported in Table 1.

We corrected the arrival times of all the events, observed by the MECS and PDS instruments, to the center of mass of the binary system using the orbital parameters and the corrected pulse period reported in Table 1. We then re-computed, on the corrected folded lightcurves, the phase delays, and found that they are compatible with zero phase delay (see Fig. 5). This suggests that our data are compatible with a circular orbit, although an eccentricity could be present but not appreciated because of the relatively low statistics of our data. To estimate a rough upper limit to the eccentricity of the system, we fitted our phase delays substituting the circular orbital correction in equation (2) with a first-order approximation of the eccentric orbital correction (see e.g. van der Klis & Bonnet-Bidaud 1984), thus obtaining $`e0.2`$ (90% confidence level) and $`\chi _{\mathrm{red}}^21.4`$.

Using the estimated pulse period we folded the MECS lightcurves in the energy bands 1–3 keV, 3–6 keV and 6–10 keV, and the PDS lightcurves in the energy bands 13–25 keV and 25–80 keV, excluding the first 50 ks of our observation. The main pulse of the modulation presents a peak around phase 0.6. A secondary peak is visible around phase 0.1, which appears to be more prominent at higher energies. The pulse fraction, defined as $`(I_{max}I_{min})/I_{max}`$, with $`I_{max}`$ and $`I_{min}`$ the maximum and minimum count rate, are $`(43\pm 7)`$%, $`(57\pm 6)`$%, $`(62\pm 8)`$%, $`(7\pm 2)`$%, and $`(36\pm 3)`$% in the 1–3 keV, 3–6 keV, 6–10 keV, 13–25 keV, and 25–80 keV energy band, respectively. The pulse fraction in the 13–25 keV and 25–80 keV energy bands is strongly reduced by the presence of GX 9+1 in the PDS FOV (note however that the pulse fraction between 25–80 keV is somewhat higher because GX 9+1 is weaker at these energies). In Figure 6 we show the folded lightcurves in the energy bands 1–3 keV, 3–6 keV, 6–10 keV and 25–80 keV.

## 4 Discussion

We have discovered a serendipitous source, SAX J1802.7–2017, in a BeppoSAX/NFI observation of the bright atoll source GX 9+1. SAX J1802.7–2017 shows coherent pulsations at a period of 139.612 s, indicating that the compact object in this system is most likely an accreting magnetic NS. Pulse arrival times show a sinusoidal modulation at a period of 4.6 days, which we interpret as due to the orbital motion of the source. The values of the orbital period and of the pulse period, as well as the pulse fraction in the MECS energy range ($`40\%60\%`$), are consistent with the values usually found for HMXB systems (e.g. Oosterbroek et al. 1999; Bildsten et al. 1997).

The mass function of the system is

$$f(M)=\frac{M_C\mathrm{sin}^3i}{(1+q)^2}=\left(\frac{2\pi }{P_{orb}}\right)^2\frac{(a_x\mathrm{sin}i)^3}{G}17\pm 5M_{},$$

where $`M_C`$ is the companion star mass, $`G`$ is the gravitational constant and $`q`$ is the ratio of the compact object mass to the companion mass. The value of the mass function is not compatible with the system being a cataclysmic variable in which a white dwarf is orbiting a $`1M_{}`$ companion. For a NS mass of 1.4 M we have estimated a lower limit on the companion star mass, $`M_C11`$ $`M_{}`$ (90% confidence level), typical of a HMXB. The epoch of the NS superior conjunction (which would correspond to the mid-eclipse time) derived from the fit of the orbital parameters falls during the first 50 ks of our observation, when the count rate of SAX J1802.7–2017 is compatible with the background count rate. This is compatible with the possible presence of an eclipse during the first 50 ks. Assuming that the eclipse ends $`0.58`$ days after the beginning of our observation (52168.69 MJD, see Fig. 2) and taking into account that the zero epoch is $`0.11\pm 0.10`$ days, the half duration of the eclipse should be $`0.47\pm 0.10`$ days. The half-angle of the eclipse subtending the portion of circular orbit covered by the source is $`\theta =0.64\pm 0.14`$ rad. The eclipse duration is related to the inclination angle $`i`$ of the system by: $`R_C/a=(\mathrm{cos}^2i+\mathrm{sin}^2i\mathrm{sin}^2\theta )^{1/2}`$, where $`R_C`$ is the radius of a spherical primary star and $`a`$ is (approximately) the separation between the two stars (see Primini et al. 1976). Then a rough lower limit on the companion star radius is $`R_C14R_{}`$.

To estimate the luminosity of the source we performed a spectral analysis in the MECS range using the correct effective area for the position of the source in the MECS FOV (G. Cusumano, private communication). Although the low statistics of the data did not allow an accurate spectral analysis, we inferred a luminosity of $`5.6\times 10^{35}`$ erg/s over the 1.8–10 keV range by fitting the data to a power law (photon index $`0.1\pm 0.1`$ and $`N_H=1.7_{0.9}^{+0.8}\times 10^{22}`$ cm-2) and assuming a source distance of 10 kpc. The lightcurve of the system shows a large variability suggesting that the system could be accreting via a stellar wind. This is also suggested by its long pulse period and relatively short orbital period, typical of HMXBs accreting via stellar wind such as Vela X–1 and 4U 1538–52 (see e.g. Corbet 1986; Nagase 1989). Indeed SAX J1802.7–2017 shows similarities to 4U 1538–52, which has a long spin period of $`529`$ s, a short orbital period of $`3.7`$ days, and a mass of $`19.8\pm 3.3M_{}`$ for the companion star (Reynolds, Bell, & Hilditch 1992). However, while the averaged luminosity of 4U 1538–52 and Vela X–1 is around $`5\times 10^{36}`$ ergs s-1, the averaged luminosity of SAX J1802.7–2017 measured during our observation is at least one order of magnitude lower (assuming an upper limit to the distance of 10 kpc). The source seems therefore to be quite underluminous with respect to massive binaries accreting via stellar wind, although more observations are needed to confirm this conclusion.

Note that, as the source has not been reported by other X-ray astronomy satellites before (in particular we did not find any detected source at the position of SAX J1802.7-2017 in ROSAT images, although the field of GX 9+1 was not observed by ASCA) it might well be a transient HMXB. Most of the transient HMXBs are Be-star systems (Liu et al. 2000). However, if SAX J1802.7–2017 belongs to this class it should be atypical considering that its orbital period and pulse period do not follow the Corbet correlation (Corbet 1986). Note also that the presence of an eccentricity cannot be excluded by our data, and more observations are needed to confirm the orbital parameters of the source.

We thank Dr. Giancarlo Cusumano for providing us the corrected effective area for our source and the anonymous referee for useful suggestions. This work was partially supported by the Italian Space Agency (ASI), by the Ministero della Istruzione, della Università e della Ricerca (MIUR) and the Netherlands Organization for Scientific Research (NWO).