text:
# Exploring the Gamma Ray Horizon with the next generation of Gamma Ray Telescopes. Part 1: Theoretical predictions.

## 1 Introduction

Imaging Cherenkov Telescopes have proven to be the most successful tool developed so far to explore the cosmic gamma rays of energies above few hundred GeV. A pioneering generation of installations has been able to detect a handful of sources and start a whole program of very exciting physics studies. Now a second generation of more sophisticated Telescopes is starting to operate and is providing already with new exciting observations. One of the main characteristics of some of these new Telescopes , is the potential ability to reduce the gamma ray energy threshold below $`1020`$ GeV, helping to fill the existing observational energy gap between the detector on satellites and the ground-based installations.

In the framework of the Standard Model of particle interactions, high energy gamma rays traversing cosmological distances are expected to be absorbed by their interaction with the diffuse background radiation fields, or “Extragalactic Background Light” (EBL), producing $`e^+e^{}`$ pairs. The $`\gamma _{HE}\gamma _{EBL}e^+e^{}`$ cross section is strongly picked to $`E_{CM}1.8\times (2m_ec^2)`$ and therefore, there is a specific range in the EBL energy which is “probed” by each gamma ray energy .

This effect should lead to the existence of a “Gamma Ray Horizon”, limiting the feasibility of observing very high energy gamma rays coming from very far distances. The actual value of this horizon distance for gamma rays of a given energy, depends on the number density of the diffuse background radiation of the relevant energy range, which is traversed by the gamma rays. In the range of gamma ray energies which can be effectively studied by the next generation of Gamma Ray telescopes (from, say, 10 GeV to 50 TeV), the most relevant EBL component is the ultraviolet (UV) to infrared (IR) contribution.

Several models have been developed to try to predict that EBL density . These models do a quite complex convolution of the measurements of star formation rate, initial mass function and dust and light recycling history. The result is a set of relatively model-independent predictions which accuracy is improving as the quality of their astrophysics inputs improves with the new deep-field observations and which fits reasonably well the existing data.

Therefore, quantitative predictions of the Gamma Ray Horizon have already been made but, unfortunately, so far no clear confirmation can be drawn from the observations of the present generation of Gamma Ray Telescopes.

On the one hand, some very high energy gamma ray events might have been observed from Mkn 501, a blazar at redshift $`z0.03`$ . The mere observation of these events would somehow contradict the above predictions indicating, might be, the presence of a new mechanism violating the forementionned gamma-gamma reaction threshold, for which, for instance, Lorenz-Invariance violation has been advocated, as we’ll discuss later. Unfortunately, the statistics is scarce and for these events the actual systematic uncertainty in the energy determination might be large and hence the situation remains somewhat unclear. On the other hand, for the handful of presently well established extragalactic sources (all of them at modest redshifts), no clear observation of a common energy cutoff which could be attributed to the gamma absorption in the intergalactic medium instead of simply to internal source characteristics, has been established so far. Nevertheless, for Mkn 501 a clear exponential energy spectrum cutoff has been observed and, under the assumption that its origin is the EBL absorption, upper limits on the EBL density in agreement with the expectations have been placed .

The fact that the new generation of Cherenkov Telescopes will reach a considerably lower energy threshold than the previous one should be of paramount importance in improving the present experimental situation for, at least, two reasons:

* Lower energy points with a much smaller uncertainty, due to the steep spectra, will be added to the spectra of the already observed sources allowing to disentangle much better the overall flux and spectral index from the cutoff position in the spectrum fit.
* Sources at higher redshift should be observable, giving a stronger lever arm in constraining the predictions and the possibility of observing a plethora of new sources that will allow unfolding the emission spectra and the gamma absorption.

The goal of this work is to analyse the physics potential of this new generation of telescopes in the measurement of the Gamma Ray Horizon and more specifically its impact in the understanding of the various models and parameters involved in its predictions. For this, the work is structured in two parts.

The first part, which is the one covered in the present paper, concerns the theoretical predictions. In this part, first the definition of the terms used in this work and their calculational procedure are reviewed in detail. Then, the theoretical predictions obtained for the Gamma Ray Horizon for different EBL approaches and also for different cosmological parameter sets are presented. Finally plausible scenarios that could effect the Gamma Ray Horizon predictions are commented. The actual sensitivity to these models and parameters is discussed.

The second part, which will be presented in a fore-coming paper, will deal with the prospective on what the experimental scenario might look like and a discussion on how much one can expect to pinpoint the parameters of the theory (with special emphasis on the cosmological ones)in the extrapolated data scenario.

## 2 Description of the calculation

In this section the detailed calculation of the Gamma Ray Horizon in terms of the predicted EBL density spectra is presented. Our strategy has been performing the complete calculation without any approximation by using a numerical integration approach. Different ansatzs for the calculation of the EBL predictions have been also analysed to see their impact on the Gamma Ray Horizon prediction. Also, the dependence on all the intervening parameters has been kept explicit to be able to track the effect of the different hypotheses in the final prediction.

The optical depth can be written with its explicit redshift and energy dependence as

$$\tau (E,z)=_0^z𝑑z^{}\frac{dl}{dz^{}}_0^2𝑑x\frac{x}{2}_{\frac{2m^2c^4}{Ex(1+z^{})^2}}^{\mathrm{}}𝑑ϵn(ϵ,z^{})\sigma [2xEϵ(1+z^{})^2]$$
(1)
where $`x1\mathrm{cos}\theta `$, $`E`$ is the energy of the $`\gamma `$-ray,$`ϵ`$ is the energy of the EBL photon, z is the redshift of the considered source and $`n(ϵ,z^{})`$ is the spectral density at the given z’.

The predicted value of the optical depth depends on several physical parameters. A part from the dependence on the actual absorption process, which enters through the gamma-gamma cross section, and the direct dependence on the cosmological parameters $`H_0`$, $`\mathrm{\Omega }_M`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ introduced by the geodesic radial displacement function, the spectral energy density is also an input parameter.

There exists observational data with determinations and bounds of the background energy density at $`z=0`$ for several energies . The determinations come from direct measurements of the EBL density using instruments on satellites whereas the bounds, happen mostly in the infrared part of the EBL and come from extrapolations using galaxy counting. Given the difficulty of observing “cold galaxies” due to the zodiacal light background, they provide just lower limits.

Several models, which fit the observational data of $`n(ϵ,z=0)`$, have been suggested  (a set of predictions for the most significant models can be seen in figure 1). These models do not provide all the necessary information for our calculation: they provide a description of spectral density at $`z=0`$ while we need to know also the evolution of $`n(ϵ)`$ as a function of $`z`$.

In this note three different approaches, which represent somehow limiting cases in the complexity of the ansatz assumed, for the $`z`$ evolution of the EBL have been used. Their comparison should give a feeling on how much the predictions change with the complexity of the theoretical assumptions and, hence, they might provide a tool to estimate how large the theoretical uncertainties in the final predictions might be. Ordered in increasing complexity, these approaches are:

> 1. Burst of star formation at high redshift .
> 2. Parameterisation of the measured star formation rate .
> 3. Star formation rate and star evolution .

Computing the GRH for these extreme scenarios allows us to get an estimation of its maximum uncertainty due to the unprecise EBL knowledge.

### 2.1 Gamma Ray Horizon

For any given gamma ray energy, the Gamma Ray Horizon is defined as the source redshift for which the optical depth is $`\tau (E,z)=1`$. Therefore, the Gamma Ray Horizon gives, for each gamma ray energy, the redshift location $`z`$ of a source for which the intrinsic gamma flux suffers an e-fold decrease when observed on Earth $`z=0`$ due to the gamma-gamma absorption.

In practice, the cut-off due to the Optical Depth is completely folded with the spectral emission of the gamma source. But on the other hand, the suppression factor in the gamma flux due to the Optical Depth depends only (assuming a specific cosmology and spectral EBL density) on the gamma energy and the redshift of the source. Therefore, a common gamma energy spectrum behaviour of a set of different gamma sources at the same redshift is most likely due to the Optical Depth.

The goal of this note has been studying the effect of the cosmological parameters and the different spectral density models in the Gamma Ray Horizon predictions for the gamma ray energy region covered by the next generation of Gamma Ray Telescopes. The results of this study are presented in the next section.

## 3 Results

### 3.1 Optical Depth and Gamma Ray Horizon

As already mentioned, for any given energy of the gamma ray that travels through the universe, the probability of interaction with the EBL photons to create $`e^+e^{}`$ pairs has a strong dependence with the energy of the background photons. Roughly speaking, each gamma energy “probes” a different EBL photon energy and therefore, the trends of the EBL spectrum as a function of the photon energy $`ϵ`$ as well as its redshift evolution are reflected in the Optical Depth as a function of gamma energy $`E`$.

In figure 2 the Optical Depth for gamma rays coming from a set of different redshifts are shown as a function of the gamma ray energy. As already mentioned, in the comoving frame, the $`\gamma _{HE}\gamma _{EBL}e^+e^{}`$ reaction has the maximum probability when $`E_{CM}=Eϵ(1\mathrm{cos}\theta )1.8\times (2m_ec^2)`$. This means that the flat zone seen in figure 2 corresponds to gamma rays that interact mainly with EBL between roughly $`0.2`$ eV and $`1`$ eV (depending on the source redshift), where the density of EBL photons has a sharp break down (figure 1). On the other hand, while the gamma rays explore the peaks due to the star radiation and the absorption and reemission in the interstellar medium, the Optical Depth keeps increasing but with a non-constant slope.

In figure 3 (dotted line) the GRH that we get solving numerically the equation $`\tau =1`$ is shown. On the one hand, it is clear that from redshift $`z=1`$ onwards, it is quite flat, so that gammas of energy about $`<30`$ GeV could reach the Earth from any distance in the observable universe. On the other hand the GRH depends strongly on the redshift for $`z<1`$.

### 3.2 Spectral density

The Gamma Ray Horizon has been calculated for the three different evolutions of $`n(ϵ,z^{})`$ already mentioned. For the first and second approaches (“Burst of star formation at high redshift” and “Parameterisation of the star formation rate”) a model which defines the $`n(ϵ,0)`$ has to be chosen. In figure 3, a specific model for $`n(ϵ,0)`$  and specific values for the Madau curve ($`\alpha _M=3.8`$, $`\beta _M=1`$, $`z_b=1.5`$ and $`z_f=10`$), which agree with recent data , have been used.

The third model, namely the “Star formation rate and star evolution” assumption, is likely the closest to reality, so it is going to be used for all further studies in this work and, in fact, has already been used as a particular example in the previous sections.

The fact that we will stick to this approach and that no error bars are shown in figure 3, does not mean that this prediction is free from theoretical uncertainties. This model has a lot of inputs that come from cosmological measurements which have, in fact, quite large uncertainties . For instance, if one would assume a fit to the star formation rate in the redshift region for $`z>z_b`$ as the classical Madau curve  instead of a slowly decreasing rate, this would produce a sizeable change in the GRH ($`2040`$ GeV) prediction at large redshift. The uncertainty for low redshifts can be estimated by computing the GRH for several models of $`n(ϵ,0)`$, which produce a factor $`5`$ difference in the GRH energy prediction for $`z<<1`$ independently of the model used for its evolution.

### 3.3 Cosmological parameters

As we have already commented, some fundamental cosmological parameters such as the Hubble constant and the cosmological densities play also an important role in the calculation of the Gamma Ray Horizon since they provide the bulk of the $`z`$ dependence of our predictions.

Over the last few years, the confidence in the experimental determinations of these cosmological parameters has increased dramatically. To understand the effect of moving these parameters, the values quoted in table 1 have been used.

Before we discuss the impact of each one of these parameters in our predictions we would like to see how actually the observables that we will measure (Optical Depths and GRH) depend on the redshift $`z`$ to compare it with the redshift dependence of other observables. For that, we have plotted the prediction for their $`z`$ evolution in figure 4. In that figure, for each observable it is shown the prediction normalized to the value at $`z=0.01`$. For comparison, the $`z`$ variation of the Luminosity-Distance, used for the determination of the cosmological parameters using Supernova 1A observations and of the Geodesical-Distance, giving the gamma ray path length, are shown. One can see that the Optical Depth has a quite different behaviour depending on the gamma ray energy explored. The $`z`$ dependence is very pronounced at large redshifts for $`20`$ GeV gammas and approaches a “Geodesical-Distance”-like shape for $`2`$ TeV gamma rays. The reason for that is the actual shape of the EBL spectrum and its redshift evolution. To give a feeling of the actual average $`z`$ dependence of the Optical Depth, the prediction for a flat $`\nu I_\nu `$ EBL spectrum is also shown. Finally, the $`z`$ dependence of the inverse of the GRH energy is also shown.

Now, to understand up to which level the measurement of the GRH would allow to get information on $`H_0`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$ and $`\mathrm{\Omega }_M`$, the actual prediction of the GRH with the most sophisticated EBL approach has been repeated for a set of different values of these cosmological parameters.

For that, first each one of the parameters was changed $`\pm 3\sigma `$ from its best fit value, keeping the rest at their best fit value. The results are shown in figures 5 and 6. In figure 5, one can see that a $`3\sigma `$ variation leads to a change in the GRH prediction at high redshift which is of $`8`$ GeV for $`\mathrm{\Omega }_M`$ and $`4`$ GeV for $`\mathrm{\Omega }_\mathrm{\Lambda }`$, while keeping the GHR prediction unchanged at low redshifts as it was expected since for $`z<<1`$ the lookback time curve does not depend on $`\mathrm{\Omega }_M`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$. Figure 6 shows that a $`3\sigma `$ variation on $`H_0`$ also leads to $`5`$ GeV difference at high redshift but there is now also a sizeable difference for low redshifts, in contrast to the behaviour in the previous case.

The Hubble constant enters in the Optical Depth calculation as global factor and therefore its variation produces a global shift of the Optical Depth. Then, the flatest zone of the Optical Depth crosses the $`\tau =1`$ line at different redshift, which is seen in the GRH as a region where the logarithm of the GRH energy as a function of the logarithm of the redshift shows a hard slope (figure 6). Therefore the zone close to the hard slope region is very sensitive to $`H_0`$, since a 3 $`\sigma `$ variation changes $``$ $`50\%`$ the GRH.

The fact that the variations in the GRH due to the Hubble constant and due to the cosmological densities are qualitatively different leaves some room to disentangle both kind of parameters. Actually, in figure 7 it can be seen that a $`3\sigma `$ difference in each parameter produces a change of around $`10\%`$ in both cases at large redshift. But while decreasing redshift the effect due to changing the cosmological densities goes to zero, the effect due to the Hubble constant remains at around $`8\%`$. Therefore the precise determination of the GRH for $`z<0.1`$ and for $`z>0.1`$ will allow to perform independent measurement of both sets of parameters. In fact, above $`z>0.1`$ also the dependence on $`\mathrm{\Omega }_\mathrm{\Lambda }`$ and $`\mathrm{\Omega }_M`$ is different and therefore, precise measurements may provide a handle to measure both independently.

Finally, the sensitivity of the measurement of the GRH energy as a function of the redshift $`z`$ on each one of the parameters varied independently while keeping the rest at their best fit value, has been computed and is plotted in figure 8. In that figure the sensitivity for each parameter $`p`$ is actually defined as

$$S_p(z)p\frac{dE_{GRH}(z)}{dp}$$
(2)
in such a way that for a given uncertainty in the estimation of the GRH energy $`\mathrm{\Delta }E_{GRH}`$ the relative precision in the single-parameter determination of $`p`$ would be

$$\frac{\mathrm{\Delta }p}{p}=\frac{1}{S_p}\mathrm{\Delta }E_{GRH}$$
(3)
It is clear that the maximal relative sensitivity is for the $`H_0`$ parameter while for $`\mathrm{\Omega }_M`$ the relative sensitivity is, depending on the $`z`$ region, between one and two orders of magnitude smaller and around a factor 5 even smaller for $`\mathrm{\Omega }_\mathrm{\Lambda }`$. In this figure it is also clear that the sensitivities evolve differently with $`z`$ and therefore, if precision measurements are obtained it should be possible to fit simultaneously the three parameters. This possibility will be explored in detail in the second part of this work already mentioned at the introduction.

### 3.4 Beyond the “standard” calculation.

In the calculation presented above, the assumptions taken are based on our present knowledge of fundamental interactions, astrophysics and cosmology. Nevertheless, at such high energies and cosmological distances, for instance the effects from physics beyond the “Standard Model”, such as Quantum Gravity or Supersymmetry, could be important.

There are plausible scenarios “beyond” the present knowledge which could affect the GRH prediction and hence, should be considered. In the following we would like to comment on our understanding on how these effects, and other effects not considered in our calculations, could change the GRH predictions presented above.

#### 3.4.1 The absorption mechanism.

So far the only considered absorption mechanism has been the gamma-gamma interaction. As we have already seen, the gamma-gamma reaction has a strong dependence on the final state fermion mass and we have checked explicitly with our calculation that the contribution coming from Standard Model fermions other than the electrons adds a negligible absorption. As far as we know no extension of the Standard Model provides any alternative light final state particle not excluded already by the present accelerators that could add any significant amount of gamma-gamma absorption. Therefore, no sizable change in the GRH prediction can be expected in Standard Model extensions such a Supersymmetry due to modifications in the gamma-gamma cross section.

It is clear that the target for the high energy gammas could also be any other particle filling the intergalactic space. Therefore, it could be neutrinos, visible matter and barionic and non-barionic dark matter. Given the expected density for these targets and the present constraints in the dark matter candidates, we are not aware of any absorption mechanism with these targets that could add any sizable absorption contribution to the one of the gamma-gamma reaction for the gamma ray energy range considered in this paper and hence, give any sizable correction to the GRH prediction.

#### 3.4.2 Lorentz Invariance Violation.

High energy gamma rays traversing cosmological distances should notice the quantum fluctuations in the gravitational vacuum which unavoidably should happen in any quantum theory of gravitation. These fluctuations may occur on scale sizes as small as the Planck length $`L_P10^{33}`$ cm or time-scales of the order of $`t_P1/E_P`$ ($`E_P10^{19}`$ GeV).

These gammas will therefore experience a “vacuum polarization” correction which should be very small ($`O(E/E_{QG})`$ where $`E`$ is the gamma energy and $`E_{QG}`$ is an effective scale for Quantum Gravity, which might be as large as $`E_P`$ but might become measurable after the gamma has traversed cosmological distances. In this Quantum Gravity scenario emerges naturally the requirement of local ”violation” of the Lorenz-Invariance symmetry  providing as a direct effect an energy-dependent propagation speed for electromagnetic waves.

This local ”violation” of the Lorenz-Invariance symmetry changes the threshold condition for the $`\gamma \gamma f^+f^{}`$ reaction in a way that depends on the Quantum Gravity model considered and its effective scale . For plausible models, the correction to the GRH predictions turns out to be quite important and hence, deserves a detailed discussion, which we presented in reference .

#### 3.4.3 Astrophysical considerations.

The gamma-gamma cross section depends strongly on the gamma polarization state. The calculation made in this paper assumes unpolarized gammas but it might happen that the specific gamma ray source producing the high energy gammas under study produces them with a non-negligible degree of polarization. If that is the case, the cross section could change in such a way that the GRH could differ from the above predictions for that specific source.

Similarly, in the whole calculation it has been assumed that the distribution of the EBL was uniform and isotropic at any scale. Given the fact that we consider cosmological distances this assumption is quite plausible. Nevertheless, for any specific gamma ray source, it might happen that the “local” EBL density might differ sizably from the “average” one and therefore, the GRH observed from that source could be sizably different from our prediction.

These aspects and other of similar kind depending on the specific characteristics of the source and its environment should be easy to disentangle from the fundamental predictions if enough sources are observed at each redshift location range.

## 4 Conclusions

A complete calculation of the Gamma Ray Horizon (GRH) in the gamma ray energy range which will be covered by the next generation of Gamma Ray Telescopes has been presented and discussed in detail.

Several approaches for the calculation of the extragalactic background light (EBL) density ranging in complexity have been compared. That comparison shows that the uncertainties due to the EBL modeling might be quite large both at low and high redshift. Nevertheless, the results for the most realistic approaches agree in predicting that the GRH energy at large redshifts is of $`30`$ GeV and, hence, should be on the reach of the next generation Cherenkov Telescopes.

Following these predictions, the observable universe should become transparent to gamma rays of below $`30`$ GeV and then new, high redshift, high energy gamma ray sources should be observable by the next generation Cherenkov telescopes <sup>1</sup><sup>1</sup>1A different scenario but with a similar spirit was analysed in reference . There, the energy threshold for the new generation gamma-ray detectors was assumed to be $`>100`$ GeV, well above the asymptotic horizon given by the GRH predictions discussed here. Given this fact, to explore the cosmological potential that work does not use a direct measurement of the GRH but, instead, the observation of the halo radiation coming from secondary gamma emission in the electromagnetic cascade generated by the absolved primary VHE gammas.

If these new sources are abundant enough to make possible a precise measurement of the GRH energy as a function of the redshift, then either they can be used to place strong constraints in the EBL modeling or as a new technique allowing an independent determination of the cosmological parameters.

Exploring deeper this second scenario, the actual dependence of the GRH predictions on the cosmological parameters has been discussed in detail. This study shows the potential capability of a precise GRH energy determination as a function of the redshift $`z`$ to disentangle the relevant cosmological parameters and provide competitive determinations.

A more quantitative study on the actual experimental possibilities to fit the cosmological parameters with the foreseen observations is presented in the second part of this work.

## Acknowledgements

We want to thank T.Kneiske and K.Mannheim for many discussions, for providing us with their EBL spectra and for helping us to cross-check some preliminary results of this work. We are indebted to N.Magnussen for his cooperation in the early stages of this study. We want to thank our colleagues of the MAGIC collaboration for their comments and support. We want to thank also G.Goldhaber, P.Nugent and S.Perlmutter for their encouraging comments.


format: text

text:
# Longitudinal spin-fluctuations and superconductivity in ferromagnetic ZrZn2 from ab initio calculations

## Abstract

The recent discovery of superconductivity coexisting with weak itinerant ferromagnetism in the $`d`$-electron intermetallic compound ZrZn2 strongly suggests spin-fluctuation mediated superconductivity. Ab initio electronic structure calculations of the Fermi surface and generalized susceptibilities are performed to investigate the viability of longitudinal spin-fluctuation-induced spin-triplet superconductivity in the ferromagnetic state. The critical temperature is estimated to be of the order of 1 K. Additionally, it is shown that in spite of a strong electron-phonon coupling ($`\lambda _{\mathrm{ph}}=0.7`$), conventional $`s`$-wave superconductivity is inhibited by the presence of strong spin-fluctuations.

The generalization of the Bardeen-Cooper-Schrieffer (BCS) theory to electron-electron interactions by Kohn and Luttinger , paved the way for speculation about the possibility of non-$`s`$-wave, or “unconventional” superconductivity. Following the suggestion that a magnetically mediated interaction plays an important role in the superfluidity of liquid 3He , the search began for superconductivity in nearly magnetic metals where strong spin fluctuations might provide the pairing mechanism . Recent experiments on Sr2RuO4  have made it a strong candidate for exhibiting spin-triplet, possibly $`p`$-wave superconductivity. For a spin-singlet Cooper pair, where the electrons have anti-parallel spins, the presence of ferromagnetic spin-fluctuations will be antagonistic towards the development of such a superconducting state. However, the recent reports of the coexistence of ferromagnetism with superconductivity in UGe2  and ZrZn2  suggest a spin-triplet Cooper pairing, probably driven by such spin-fluctuations. Moreover, in ZrZn2, the disappearance of superconductivity at the same point as magnetism, and the sensitivity of its occurrence to sample purity  are perhaps the strongest indications yet that the superconductivity is intimately connected with the magnetism in this material.

Unlike other “magnetic” superconductors (e.g. borocarbides  , RuSr2GdCu2O8 ) where the magnetism and superconductivity occur in different parts of the unit cell, in both UGe2 and ZrZn2 it is the same itinerant electrons that are thought to form the Cooper pairs as well as produce ferromagnetism. Moreover, whereas some questions regarding the itineracy of 5$`f`$ electrons and the roles of the strong magnetocrystalline anisotropy and quasi-two-dimensional electronic structure can be raised with respect to UGe2, ZrZn2 is a three-dimensional intermetallic compound free of such effects. Discovered by Matthias and Bozorth  in the 1950s, it was initially of interest because of the presence of weak ferromagnetism, in spite of the fact that neither constituent was ferromagnetic. It has the C15 cubic Laves crystal structure, with a lattice constant of 7.393Å (13.97 a.u.) . The possibility of there being triplet pairing in high-purity C15 compounds like TiBe2 and ZrZn2 was first suggested by Fay and Appel . In this Letter we investigate the viability of their suggestion in the case of ZrZn2.

Firstly, we have calculated the electronic structure of ZrZn2 using the LMTO method . Exchange-correlation effects are described within the local spin density approximation (LSDA). Self-consistency was attained using 505 k-points within the irreducible wedge of the face-centered cubic Brillouin zone (BZ). The basis included $`s`$, $`p`$, $`d`$ and $`f`$ states for all atoms. Our results agree with the calculations of Jarlborg and Freeman  and Huang et al.  in the non-magnetic state and those of de Groot et al.  and Jarlborg et al.  in the spin-polarized one. The Fermi surface (FS), comprising four sheets, is shown in Fig. 1. Both non-magnetic (NM) and ferromagnetic (FM) calculations were performed at a series of different lattice parameters, and the magnetic moment was found to disappear near 13.47 a.u., where the calculated pressure is 45 kbar. This is just below the total energy minimum, indicating a calculated equilibrium lattice constant of about 13.6 a.u., in excellent agreement with the recent FLAPW calculation of Bruno et al. . This underestimation ($``$2.5%) is typical of the LDA, particularly when including $`f`$ states. Experiments confirm that the system is near the ferromagnetic instability, since the critical pressure for the disappearance of magnetism has been reported in the range 8.5—22 kbar .

Polarized neutron studies of the magnetization density  have shown that there is a significant spin-density along the Zr–Zr bond directions. Our calculation also shows this delocalization, but typically with 85% of the total moment on the Zr. The dominance of Zr is also reflected in the DOS at $`E_F`$, of which 70% is of Zr-$`d`$ character. Mattocks and Dixon  inferred an exchange splitting of 4.5 mRy (in a field of 8T) from their de Haas-van Alphen data for orbits on the $`\mathrm{\Gamma }`$-centered spheroid (band 30). This compares favorably with the value obtained from the current calculation at the equilibrium lattice constant (see Table I).

Secondly, we have studied the FS and its nesting properties. As seen in Fig. 1, it is rich in details. Four bands (27–30) cross the Fermi level. We used the relaxed lattice parameter given by LSDA ($`a=13.6`$ a.u.) corresponding to a total moment of 0.17 $`\mu _B`$/Zr. All these sheets show a strong Zr-$`d`$ character, though bands 27 and 30 exhibit a significant hybridization with Zn-$`p`$ (25% and 30% respectively). Furthermore, the DOS at $`E_F`$ is dominated by the contribution from band 29 (50%) and band 28 (32%), whereas band 30 contributes less than 1%. The effect of the spin-splitting is most noticeable in the change in the topology of the band 29 sheet, where the neck near the L point, present for majority spin electrons ($``$), disappears for minority ones ($``$). From Fig. 1, we expect strong intra and inter-band nesting features, especially along the $`<`$100$`>`$ direction. In order to understand how the nesting will actually affect the response of the electrons in this system, we have calculated the bare-band generalized static susceptibility,

$$\chi _0^{\sigma \sigma ^{}}(𝐪,0)=\underset{nn^{}𝐤}{}\frac{f_{\sigma n𝐤}(1f_{\sigma ^{}n^{}𝐤+𝐪})}{E_{\sigma ^{}n^{}𝐤+𝐪}E_{\sigma n𝐤}+i\delta },$$
(1)
where $`n`$ denotes the band index, $`\sigma `$ the spin, and $`f_{\sigma n𝐤}`$ the Fermi-Dirac functions. This sum was calculated on a mesh of more than $`5\times 10^5`$ k-points in the cube (shown in Fig. 1) using a tetrahedron interpolation technique similar to that of Rath and Freeman . These calculations were performed for both NM and FM cases along $`<`$100$`>`$, $`<`$110$`>`$ and $`<`$111$`>`$ in $`𝐪`$-space, at the LSDA equilibrium lattice constant. All these are peaked at $`q=0`$ and (rather surprisingly, given the FS topology) show very little structure at $`q>0`$. This was expected for the NM case because of the ferromagnetism of ZrZn2, but the absence of finite $`q`$ peaks in the FM case shows that this compound does not favor antiferromagnetic (AF) spinwaves as confirmed by our frozen spinwave calculations. In fact, the only significant $`q>0`$ peak is found in the $`<`$100$`>`$ direction at $`q=0.2`$ (in units of $`2\pi /a`$) in $`\chi _0^{}(q)`$. This peak, originating from intra-band contributions in band 29 near the corner of the cuboidal FS sheet (i.e. near the L point) is (obviously) not present in the NM case. The importance of this peak (and correspondingly the absence of the other expected ones) is due to the concentration of the DOS at $`E_F`$ near the border of the BZ (points X, K and L) where band 29 flattens considerably. In other words, this means that while nesting is present elsewhere, it is inhibited by the low DOS.

We now turn to the question of longitudinal spin-fluctuation-driven superconductivity in ZrZn2 as proposed by Fay and Appel . On the FM side of the transition, when the band structure is different for the two spins, we calculate the longitudinal coupling constant, $`\lambda _{\mathrm{sf}}^\mathrm{L}`$, from the generalized susceptibilities. In the FM region, within the RPA, the pairing potential can be written as

$$V^{\sigma \sigma }(𝐪)=\frac{I^2(𝐪)\chi _{0𝐪}^{\sigma ,\sigma }}{1I^2(𝐪)\chi _{0𝐪}^{\sigma ,\sigma }\chi _{0𝐪}^{\sigma \sigma }}.$$
(2)
In the approximation of a spherical FS (not unreasonable, given the dominant influence of band 29 and its FS topology), the longitudinal coupling parameter is given by ,

$$\lambda _{\mathrm{sf},l}^{\mathrm{L}\sigma }=N_\sigma (E_F)_0^{2k_{F\sigma }}𝑑q\frac{qV^{\sigma \sigma }(q)}{2k_{F\sigma }^2}P_l\left(1\frac{q^2}{2k_{F\sigma }^2}\right),$$
(3)
where $`P_l`$ is the Legendre polynomial. The objective is now to estimate the $`s`$ and $`p`$ components of $`\lambda _{\mathrm{sf}}`$ from our band structure calculation. The exchange integral, $`I`$, is obtained through the calculation of the Stoner enhancement, $`S`$, defined as the increase of the exchange splitting of the Zr potential divided by the energy of the applied magnetic field. The corresponding Stoner factor $`\overline{S}=11/S`$ is simply related to $`I`$ through $`I=\overline{S}/N`$, where $`N`$ is the DOS at $`E_F`$. We have calculated $`S(q)`$ for the FM case (see Table I) and for two AF spinwaves ($`q=0`$ and $`q=2\pi /a`$). We find that $`S`$ is quickly suppressed ($`S(q=2\pi /a)1.5`$) for AF spinwaves. This shows that ZrZn2 does not support AF fluctuations (which is consistent with the absence of peaks for $`q>0`$ in our calculated generalized susceptibilities) and that the $`q`$-dependence of $`I`$ cannot be neglected. From $`S(q)`$, we model $`I(q)=I_0/(1+b^2q^2)`$  with $`I_0=0.04`$ Ry and $`b^2=0.33`$ $`(a/2\pi )^2`$. Since the contribution to $`\lambda _{\mathrm{sf}}`$ from each spinwave mode is $`\frac{1}{2}S(q)\overline{S}(q)^2`$ , we get, as our first estimate, $`\lambda _{\mathrm{sf}}=1.2`$ by averaging over these three modes. We make the further approximation $`\chi _0^{\sigma \sigma }(q)=N_\sigma (E_F)`$ $`q`$, and calculate the longitudinal $`\lambda _{\mathrm{sf},l}^{\mathrm{L}\sigma }`$ in the $`s`$ ($`l=0`$) and $`p`$ ($`l=1`$) channels from Eqs. (2) and (3), taking a $`k_{F\sigma }0.6`$ $`2\pi /a`$, appropriate to the band 29 sheet. For the relaxed FM case in which $`\overline{S}=1.12`$, $`\lambda _{\mathrm{sf}}^\mathrm{L}`$ is negligible in the $`p`$ channel. However, moving closer to the FM transition, i.e. for $`\overline{S}=1.01`$, we get $`\lambda _{\mathrm{sf},0}^{\mathrm{L}\sigma }=1.9,2.0`$ (for $`\sigma =,`$) in the $`s`$ channel and $`\lambda _{\mathrm{sf},1}^{\mathrm{L}\sigma }=0.81,0.76`$ (for $`\sigma =,`$) in the $`p`$ channel. This shift is justifiable given that $`S`$ is extremely sensitive close to the transition. As noted by Fay and Appel, the $`s`$ component is much larger than the $`p`$ one and both diverge when $`\overline{S}1`$. These values are consistent with our previous estimate.

The electron-phonon interaction cannot be ignored and can even be expected to be rather large owing to the occurrence of conventional superconductivity in both Zr and Zn, and the large DOS at $`E_F`$. This suggests that the electron-phonon coupling could be sufficient to overcome the pair-breaking effects due to spin-fluctuations. The electron-phonon coupling constant, $`\lambda _{\mathrm{ph}}`$, can be expressed as $`\lambda _{\mathrm{ph}}=_i\frac{\eta _i}{M_i\omega _i^2}`$, where the sum runs over all atoms, $`i`$, with masses, $`M_i`$, and phonon frequencies, $`\omega _i`$, while the numerator, $`\eta _i=N_i(E_F)(V_i)^2`$, is the Hopfield parameter that describes the electronic contribution . Here, $`\eta `$ was calculated in the rigid muffin-tin approximation , i.e. retaining only dipolar terms and neglecting electronic screening of the ionic displacements. The values for $`\omega _i^2`$ are taken as one half of the Debye frequency of the atom $`i`$. Furthermore, we assume that the volume dependence of $`\omega _i^2`$ follows $`\sqrt{aB}`$, where $`B`$ is the (calculated) bulk modulus and $`a`$ the lattice parameter, which is reasonable as long as all small-$`q`$ phonons behave identically with pressure . As shown in Table I, $`\lambda _{\mathrm{ph}}`$ is of the order 0.7 near the calculated equilibrium volume and twice as large at the experimental lattice constant. Ignoring completely the destructive effects of spin fluctuations, and using $`\mu ^{}=0.13`$ in the McMillan formula , these correspond to respectable Tc’s of about 8 and 21 K, at the respective lattice constants. The decrease of $`\lambda _{\mathrm{ph}}`$ in the FM state can be attributed to the smaller total DOS at $`E_F`$ and its pressure dependence can be ascribed almost entirely to the behavior of the Debye frequency. In the vicinity of the ferromagnetic transition, this $`\lambda _{\mathrm{ph}}`$ is nevertheless insufficient to overcome the dominance of the spin-fluctuations as indicated by the large $`S`$ enhancements near the critical pressure (see Table I). A possibility exists for phonon-mediated superconductivity at larger pressures, i.e. well outside of the magnetic region when the Stoner factor would be further decreased to completely suppress spin fluctuations. However, our calculations at $`a=13.17`$ a.u. (equivalent to 160 kbar) indicate that $`\lambda _{\mathrm{ph}}`$ drops to 0.4 while $`S`$ is still 4.1 (giving a $`\lambda _{\mathrm{sf}}`$ of the order of 0.4), which combine to make the conditions non-favorable for phonon-mediated superconductivity. Note that the persistence of such a large Stoner enhancement over this range of pressures shows again the importance of the spin fluctuations and that the large values for $`S`$ in the FM region can be related to the observed absence of saturation of the magnetic moment  and the weak ferromagnetism. Since the magnetic moment, Stoner enhancement and spin fluctuations are associated with the Zr sublattice, phonon-mediated superconductivity might be envisaged to take place within the Zn sublattice, but such an explanation can be ruled out because the Zn contribution to $`\eta _i`$ is negligible.

Having presented evidence against the possibility of electron-phonon driven superconductivity, we now try to estimate $`T_c`$ from the longitudinal spin fluctuations. The typical spin-fluctuation cut-off frequency, $`\omega _{\mathrm{sf}}`$, can be estimated from the Stoner factor by $`\omega _{\mathrm{sf}}=1/(4NS)`$ , giving about 90 K at the relaxed lattice parameter. Using the Allen-Dynes formula, we arrive at a simplified expression for the superconducting transition temperature :

$$k_BT_c=\frac{\mathrm{}\omega _{\mathrm{sf}}}{1.2}\mathrm{exp}\left(\frac{1+\lambda _{\mathrm{ph}}+\lambda _{\mathrm{sf},0}^{\mathrm{L},\mathrm{T}}}{\lambda _{\mathrm{sf},1}^\mathrm{L}}\right).$$
(4)
Note that the rather strong electron-phonon interaction $`\lambda _{\mathrm{ph}}`$ contributes to the mass renormalization (numerator) and is detrimental to superconductivity in this case. Furthermore, the $`s`$-wave $`\lambda _{\mathrm{sf},0}^{\mathrm{L},\mathrm{T}}`$ contains both the longitudinal (L) and transverse (T) contributions. From the large measured electronic specific heat coefficient $`\gamma _{\mathrm{exp}}=47`$ $`\mathrm{mJ}\mathrm{mol}^1\mathrm{K}^2`$ , and our calculated values for $`\lambda _{\mathrm{sf},0}^{\mathrm{L}\sigma }`$ and $`\lambda _{\mathrm{ph}}`$, we infer a transverse contribution $`\lambda _{\mathrm{sf},0}^\mathrm{T}`$ of about 0.8. Using the values for the case $`\overline{S}=1.01`$, we get a $`T_c^\sigma =1.0,0.8`$ K for $`\sigma =,`$ respectively. These estimates are very approximate but they confirm, in our opinion, the viability of triplet $`p`$-wave superconductivity in ZrZn2.

In conclusion, we have shown that calculations based on our electronic structure results strongly support the idea that the recently observed superconductivity in ZrZn2  is indeed a result of triplet pairing, as suggested by Fay and Appel . This would lead to $`p`$-wave superconductivity, and since impurity scattering acts as a pair-breaker for pairing in the $`l0`$ channels, the high purity of samples is crucial. However, the experimental absence of superconductivity in the paramagnetic phase just above the critical pressure  is still unanswered by this theory which predicts an even larger $`T_c`$ in the NM region. The answer may lie in the peak of the transverse susceptibility, $`\chi _0^{}(q_{100})`$, at $`q=0.2`$ which could provide an attractive coupling that would naturally disappear outside of the FM phase. Finally, it might be worthwhile revisiting the properties of C15 compound TiBe2 under pressure since its electronic structure is very similar to that of ZrZn2 and conventional superconductivity would be more favored owing to the lighter masses of its constituents.

We are grateful to S. M. Hayden for communicating his results prior to publication and to B. L. Györffy and J. B. Staunton for helpful discussions. We also want to acknowledge the financial support of the Swiss National Science Foundation (GS) and the Royal Society (SBD).


text:
# Disordered Environments in Spatial Games

## I Introduction

Game theoretical methods have been applied, in a quite sucessfull way, to several different fields , biology being one of the most successful branches . In these evolutionary games, rewards are translated in terms of subsequent reproductive success, a more natural and clear scale than the rationality dependent ones used in, for example, economy. Of particular interest is the emergence and sustainability of cooperation, what has attracted a lot of attention as it poses a difficult problem from an evolutionary point of view since cooperators, by their nature, have their fitness decreased while interacting with defectors. Theories advanced to explain the evolution of cooperation usually consider either kin selection , reciprocal altruism  or group selection , the separation between them not always being clear. Group foraging and young raising , alarm calls , bacteria-infecting viruses , predator inspection by fish , birds female-female cooperation , cleaner fishes , are but a few examples of such acts in biological populations, ranging from simple to complex organisms. A number of other examples may be found in ref. .

The Prisoner’s Dillema game is generally studied as an archetypical model for reciprocal altruism. Each of the two players either cooperates ($`C`$) or defects ($`D`$), without knowledge of the opponent’s strategy. The result depends on the mutual choice and is given by the payoff matrix whose elements are: a reward $`R`$ (punishment $`P`$) if both cooperate (defect), $`S`$ (sucker’s payoff) and $`T`$ (temptation) if one cooperates and the other defects, respectively. Moreover, these quantities should satisfy the inequalities $`T>R>P>S`$ and $`2R>T+S`$. It is clear that, independent of the opponent’s choice, defecting is always the best bet. Thus, two basic evolutionary problems are $`i`$) the onset of cooperation, that is, under which conditions a given cooperative behavior can invade a population of defectors and, once established, $`ii`$) its stability, that is, under which conditions the population of cooperators is uninvadable. In other words, we are looking for Evolutionary Stable Strategies (ESS)  where there is stable mutual cooperation between individuals. In a random mating, infinite population of asexual (haploid) elements, where two pure strategies are present (cooperators $`C`$ and defectors $`D`$), it can be shown that defecting will be the most rewarding strategy (an ESS). If instead of only interacting once or a known finite number of times, the players have a probability $`w>0`$ of meeting again in the next round and remember the chosen strategies in previous encounters, more complex rules may be devised . For instance, in a round robin tournament, Axelrod  found that tit-for-tat (TFT) (starts cooperating and, after that, do what the opponent did in the previous step) was the most successful strategy, but other successful strategies have also been found . Nevertheless, TFT cannot invade a population of defectors before reaching a minimum population, what can be achieved in several ways .

Following a suggestion of Axelrod , Nowak and May  have shown how, cooperation can arise even with only pure, memoryless strategies, like always cooperate or always defect in the presence of spatial structure. They considered a deterministic cellular automaton consisting of a square lattice with near and next nearest neighbors interactions and self interaction. Since the agents are spatially localized, they are more likely to interact only with their nearest neighbors, differently from the standard, mean field-like approach, that consider an infinite, random mixing population. The reasons for this are manyfold: individuals usually occupy well-defined territorial regions, individuals do not move far from their places of birth (population viscosity ), interactions occur in places where animals usually meet such as water ponds, etc. After the combats, each player compares its total payoff with the ones of its neighbors and changes strategy, following the one with the greatest payoff among them. For a certain range of values of the payoff matrix, very complex spatial patterns show up with cooperators and defectors coexisting (spatial chaos). In these structured populations, cooperative strategies can build clusters in which the benefits of mutual cooperation can outweight losses against defectors, maintaining the population of cooperators stable. Irrespective of the initial state that may be chosen with either only one initial defector or a fraction of randomly distributed defectors, the asymptotic density depends only on the payoff matrix parameters. The actual values depend on the neighborhood chosen for the dynamics and whether self-interactions are included or not . This cooperation enhancement effect due to the spatial structure is also seen in other games. For instance, in the spatial version of the Hawk-Dove (HD) game , the density of doves (analogous to cooperative individuals) is increased by the spatial structure, although in this game there would exist polymorphism even in the absence of spatial effects. In a similar way, several strategies are more successful in a spatially structured population than otherwise .

In real populations, however, not all individuals interact the same number of times either due to the non sincronous character of the interaction or to the environment structure that prevents some of the contacts. Both the spatial and temporal aspects of the environment affect the interactions between individuals, being central issues in ecological and evolutionary theory. Thus, a natural question arises: how is cooperation affected by this inherent inhomogeneity? In other words, how robust is cooperation in the presence of disorder? In this paper we try a simple approach to this question by considering a regular lattice where some sites are empty. Dilution can be either quenched (fixed) or annealed (evolving). In the latter case, individuals may diffuse in the lattice, what will be considered in detail in a future work. Quenched vacancies (or defects) may account for the presence of environmental features (e.g., geographical) in the game, making some individuals have less neighbors than others. Deviations from the ordered lattice can also be achieved in several other ways. For instance, by allowing that some of the individuals also interact with distant ones, in a small-world network fashion . In , a diluted lattice was used with individuals interacting inside a region of radius $`r`$, and persistent polymorphism of $`C`$’s and $`D`$’s was found, unless $`r`$ was made too big, reaching the long range connectivity associated to mean field behavior, where the defector population dominates. Nevertheless, their treatment was quite brief and several interesting dynamical behaviors have passed unnoticed, as well as the important issue of whether disorder enhances or not the fraction of cooperators in the population.

Here we show that depending on the amount of disorder, cooperation can be enhanced, there being a point where a dynamical transition settles in, separating the region with spatial chaos from the one with localized groups of cooperators. The paper is organized as follows. In section II, the diluted version of the spatial prisoner’s game is presented and the order parameters are defined. In section III the main results are presented and, finally, in section IV, we present the conclusions and comments.

## II The Model

We consider the spatial version of the prisoner’s dilemma , placing the interacting elements in the vertices of a $`d`$-dimensional array, usually a hypercube, with periodic boundary conditions. The results presented here are for $`d=2`$. Differently from the original case, we allow that some of the sites may be empty. To describe the occupation of the $`i`$-site ($`1<i<N=L^d`$, where $`L`$ is the system linear dimension) we take $`n_i`$ to be either 1, if the site is occupied or 0, otherwise. In every generation, each individual assumes an unchangeable strategy from a set $`\mathrm{\Omega }`$ and, following ref. , we only consider the simplest case of pure strategies, $`C`$ and $`D`$, represented by the variable $`S_i=\pm 1`$, respectively. In each step (generation), the $`i`$-th individual ($`n_i0`$) combats with all other elements inside a given neighborhood $`𝒱_i`$, and accumulates a payoff $`p_i`$, depending on the chosen strategies, according to the reduced payoff table for the prisoner’s game : $`R=1`$, $`P=S=0`$ and $`T=b>1`$, reducing the problem to only one parameter (besides the density). For the case considered here, where the neighborhood $`𝒱_i`$ is restricted to the nearest neighbours of site $`i`$ and no self interactions are included, both $`C`$ and $`D`$ coexist in the region $`4/3<b<3/2`$ in which also the number of active sites (see below) is large. For $`1<b<4/3`$ cooperators dominate while for $`b>3/2`$, defectors are dominant. The player’s payoff is a measure of its reproductive success: when reproducing, the $`i`$-th element compares its own payoff with all $`j𝒱_i`$ and changes to the strategy of the site that has the greatest payoff in $`\{i\}𝒱_i`$. In this way, the global density is kept fixed since no empty site will be ever filled.

To characterize the macroscopic behavior of the system we introduce two order parameters. Let $`\rho _c(t)`$ represent the fraction of cooperators at a given time:

$$\rho _c(t)=\frac{1}{2N}\underset{i=1}{\overset{N}{}}(S_i+1)n_i$$
(1)
where $`N=L^d`$ is the total number of sites. Clearly, $`\rho _d(t)`$, the defectors density, is $`\rho _d(t)=\rho \rho _c(t)`$, where $`\rho `$ is the total density. Since we are interested in the long time regime, and the results depend on the choice of the frozen empty sites, we define the order parameter as the average over time ($`\mathrm{}`$) and over the realizations of the disorder ($`\overline{\mathrm{}}`$) of the asymptotic cooperators density, $`\rho _c=\overline{\rho _c(\mathrm{})}`$, for large $`N`$. Sometimes it is more useful to have the relative cooperators density, $`\rho _c/\rho `$. Thus, $`\rho _c=0`$ means that the population was fully invaded by defectors and $`\rho _c=\rho `$, by cooperators. An intermediate case, $`0<\rho _c<\rho `$ in which both strategies coexist is also possible. Moreover, it is interesting to know the fraction of active sites, that is, the fraction of elements that change strategy in time:

$$\rho _a(t)=\frac{1}{2N}\underset{i=1}{\overset{N}{}}(1S_i^tS_i^{t1})n_i$$
(2)
This defines our second order parameter, $`\rho _a=\overline{\rho _a(\mathrm{})}`$. Thus, $`\rho _a=0`$ means that all sites are frozen and $`\rho _a=\rho `$ that all elements are changing strategy.

## III Results

In fig. 1, the asymptotic cooperator density $`\rho _c/\rho `$ is plotted against the total occupation of the lattice, $`\rho `$. If the occupation fraction is near zero, almost all sites are isolated and do not change strategy since there is no combat at all, the asymptotic density being the same as the initial one: $`\rho _c=\rho _c(0)=\rho /2`$. Indeed, for small $`\rho `$, the curves for different values of $`b`$ merge. As the density increases, the probability of occurring pairs of occupied sites increases and, irrespective of the value of $`b`$, all $`CD`$ pairs will become $`DD`$, what can be seen as a decreasing curve of $`\rho _c`$ from the origin. Around $`\rho 0.3`$, as the number of interacting individuals increases due to the increasing clusters sizes, the dynamics will define the fate of each cluster and the curves for different values of $`b`$ depart: the higher $`b`$ is, the better is for defectors and the lower is the corresponding cooperator curve. Still further, when clusters of cooperators have a reasonable probability to be formed, their density start to increase recovering, as the total density approaches unity, the Novak-May results . The most interesting case happens if $`b`$ is in the active region, $`4/3<b<3/2`$, where the behavior is not monotonic and a sharp decrease in $`\rho _c`$ appears near $`\rho ^{}0.95`$. At the same point, the fraction of active sites presents a sharp increase, as can be seen in the inset of fig. 1. Notice that for the others regions, the behavior is almost uniform, with a very small number of active sites. Interestingly enough, depending on the region of $`b`$, the optimum density that maximizes the fraction of cooperators occurs at different values: for $`b<4/3`$, the more individuals, the better, and their maximum occurs at $`\rho =1`$; above $`b=3/2`$, on the other hand, the less occupied the network, the better, due to the exploitation by defectors; the maximum occurs in the limit $`\rho 0`$. In these cases, respectively, cooperators and defectors are in advantage when interacting. In the intermediate, active region, the behavior is non trivial and the maximum happens at the transition point $`\rho ^{}`$. Remarkably, this point is much above the site percolation transition, that for a square lattice is located at $`\rho _{perc}0.59`$ : although there is a connected, infinite cluster, regions of active sites are bounded (pinned) to small regions due to the presence of defects, as will be shown below. Below the transition, the approach to equilibrium is exponentially fast, $`\rho _c(\mathrm{})\rho _c(t)\mathrm{exp}(t/\tau )`$, $`\tau `$ diverging as one approaches the point $`\rho ^{}`$. This is because most of the interactions occur inside of the localized groups. Also, if the lattice is completely filled, the approach is fast. On the other hand, just above the critical value, the presence of defects, and the fact that only some of the groups are depinned makes the approach to equilibrium extremely slow, power-law like, what is reminiscent of disordered systems.

In the active region, for different initial concentrations of cooperators, it can be seen that almost up to $`\rho ^{}`$, the asymptotic density depends on the initial state, as shown in fig. 2. For small $`\rho `$, since the clusters are independent, obviously the higher the initial fraction of cooperators is, the higher $`\rho _c`$ will be. Notice that the optimum density for cooperation is a function of the initial density of cooperators: for initial concentrations below the height of the peak at $`\rho ^{}`$, $`\rho _c(0)0.54\rho `$, $`\rho ^{}`$ is the ideal total concentration for cooperators, while above, the smaller is $`\rho `$, the better. Slightly below $`\rho ^{}`$ all curves merge because the cluster boundaries are no longer completely pinned and the memory of the initial state is washed out.

Although both cooperators and active sites densities are useful parameters to describe the system asymptotic behavior, they are not sufficient in order to understand the complex dynamical behavior that arises in the presence of disorder. To clarify what is going on at the transition point, we also measured the persistence, the fraction $`P(t,t_w)`$ of sites that do not change strategy between an initial waiting time $`t_w`$ and the time $`t>t_w`$ , as can be seen in figs. 3 and 4. Differently from the fraction of active sites, the persistence is a very complex measure that depends on the whole time history of the system since $`t_w`$. For example, if the persistence does not go to zero, we know that there is a fraction of sites that flips only finitely many times  (blocking) and domain wall movements are constrained (pinning). That is precisely what happens for $`\rho <\rho ^{}`$, as shown in fig. 3: after an initial decrease, the persistence attains, for large times, a plateu whose value depends both on $`\rho `$ and $`t_w`$. Denoting this plateu by $`P(\mathrm{},t_w)`$, we notice that it goes to zero at $`\rho ^{}`$, as shown in the inset of fig. 3 for $`t_w=0`$, signalling a depinning transition. In the critical region, the behavior is power-law: $`P(\mathrm{},0)(\rho ^{}\rho )^{2.2}`$. It is important to notice that the contribution from isolated sites to the plateau is small, and most of the sites forming the plateau come from the infinite cluster. At $`\rho ^{}`$ and above, the interfaces are no longer constrained, the number of active sites suffer a sudden increase and all sites eventually change strategy. It has been shown for several models  that in the presence of disorder (but also in some homogeneous systems), the decay to the plateau, $`P(t,0)P(\mathrm{},0)`$, is exponential  at large times. The same behavior is found here for a non Hamiltonian model. In fig. 5 we present a snapshot of the system below $`\rho <\rho ^{}`$, showing both the empty and the pinned sites, as well as their strategies. This configuration is almost stable, many of the cooperator groups shown no longer change, and the small amount of active sites are either confined to a few regions or belong to one of the blinkers present. This has to be compared with the spatial chaos region above $`\rho ^{}`$, where there is no blocking and every spin flips infinitely many times. Since we are near the transition, the quantity of pinned sites is quite reduced.

For densities $`\rho >\rho ^{}`$, considering only elements that belong to the infinite cluster, the persistence goes to zero, as is shown in fig. 4. For $`\rho =1`$, the behavior is exponential, as shown in the inset of fig. 4. Notice, also in the inset (for $`\rho =0.99`$), that as we dilute the lattice, only the initial behavior is exponential. After this initial exponential decrease, whose length decreases as $`\rho `$ departs from 1, the behavior follows a power law, $`P(t,0)t^{\theta (\rho )}`$, and as we approach $`\rho ^{}`$ from above, the plateau starts developing. The exponent $`\theta (\rho )`$ is the persistence exponent and has non trivial values, depending on the total density. For example, $`\theta (0.98)3.6`$ and $`\theta (0.99)5.7`$.

Below the critical value, the groups of cooperators are localized and their borders cannot move because of the presence of the defects. We say the groups are pinned by the environment. This explains why the fraction of cooperators is larger in this regime. As the number of defects decreases, the groups start to become depinned, the interfaces start moving and interference effects between the groups appear, explaining the sudden increase in the number of active sites. From the persistence data, we can see that below the critical value, there is a constant fraction of sites in the infinite cluster that never change strategy, meaning that the clusters are localized. As we approach the deppining transition, some of the groups are depinned and the persistence decreases. Since some of the groups are still pinned, we observe a constant plateau that decreases as $`\rho `$ approaches $`\rho ^{}`$. Notice that below the transition, the number of active sites is not zero, but small, due to blinking sites.

## IV Conclusions

One of the main features of the Prisoner’s dilemma, responsible for its widespread use in the problem of cooperation, is its robustness. Here we have shown that the system is able to sustain cooperation even under the presence of lattice disorder, what was already suggested in ref. . However, what had been unnoticed was the dynamical transition as a function of the amount of quenched disorder as well as the fact that the disorder may enhance cooperation. Dilution changes the scenario presented by Nowak and May  for the filled lattice in a dramatic way. The spatial chaos is absent when the disorder is above a given amount, what is reflected by the small number of active sites and the pinning of the interfaces. It must be noticed that even without disorder, we are not dealing with a system presenting coarsening as the movement of the interfaces is not ruled by surface tension as, e.g. in the Ising model.

Besides presenting a dynamical transition in the presence of disorder, dilution can even enhance the fraction of cooperators in the population. Indeed, in the region $`4/3<b<3/2`$, the relative density of cooperators in the population can reach a value almost $`40\%`$ higher than in the undiluted case, the maximum happening at the transition point.

Some questions, however, are still open. For example, if the condition of strong population viscosity is relaxed, disorder is no longer quenched and diffusion is allowed, what happens with the cooperator groups? Obviously, the results depend on the chosen rule for the diffusing elements, and there are several biologically motivated rules. Work is in progress in this direction. Moreover, from the point of view of the study of persistence, it would be interesting to study the behavior of the $`P(t,t_w)`$ above $`\rho ^{}`$, for $`1t_wt`$ (instead of $`t_w=0`$) as well as to compare the results presented here with the site diluted Ising model.

Acknowledgments: We are grateful to H. Chaté, N. Lemke and L. Peliti for discussions and suggestions and Y. Levin for a critical reading of the manuscript. Work partially supported by CNPq and PROPESQ-UFRGS.


text:
# Deep Radio, Optical, and Infrared Observations of SGR 1900+14

## 1 Introduction

Soft $`\gamma `$-ray Repeaters (SGRs; see Hurley 2000 for a recent observational review) are sources both of bursts of hard X-rays/$`\gamma `$-rays and softer quiescent X-ray emission. They are generally thought to be young ($`<10^4`$ yr) neutron stars with extremely strong magnetic fields, or magnetars (e.g. Duncan & Thompson 1992; Thompson & Duncan 1993). This belief was motivated by the energetics of their bursts (Paczyński 1992; Thompson & Duncan 1995), their associations with young supernova remnants (SNRs; Evans et al. 1980; Kulkarni & Frail 1993; Vasisht et al. 1994), and the detection of X-ray pulsations with long (5–10 s) periods and large ($`10^{11}\text{ s s}^1`$) spin-down rates (Kouveliotou et al. 1998; Hurley et al. 1999d). However, more recent observational evidence (Lorimer & Xilouris 2000; Hurley et al. 1999b; Kaplan et al. 2002) has shown that in most cases the SGR/SNR associations are false and can be attributed to random chance (Gaensler et al. 2001). Instead of being constrained by the SNR associations, the small ages of the SGRs are now limited by the proximity of the SGRs to sites of very massive star formation (Fuchs et al. 1999; Vrba et al. 2000).

An alternate theory, first proposed (van Paradijs, Taam, & van den Heuvel 1995; Chatterjee, Hernquist, & Narayan 2000) for the emission from the related Anomalous X-ray Pulsars (AXPs; see Mereghetti 2000), and now for the SGRs as well (Marsden et al. 2000, 2001; but see Duncan 2001), is accretion from a “fallback” disk of material produced during the supernova explosion onto a central neutron star. In this model, the neutron stars are spinning near the equilibrium period which is proportional to the dipole field strength. This model then explains the spin-down and period clustering of the AXPs: for a typical field strength of $`10^{12}`$ G, the resultant spin period is $`8`$ s.

To date, no confirmed persistent counterparts have been found for SGRs at optical/IR or radio wavelengths (e.g. Kaplan et al. 2001), although Eikenberry et al. (2001) may have found an infrared counterpart to SGR 1806$``$20. It generally has been the case that when counterparts are identified at other wavelengths, especially optical or IR, there is considerable progress in our understanding of these objects. For example, the identification of an optical/IR counterpart of the AXP 4U 0142+61 (Hulleman, van Kerkwijk, & Kulkarni 2000) has posed the strongest challenge to accretion models. In view of this, it is important to press upon sensitive multi-wavelength observations of SGRs. Here we present optical, IR, and radio observations of the very accurately localized SGR 1900$`+`$14 (Frail, Kulkarni, & Bloom 1999; Hurley et al. 1999d).

### 1.1 SGR 1900$`+`$14

SGR 1900$`+`$14 had been a relatively innocuous source until late 1998. The source of sporadic bursting emission (Mazets, Golenetskij, & Guryan 1979; Kouveliotou et al. 1993), it had been associated with a soft, persistent X-ray source (Hurley et al. 1996) and a SNR, G042.8$`+`$00.6 (Vasisht et al. 1994). On 1998 August 27, the spacecraft Konus detected an intense burst of $`\gamma `$-rays (Hurley et al. 1999a). Arrival-time localization (Hurley et al. 1999c) soon identified SGR 1900$`+`$14 as the source. Prompt radio observations detected a fading, non-thermal source (Frail et al. 1999) coincident with the X-ray source (Hurley et al. 1999d), and timing of the X-ray emission revealed a 5.16-s pulsar (Hurley et al. 1999d) with magnetar-like spin-down (Kouveliotou et al. 1999). The radio source associated with the 1998 August 27 burst is located at $`\alpha _{2000}=19^\mathrm{h}07^\mathrm{m}14\stackrel{\mathrm{s}}{\mathrm{.}}33`$, $`\delta _{2000}=+09\mathrm{°}19\mathrm{}20\stackrel{}{\mathrm{.}}1`$, with uncertainty $`\pm 0\stackrel{}{\mathrm{.}}15`$ in each coordinate (Frail et al. 1999).

Searches for SGR 1900$`+`$14 at wavelengths other than X-ray and $`\gamma `$-ray have not detected any quiescent counterpart (Hurley et al. 1996; Lorimer & Xilouris 2000; Eikenberry & Dror 2000; Vrba et al. 2000). However, Lorimer & Xilouris (2000) did detect a young ($`10^{45}`$ yr old) radio pulsar (PSR J1907$`+`$0918) located in the vicinity111PSR J1907$`+`$0918 is located at $`\alpha =19^\mathrm{h}07^\mathrm{m}22\stackrel{\mathrm{s}}{\mathrm{.}}441`$, $`\delta =+09\mathrm{°}18\mathrm{}30\stackrel{}{\mathrm{.}}76`$, with uncertainty $`\pm 0\stackrel{}{\mathrm{.}}05`$ in each coordinate (Lorimer & Xilouris 2000). of both SGR 1900$`+`$14 and G042.8$`+`$00.6, suggesting that G042.8$`+`$00.6 could be associated with PSR J1907$`+`$0918 just as easily as with SGR 1900$`+`$14.

Vrba et al. (2000) have suggested that SGR 1900$`+`$14 is somehow related to a compact cluster found in its vicinity (about $`12\mathrm{}`$ away). A similar cluster has been identified near SGR 1806$``$20 (Fuchs et al. 1999). We note, though, that no such cluster is found in the vicinity of SGR 0526-66 in the HST images presented by Kaplan et al. (2001). Furthermore, SGR 1806$``$20 and SGR 1900$`+`$14 are found at low Galactic latitudes in the inner Galaxy and the chance coincidence probability with star clusters, the rare Luminous Blue Variables (LBVs) and SNRs is not negligible (e.g. Gaensler et al. 2001). Indeed, there is now significant doubt about the association of SGR 1806$``$20 with the SNR G10.0$``$0.3 (Hurley et al. 1999b; Kaplan et al. 2002). Nonetheless, the presence of star clusters in the vicinity of SGR 1900$`+`$14 appears intriguing, but while our radio maps are deeper than previous images (e.g. Vasisht et al. 1994) we find no radio source or nebulae associated with this star cluster (see Figures 35).

For the X-ray properties of SGR 1900$`+`$14 in what follows, we use the power-law $`+`$ blackbody fit to the second epoch (the one closer to normal quiescence) of the recent Chandra observations by Fox et al. (2001): $`N_H=2.3\pm 1\times 10^{22}\text{ cm}^2`$, power-law index $`\mathrm{\Gamma }=1.9\pm 0.3`$, blackbody temperature $`kT0.49\pm 0.03`$ keV, and an un-absorbed 0.5–10 keV X-ray flux $`f_\mathrm{X}^\mathrm{U}=2.0\times 10^{11}\text{ erg s}^1\text{ cm}^2`$ (also see Kouveliotou et al. 2001). This is consistent with the ASCA measurements of Hurley et al. (1999d). We adopt a distance $`d=5`$ kpc (Hurley et al. 1999d), which is based on hydrogen absorption. The hydrogen column density implies a visual extinction $`A_V=12.8\pm 0.8`$ mag, assuming $`N_H=1.79\pm 0.03\times 10^{21}A_V\text{ cm}^2`$ (Predehl & Schmitt 1995). These values differ significantly from those assumed for the star cluster ($`d=12`$–15 kpc, $`A_V19.2`$ mag; Vrba et al. 1996), and this difference may make the case for associating SGR 1900$`+`$14 and the cluster less compelling.

## 2 Observations

We have carried out a campaign at optical/IR and radio wavelengths to detect SGR 1900$`+`$14. The radio data can also determine the existence of a low surface-brightness nebula surrounding SGR 1900$`+`$14, and localize SGR 1900$`+`$14 relative to nearby SNRs. See Table 1 for a log of all the observations.

### 2.1 Optical

The optical data consist of 5 exposures taken by the Space Telescope Imaging Spectrograph (STIS) aboard the Hubble Space Telescope (HST), totaling 5195 s. The data were taken in imaging mode using the clear (50CCD) filter, giving a bandpass centered at roughly $`585.2`$ nm and a full-width at half-maximum (FWHM) of $`400`$ nm. The first exposure had a duration of 30 s, and was used for astrometry (see below). The remaining exposures were $`1300`$-s long and were dithered in the usual manner to improve resolution and mitigate the effects of CCD irregularities.

The reduction followed the standard procedure for STIS images. We followed the standard drizzling reduction procedure (Fruchter & Mutchler 1998) expanding the individual pixels by a factor of two. This drizzling changed the pixel scale from $`0\stackrel{}{\mathrm{.}}05\text{ pixel}^1`$ to $`0\stackrel{}{\mathrm{.}}025\text{ pixel}^1`$. Typical stellar sources had a FWHM of $`0\stackrel{}{\mathrm{.}}095`$; see Figure 1.

For the astrometry, we selected from the USNO-A2.0 catalogue (Monet 1998) all 76 stars that overlapped with a 30-s R-band image taken on 28 August 1998 with the Low-Resolution Imaging Spectrograph (LRIS; Oke et al. 1995) on the Keck telescope. For the 54 stars that were not overexposed, we measured centroids and corrected for instrumental distortion using a bi-cubic function determined by J. Cohen (1997, private communication). We used these to solve for the zero-points in each coordinate, the plate scale, and the position angle on the sky. After rejecting 8 outliers (residual larger than $`0\stackrel{}{\mathrm{.}}8`$), the root-mean-square residuals were $`0\stackrel{}{\mathrm{.}}29`$ in each coordinate. This is somewhat larger than we found for other projects, but not atypical for relatively crowded Galactic fields on the Northern hemisphere. Given this, the astrometry is tied to the USNO-A2.0 system at the $`0\stackrel{}{\mathrm{.}}05`$ level. The typical accuracy with which the USNO-A2.0 catalogue is on the International Celestial Reference Frame (ICRF) is about $`0\stackrel{}{\mathrm{.}}2`$ (Deutsch 1999).

The astrometry was transferred to the STIS image using 43 stars close to the radio position, solving again for plate scale, rotation and zero-points. The rms residuals were $`0\stackrel{}{\mathrm{.}}02`$ in each coordinate, giving an overall astrometric uncertainty of $`0\stackrel{}{\mathrm{.}}21`$ in each coordinate. The 3-$`\sigma `$ limiting magnitude, based on calibration using the photometric keywords in the image header, is $`29.0`$, or $`f_{\nu ,\mathrm{min}}=10`$ nJy 222Note that STIS magnitudes are in the STMAG system, where $`m=21.12.5\mathrm{log}F_\lambda `$, with $`F_\lambda `$ in units of $`\text{ erg s}^1\text{ cm}^2\text{ Å}^1`$. All other magnitudes are Vega-based..

### 2.2 Infrared

The infrared data were taken with the Near Infrared Camera (NIRC; Matthews & Soifer 1994) mounted on the 10-m Keck I telescope. The observations had photometric conditions and very good seeing ($`0\stackrel{}{\mathrm{.}}35`$ at $`K_s`$-band). The data were reduced with custom IRAF software. The images were dark subtracted, flat-fielded, and corrected for bad pixels and cosmic rays. We then made object masks, which were used to construct improved flat fields for a second round of reduction. The data were finally registered, shifted, and co-added; see Figure 2.

The astrometry was performed with 50 stars from the deep STIS image. This solution had errors of $`\pm 0\stackrel{}{\mathrm{.}}04`$ (in each coordinate) relative to the STIS image, or again $`\pm 0\stackrel{}{\mathrm{.}}21`$ relative to the ICRF. For photometric calibration we used the standard star SJ 9101 (Persson et al. 1998). The 3-$`\sigma `$ limiting magnitudes are $`K_s20.8`$ and $`J22.8`$, corresponding to $`f_{\nu ,K_s}3.1\text{ }\mu \text{Jy}`$ and $`f_{\nu ,J}1.2\text{ }\mu \text{Jy}`$, based on the calibration of Bessell, Castelli, & Plez (1998).

### 2.3 Radio

The 1.4-GHz data were taken at the Very Large Array (VLA) in its C and D configurations using the standard continuum mode, with $`2\times 50`$ MHz bandwidth. The separate observations were calibrated independently, following the standard procedure in AIPS. The phase calibrator was JVAS J1925+2106, and the flux calibrator was 3C 48. The data were then self-calibrated and combined for imaging. The nominal beam size was $`20\mathrm{}`$, and the rms noise is $`0.5\text{ mJy beam}^1`$; see Figure 3.

To aid in the detection of diffuse structure we also imaged using a $`uv`$-taper of $`2.5\text{ k}\lambda `$, giving an effective beam size of $`60\mathrm{}`$. The rms noise is $`2\text{ mJy beam}^1`$, largely due to confusion from unresolved sources in the Galactic plane; see Figure 4.

The 332 MHz data were taken at the VLA in spectral-line mode with 32 channels of 97 kHz each to aid in the excision of radio-frequency interference (RFI). We reduced the observations in AIPS, following the procedure for low-frequency data reduction333http://rsd-www.nrl.navy.mil/7213/lazio/tutorial/index.html (see also LaRosa et al. 2000). First we removed any egregious RFI. Then we phase (using 3C 380) and gain (using 3C 286 and 3C 48) calibrated the fringe visibilities. We again searched for and removed RFI, then averaged the channels together. We mapped and self-calibrated the data iteratively, using 3-dimensional Fourier transforms in IMAGR, the basic AIPS deconvolution and imaging task. While imaging we used 25 planar facets arranged in a $`5\times 5`$ grid, each of which had $`256\times 256`$ $`13\mathrm{}`$ cells. These facets were used to mitigate the effects of non-coplanar arrays that are present at low frequencies and with wide fields (Cornwell & Perley 1992; Perley 1999). Finally, we re-assembled the facets and flattened them into a single planar map using the task FLATN. The rms noise is $`2\text{ mJy beam}^1`$, limited primarily by confusion due to unresolved emission in the Galactic plane and to sidelobes from the bright W49 complex. The final beam size is $`50\mathrm{}`$, comparable to that of the tapered 1.4-GHz images; see Figures 5 and 6.

## 3 Analysis & Discussion

### 3.1 Optical/IR

For the infrared and optical observations the astrometric uncertainty is $`0\stackrel{}{\mathrm{.}}21`$. Combined with the $`0\stackrel{}{\mathrm{.}}15`$ uncertainty from Frail et al. (1999), this gives us a 99% confidence radius of $`0\stackrel{}{\mathrm{.}}79`$. We examined all sources in or near a circle of this radius around the position of Frail et al. (1999). In the deepest observation, the $`K_s`$-band image, we detect 7 sources that we label A–G on Figures 1 and 2 (only sources A–D are visible in Figure 1). Of these, C and D are significantly blended on all the images, as are sources C and G. Sources E, F, and G are too faint to measure on the STIS image, and sources E and F are too faint to measure on the $`J`$-band image. Because of the degree of blending, we used a point-spread function (psf) fitting and subtraction technique to perform photometry on the infrared data, using the DAOPHOT package. First we measured the psf’s of a number ($`30`$) of stars in the field without close neighbors, and constructed an average psf for each image. We then fit this psf to each source in Figures 1 and 2, subtracted the fitted stellar profile from the images, and used the subtracted images to fit the sources that were hidden by the original source. After iterating this procedure, we were able to obtain consistent magnitudes for most of the sources. For the sources that we were able to fit, all had sizes consistent with unresolved sources. For the STIS data, we adjusted the measured magnitudes to a nominal infinite aperture with a $`0.13\pm 0.03`$ mag correction, based on the observed STIS psf444http://www.stsci.edu/instruments/stis/documents/ihb/stis\_cy11\_ihb.pdf. We present the positions and magnitudes of sources A–G in Table 2.

In Figure 7, we show a color-color diagram for the sources in Table 2 and 37 field stars for which we were able to perform reasonably accurate photometry in all bands. The sources from Table 2 all lie in the region populated by the other sources. In Figure 7 we also show model main-sequence colors (based on Cox 2000, p. 388) for $`A_V=5`$–10 mag. These fit the majority of the points, with the exception of source D. However, since source D was blended with source C (making the IR photometry less reliable) and it is still in the region of the color-color diagram populated by other stars (some of whom were slightly blended themselves), we consider it unlikely that D has significantly non-stellar colors, and therefore is probably not a counterpart to SGR 1900$`+`$14 (both the disk models for AXPs \[Perna, Hernquist, & Narayan 2000\] and the optical/IR detection of an AXP \[Hulleman, van Kerkwijk, & Kulkarni 2000\] have non-stellar colors, and we expect something similar for SGRs). In Figure 8 we show color-magnitude diagrams for the sources in Table 2 and 37 field stars. Again, the sources populate the same region (although the field stars were chosen to be somewhat brighter than the faintest of the candidate sources), and are well bounded by model main sequences with $`d=8`$ kpc and $`A_V=5`$–10 mag. Source D, in particular, has slightly strange colors (also as seen in Figure 7) but its magnitudes are well within the range of the main sequence. Main sequences with $`d=12`$–15 kpc and $`A_V19.2`$ mag, such as those plotted in Vrba et al. (2000), are too faint and too red to match these sources. Therefore the sources in Table 2 are likely to be in the foreground of the cluster found by Vrba et al. (2000).

Source G in Figures 1 and 2 and in Table 2 is too faint for accurate photometry given its proximity to other sources. It could be a counterpart of SGR 1900$`+`$14, and we require additional deep, high-resolution infrared imaging to settle this issue. Regardless, it is near our detection limit, so we can place general upper limits to the emission from SGR 1900$`+`$14 of $`m_{50\mathrm{C}\mathrm{C}\mathrm{D}}29.0`$ mag, $`J22.8`$ mag, and $`K_s20.8`$ mag, consistent with previously published limits (Eikenberry & Dror 2000; Vrba et al. 2000). We show these limits along with the CXO X-ray spectrum (Fox et al. 2001) in Figure 9. These limits translate into unabsorbed (“U”) X-ray-to-infrared flux ratios $`f_X^\mathrm{U}/(\nu _Jf_{\nu ,J}^\mathrm{U})160`$ and $`f_X^\mathrm{U}/(\nu _{K_s}f_{\nu ,K_s}^\mathrm{U})1.1\times 10^3`$ (assuming $`A_V=12.8`$ mag). This last limit is only a factor of $`10`$ less than the detected X-ray-to-$`K_s`$ limit found for the AXP 4U 0142$`+`$61 (F. Hulleman 2001, personal communication; see Figure 9) and thus we can exclude the basic disk model for AXP/SGR emission (Perna et al. 2000; see also Kaplan et al. 2001). If a disk does exist, then it must be truncated at $`R10^{10}`$ cm. As seen for SGR 0526$``$66 (Kaplan et al. 2001), the optical/IR limits require that the X-ray power-law break between the X-ray and the optical bands, but allow for a continuation of the X-ray blackbody.

We were able to detect but not photometer many luminous members of the massive star cluster that Vrba et al. (2000) proposed as the origin of SGR 1900$`+`$14, as we only had a small number of exposures where the central M giant stars were unsaturated. In the remaining exposures these stars were too saturated to allow good psf subtraction, and their scattered light, bleed trails, and readout trails prevent accurate measurement of the stars in the cluster.

### 3.2 Radio

We do not detect either SGR 1900$`+`$14 or PSR J1907$`+`$0918 in the radio images. We therefore have 3-$`\sigma `$ limits on point source emission of 6 mJy (332 MHz) and 1.5 mJy (1.4 GHz). We also have 3-$`\sigma `$ limits on extended emission of $`6.1\text{ mJy arcmin}^2`$ (332 MHz) and $`6.2\text{ mJy arcmin}^2`$ (1.4 GHz).

We were able to detect SNRs G042.8$`+`$00.6, G043.3$``$00.2, G041.1$``$00.3, and possibly G043.9$``$01.6 at 332 MHz (Figure 5). However, G041.1$``$00.3, G043.3$``$00.2, and G043.9$``$01.6 are too attenuated by the primary beam of the VLA antennas for the flux measurements to be useful. For G042.8$`+`$00.6, we measure $`S_{332\text{ MHz}}=1.24\pm 0.08`$ Jy, excluding the contributions of point sources (see Figure 6), but we are likely missing some flux on the shortest baselines so this flux may be underestimated by 5%–10%.

#### 3.2.1 Serendipitous Sources

The 332-MHz image (Figure 5) shows three sources besides the known SNRs that have shell-like morphologies but are in neither the current version of the Galactic SNR catalog (Green 2000) nor any online database. We consider these sources to be candidate SNRs, but with only limited data at other wavelengths we are unable to confirm an SNR classification for any of the three. The sources could be H II regions or other similar sources and require additional confirmation. Here we present what limited information we have on each.

The first source is G043.5$`+`$00.6 (see Figure 10). Its center is at position $`\alpha _{2000}=19^\mathrm{h}08^\mathrm{m}26^\mathrm{s}`$, $`\delta _{2000}=+09\mathrm{°}42.0\mathrm{}`$. It is circular, with radius $`8.7\mathrm{}`$. The total flux at 332 MHz is $`1.05\pm 0.08`$ Jy (excluding the point source to the south-east), with an average surface brightness of $`3.9\text{ mJy arcmin}^2`$. The edge appears to be brighter, with typical surface brightness of 6–8$`\text{ mJy arcmin}^2`$, while the interior is more typically 1–3$`\text{ mJy arcmin}^2`$.

Without a measurement at another frequency, it is hard to determine the properties of G043.5$`+`$00.6. We were not able to find any archival X-ray or radio data at this position that could definitively classify G043.5$`+`$00.6. However, there is a marginal detection in the 1.4-GHz NRAO VLA Sky Survey (NVSS; Condon et al. 1998). The source is present in the NVSS with the same general morphology although it is too faint to accurately measure a flux. This faintness, though, suggests that G043.5$`+`$00.6 is not a thermal source, as it would have a NVSS flux of $`1`$ Jy and would then be easily detectable. There are two X-ray sources from the EXOSAT medium-energy slew survey catalog (Reynolds et al. 1999) that lie in or near G043.5$`+`$00.6: EXMS B1906+095 and EXMS B1906+094. Both sources had been associated with the X-ray binary 4U 1907+097, which is within their $`10\mathrm{}`$ position uncertainties, but the positions of the EXOSAT sources are closer to G043.5$`+`$00.6 than to 4U 1907+097 (Schwartz et al. 1980). It is therefore conceivable that the EXOSAT sources are associated with regions of G043.5$`+`$00.6, but we require additional X-ray imaging to confirm this.

We can estimate the distance to G043.5$`+`$00.6 very roughly using the updated $`\mathrm{\Sigma }`$-$`D`$ relation (Case & Bhattacharya 1998). Assuming a spectral index $`\alpha =0.5`$, appropriate for many shell-type SNRs (e.g. Case & Bhattacharya 1998), we find $`\mathrm{\Sigma }_{1\text{ GHz}}=3.8\times 10^{21}\text{ W m}^2\text{ Hz}^1\text{ ster}^1`$, which implies a linear diameter $`D40`$ pc and distance $`d8.1`$ kpc. We note that these numbers are extremely imprecise (given both intrinsic uncertainties in the $`\mathrm{\Sigma }`$-$`D`$ relation and in estimates of the 1-GHz surface brightness) probably having errors approaching 100%. However, this value is consistent with G043.5$`+`$00.6 being near the tangent point of spiral arm 3 of the Taylor & Cordes (1993) electron density model.

Also intriguing are the objects in Figure 5 marked G042.0$``$00.1 and G041.5+00.4 (see Figure 11). These objects have morphologies that are less well-determined than that of G043.5$`+`$00.6 and are further out in the VLA primary beam, so we cannot make very useful measurements, but their shapes are suggestive of either SNRs or H II regions. G042.0$``$00.1 is centered at $`\alpha _{2000}=19^\mathrm{h}08^\mathrm{m}11^\mathrm{s}`$, $`\delta _{2000}=+08\mathrm{°}00.5\mathrm{}`$, is roughly circular with radius $`4.3\mathrm{}`$, and has a 332-MHz flux of $`1.8\pm 0.4`$ Jy (this value has been corrected for the attenuation of the primary beam of the VLA antennas and is therefore highly uncertain). The protruding region to the south-east of G042.0$``$00.1 was identified as an H II region on the basis of radio recombination lines by Lockman, Pisano, & Howard (1996), but their data do not mention the rest of the complex and their angular resolution should be able to distinguish between the center of G042.0$``$00.1 and the smaller region. G042.0$``$00.1 is present in the NVSS with similar morphology to the 332-MHz image and flux $`S_{\mathrm{NVSS}}=0.3\pm 0.1`$ Jy, suggesting a non-thermal spectrum and leading to a provisional classification as an SNR.

G041.5+00.4 is centered at $`\alpha _{2000}=19^\mathrm{h}05^\mathrm{m}48^\mathrm{s}`$, $`\delta _{2000}=+07\mathrm{°}46.1\mathrm{}`$. It has a complicated morphology, with a brightened rim to the north-east and a compact core to the south-west ($`\alpha _{2000}=19^\mathrm{h}05^\mathrm{m}37^\mathrm{s}`$, $`\delta _{2000}=+07\mathrm{°}44.8\mathrm{}`$), suggestive of a pulsar wind nebula. The source has an overall radius of $`4.9\mathrm{}`$ and a 332-MHz flux of $`1\pm 1`$ Jy (again, this is a highly uncertain value). Comparison of the overall 332-MHz flux with the NVSS flux of $`1.3\pm 0.2`$ Jy leads to a flat spectrum like that of core-dominated SNR or an H II region. We need more data to conclusively classify all three of these objects; further analysis and observations are in progress.

The wealth of information available in Figure 5, concerning both already known objects and newly discovered ones, is readily apparent. It serves to illustrate how productive low-frequency wide-field imaging of the Galactic plane with high resolution can be (e.g. LaRosa et al. 2000).

We would like to thank J. Lazio for his assistance with the low-frequency data reduction and D. Fox for useful comments. DLK is supported by the Fannie and John Hertz Foundation, SRK by NSF and NASA, and MHvK by a fellowship from the Royal Netherlands Academy of Arts and Sciences. The NASA/ESA Hubble Space Telescope is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract No. NAS5-26555.The VLA is operated by the National Radio Astronomy Observatory, which is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The W.M. Keck Observatory is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. We have made extensive use of the SIMBAD database, and we would like to express our appreciation of the astronomers who maintain this database.


text:
# 1 A Slight Improvement

## 1 A Slight Improvement

Consider a finite- dimensional quantum system:

$$\dot{c}=Ac+Bcϵ(t),c(0)=c_0$$
(1.1)
Here $`c`$ evolves on the complex sphere in some $`C^N`$, and $`A`$ and $`B`$ are two skew-Hermitian matrices. It is well known that the set of states reachable from $`c_0`$, via the application of admissible controls, equals the orbit through $`c_0`$ of the group, $`G`$ corresponding to the Lie algebra generated by the matrices $`A`$ and $`B`$, provided this Lie group is compact. See  for instance. This result is valid even if $`G`$ is not compact. This result has been mentioned explicitly by the author in several talks and grant proposals, though the result has been recorded only cryptically in  (primarily because it offers no significant practical improvements). Since such a nominal improvement may, nevertheless, be of interest, a discussion of it is provided below.

To that end, notice that $`A`$ may be assumed to be diagonal, $`A=\text{d}iag(i\lambda _1,i\lambda _2,\mathrm{},i\lambda _N)`$. Writing $`c_k=a_k+ib_k,k=1,\mathrm{},N`$ and viewing the given quantum system as evolving on $`R^{2N}`$, it follows that the drift portion of the system, i.e., the vector field $`Ac`$, is a classical Hamiltonian system with Hamiltonian, $`H=_{k=1}^N\lambda _k(a_k^2+b_k^2)`$ (under the standard symplectic structure of $`R^{2N}`$).

Now the quantum system evolves on the sphere $`S=\{_{i=1}^Nc_i^2=1\}`$. $`S`$ is of course also the set $`\{_{i=1}^Na_i^2+b_i^2=1\}`$. Thus $`S`$ is a bounded set and further the orbit of $`G`$ through the initial condition in $`S`$ is also bounded, even if it is not compact. Thus, the quantum system can be viewed as a nonlinear control system on a bounded manifold, whose drift is a classical Hamiltonian system. Since this implies that the drift is weakly positively Poisson stable, the desired conclusion follows from (see also ). Indeed, these articles show that, under the weak positively Poisson stability of the drift, the reachable set from the initial condition is exactly the leaf of the foliation associated to the drift and control vector fields. But this leaf is nothing but the orbit of $`G`$ through the initial condition. It also implies that the reachable set of the associated unitary generator system from the identity matrix is $`G`$ (cf., ), since if this is not the case, the reachable set of the quantum system could not be the orbit of $`G`$ through the initial condition.

In the opinion of the author, this improvement is of mathematical interest only. The result is of relevance in the situation that the quantum system is not completely controllable, i.e., when the orbit of $`G`$ is not all of $`S`$. However, for the result to have practical value it is desirable to describe the orbit explicitly. One way to do this is to use invariant theory and describe the orbit as a level set. However, invariant theory is tractable only for compact and reductive groups. The author is not aware of any practical situation in which $`G`$ is not compact but is reductive.