text:
# Future supernova probes of quintessence

## I Introduction

Type Ia supernovae (SNe Ia) can be used as standard candles to infer the luminosity distance, $`d_L`$, as a function of redshift , and such data provide a key element in the case for cosmic acceleration. The measurements provide a simple way to estimate cosmological parameters. This is the aim of at least two collaborations: the Supernova Cosmology Project  and the High-redshift Supernova Search Team . However, the analysis of some $`50`$ supernovae at redshifts around $`z0.5`$ has not yet provided strong constraints on the nature of the “dark energy”.

“Quintessence”  has been proposed as a candidate for the “dark energy” to provide a dynamical solution to the cosmological constant problem . The supernova data, although sufficient to constrain the parameters of quintessence models, is not yet abundant enough to allow a discrimination between competing quintessence models. Moreover, the unexpected dimming of the type IA supernovae at redshifts $`z0.5`$ is a phenomenon which could be readily attributed either to an accelerating Universe or to an unexpected luminosity evolution . The present data alone is not sufficient to allow a discrimination between these two possibilities with complete confidence .

Very recently the identification of a SN Ia event at redshift $`z1.7`$  has provided tantalizing evidence that the expansion of the universe could have been decelerating at that epoch. Naturally, no firm conclusions can be drawn from a single event. However, if this result holds up it would rule out the simplest models of source evolution. Furthermore, supernovae events at such redshifts are the type of observations which should enable a discrimination to be made between some of the different quintessence models.

Much more data is needed to enable an accurate estimation of the nature of the dark energy. This might be accomplished by a dedicated space telescope, the SuperNova Acceleration Probe (SNAP) , which aims to collect a large number of supernovae with $`z<2`$. In this paper, we assess the ability of the SNAP mission to determine various properties of the “dark energy”. By analysing a simulated data set, as might be obtained by the proposed SNAP satellite, we can test the ability of such experiments to distinguish among currently attractive quintessence models .

The feasibility of determining the properties of the dark energy component by using simulated data sets has already been considered by several authors . One common approach is to assume that the quintessence field is described by a perfect fluid with equation of state $`P=w\rho `$, where $`w`$ is approximately constant over epochs of interest, or else slowly varying. For example, various authors  consider models with an equation of state linear in redshift, $`w(z)=w_0+w_1z`$. However, many realistic cosmological models could fall outside the confines of these approximations. In this paper we wish to take a different approach by considering models in which the quintessence field is directly given by a Lagrangian, and in which $`w(z)`$ has the freedom to vary widely over measurable redshifts. We will fit simulated SNAP data sets to the exact $`d_L(z)`$ of different quintessence models derived from numerical integration of the coupled Einstein-scalar field equations.

For the fiducial cosmological model used in the simulation of the supernova data set, we consider a form of quintessence, an ultra-light pseudo Nambu–Goldstone boson (PNGB)  which is still relaxing to its vacuum state. Our reasons for this choice of quintessence model are twofold. Firstly, from the viewpoint of quantum field theory PNGB models are the simplest way to have naturally ultra-low mass, spin-$`0`$ particles and hence perhaps the most natural candidate for a presently-existing minimally-coupled scalar field. The effective potential of a PNGB field $`\varphi `$ can be taken to be of the form

$$V(\varphi )=M^4[\mathrm{cos}(\varphi /f)+1],$$
(1)
where the constant term ensures that the vacuum energy vanishes at the minimum of the potential. This potential is characterized by two mass scales, a purely spontaneous symmetry breaking scale $`f`$ and an explicit symmetry breaking scale $`M`$.

The second motivation for choosing a PNGB model rather than other forms of quintessence is that it provides a natural framework for studying the question of the possibility of source evolution.

Analyses of the SN Ia data performed with source evolution have been undertaken in the case of open Friedmann–Robertson–Walker (FRW) models . However, the measurement of the positions of the first acoustic peaks in the angular power spectrum of cosmic microwave background radiation anisotropies by the BOOMERANG, MAXIMA and DASI experiments now gives unequivocal evidence that the Universe is close to being spatially flat . This makes the choice of a spatially open FRW model uncompelling.

On the other hand, the PNGB cosmologies have the virtue that while they are spatially flat, there is no a priori preference for either accelerated or decelerated expansion . Both possibilities are available at modest redshifts, depending on parameter values. Ultimately, the scalar field will undergo coherent oscillations about its minimum, and the resulting luminosity distance will become indistinguishable from that of an Einstein–de Sitter model at late times, although the fraction of energy density in clumped matter, $`\mathrm{\Omega }_m`$, and the fraction of energy density in quintessence, $`\mathrm{\Omega }_\varphi `$, can take any values consistent with $`\mathrm{\Omega }_m+\mathrm{\Omega }_\varphi =1`$. Whether the scalar field is currently rolling down the potential for the first time in the history of the universe – leading to a luminosity distance relation similar to that produced by a cosmological constant – or whether it has already undergone one or more oscillations by the present epoch, is a matter of a choice of initial conditions. Beyond some bounds set by primordial nucleosynthesis these initial conditions are not very much constrained by our present knowledge of the models, resulting in diverse possibilities for cosmological evolution. The current and final values of $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\varphi `$ likewise depend on initial conditions for the scalar field, and on the parameters $`M`$ and $`f`$.

## II accelerating Universe vs luminosity evolution

In an earlier paper , we considered the observational constraints arising from SNe Ia data and gravitational lensing data on cosmological models based on Einstein gravity minimally coupled to a scalar quintessence field with a PNGB potential (1). In the case of the supernovae, we studied the constraints on the PNGB parameters $`M`$ and $`f`$ from 60 supernovae data from the Supernova Cosmology Project (hereafter P98) of Perlmutter et al. . We numerically evolved the coupled Einstein-scalar field equations of motion forward from the epoch of matter-radiation equality to obtain the luminosity distance – redshift relation, comparing the results with the data both with and without an allowance for the possibility of source evolution.

In allowing for source evolution we are acknowledging that the peak luminosities of distant SNe Ia have been normalized according to empirical “Phillips relations”  between observed peak luminosity and supernova decay time, which have been found to be valid at low redshifts. Although we would hope that such relations remain applicable at high redshifts, until the Phillips relations can be modelled and accounted for physically some doubt will always remain about applying these relations at higher redshifts.

A simple empirical model for possible source evolution was employed in ref. : following Drell, Loredo and Wassermann  we added a term $`\beta \mathrm{ln}(1+z)`$ to the distance modulus. This particular luminosity evolution function is simply chosen as an illustrative example, and is not singled out by any particular physical model. Naturally, one can criticize it on these grounds. However, the Phillips relations are also purely empirical, and the purpose of our study as with that of ref.  was simply to test the extent to which any form of source evolution was able to account for the observed data as compared to an accelerated expansion.

For the purposes of analysing the data it was assumed in refs.  that the prior for the parameter $`\beta `$ was a Gaussian with a mean $`\beta _0`$ and standard deviation $`b`$. The best–fit value of $`\beta _0`$ for the P98 data set was found to vary very slightly with initial conditions, taking a value $`\beta _0=0.414`$ for $`w_i\varphi _i/f=1.5`$ and $`\beta _0=0.435`$ for $`w_i=0.2`$, where $`\varphi _i`$ is the value of the quintessence field at the time of onset of matter domination ($`z=1100`$). In this paper we use the initial value $`w_i=0.2`$. Our reason for doing this is that in general, and especially for larger values of $`w_i`$, there is a tension between the values of $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ favoured by current observations, with larger values of $`H_0t_0`$ often corresponding to values of $`\mathrm{\Omega }_{\varphi 0}`$ rather greater than 2/3. It is still possible to find parameter values of $`M`$ and $`f`$ which give $`\mathrm{\Omega }_{\varphi 0}0.7`$ while yielding an age for the Universe consistent with the lower bound of the current observational range , even for $`w_i=1.5`$. However, one finds that the tension between the values of $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ is more easily mitigated by choosing a lower value of $`w_i`$, for which the scalar field spends more of its early dynamical history higher up the potential hill. Confidence limits on the $`M,f`$ parameter space obtained from the P98 data set, assuming no evolution, are shown in Fig. 1.

In refs.  we included all 60 data from the P98 data set. In fact, in ref.  Perlmutter et al. considered a subset of 54 SNe Ia which excluded 6 supernovae events: two that are the most significant outliers from the average light–curve width, two that have the largest differences between observed and expected magnitudes (or fluxes), and two that are likely to be reddened. It was discovered that excluding these six supernovae produced a more robust fit of the cosmological parameters $`\mathrm{\Omega }_{m0}`$ and $`\mathrm{\Omega }_{\mathrm{\Lambda }0}`$. In refs.  we took a conservative approach and included all 60 SNe Ia, since we are investigating models with different cosmological parameters from those fitted in ref.  and we did not wish to prejudge matters. However, one may of course redo the analysis of refs.  for the “fit C” data set of Perlmutter et al. For the PNGB model we find that in the best–fit case this reduces the normalized $`\chi ^2`$ parameter from 101.6, i.e., 1.75 per degree of freedom, for the full data set to 58.65, or 1.13 per degree of freedom, for the reduced “fit C” data set. For the models with source evolution the normalized $`\chi ^2`$ is similarly reduced from 101.1, or 1.74 per degree of freedom, to 58.21, or 1.12 per degree of freedom. Thus the “fit C” data provides a more robust fit in all cases we have studied, whatever the model assumptions.

Since the results in ref.  were based on the full P98 data set, we have included results based on both the full data set and the reduced data set in Fig. 1, to compare the differences. Note that although the 2$`\sigma `$ allowed area to the right disappears entirely for the reduced data set for $`w=0.2`$, this is not the case for $`w=1.5`$, and thus the discrepancy between the P98 data set and that of Riess et al.  which was noted in Ref.  remains for the PNGB model whether or not one uses the full or reduced data set.

Here we will take a somewhat different approach to ref.  in our treatment of the case with source evolution. We will assume that the parameter $`\beta `$ has a flat prior bounded over some range $`\mathrm{\Delta }\beta =2.3`$, with limits corresponding to $`0.6<\beta <1.7`$, and marginalize the likelihood function over this prior. The choice of the bounds on $`\beta `$ in explained below in section IIB.

Despite the differences from the approach used in Ref.  the conclusions of the analysis are broadly similar, and are shown in Fig. 2. In particular, two regions of parameter space appear to be singled out: region I (see Fig. 1), which corresponds to cosmological models for which the scalar field is rolling down the potential (1) for the first time (with $`\varphi `$ increasing) at the present epoch, and region II (see Fig. 2) in which the scalar field is rolling down the potential for the second time (with $`\varphi `$ decreasing) at the present epoch. In the absence of any evolution the P98 data set favours Region I, whereas if evolution is allowed for then region II appears to be favoured slightly more than region I.

Note that over the entire parameter space considered, we find the value of $`\widehat{\beta }`$ with the best fit of all is $`\widehat{\beta }0.622`$ using the entire P98 data set, which occurs at parameter values $`M=4.55\times 10^3h^{1/2}`$eV, $`f=0.884\times 10^{18}`$GeV, while for the “fit C” data set it is $`\widehat{\beta }0.720`$, which occurs at $`M=4.64\times 10^3h^{1/2}`$eV, $`f=0.901\times 10^{18}`$GeV. Averaged over the $`(M,f)`$ parameter space we find $`\widehat{\beta }0.286`$ for the full P98 data set, and $`\widehat{\beta }0.348`$ for the “fit C” data set. Thus in both cases the value of $`\widehat{\beta }`$ is peaked in region II, but for the “fit C” data the overall values of $`\widehat{\beta }`$ are somewhat larger.

The results of refs.  were inconclusive. Ostensibly the source evolution models appeared to provide a slightly better fit. However, the extent to which they were better was not quantified statistically. One relatively straightforward way to compare the relative strengths of rival theoretical models in fitting a common set of data is by the calculation of Bayes factors , even when the models involve different numbers of parameters. Given more parameters, one may always find a better fit to the data – however, the Bayes factor approach includes an effective “Ockham’s razor factor” which adds a weighting against the inclusion of extra parameters. The Bayes factor method may be readily applied to the results of  to give an explicit statistic to compare the relative strengths of a cosmological explanation of the SNe Ia data versus an explanation in terms of source evolution.

In Bayesian inference, when comparing rival models $`M_i`$, each with parameters $`𝒬_i`$, the likelihood for a model conditional on the data ($`D`$) in a model comparison calculation is equal to the average likelihood for its parameters,

$$p(D|M_i)=𝑑𝒬_i\text{ }p(𝒬_i|M_i)(𝒬_i),$$
(2)
where $`p(𝒬_i|M_i)`$ is the prior probability for $`𝒬_i`$, and $`(𝒬_i)`$ is the sampling probability for $`D`$ presuming $`M_i`$ to be true. The ratio of model likelihoods,

$$B_{ij}=\frac{p(D|M_i)}{p(D|M_j)}$$
(3)
is called the Bayes factor. When the prior odds does not strongly favour one model over another, the Bayes factor can be interpreted just as one would interpret an odds in betting. Kass and Raftery  provide a comprehensive review of Bayes factors, and the recommended interpretation is summarized in Table 1.

### A P98 data set

In order to ascertain what degree of improvement could be expected with data from SNAP, we will begin by calculating Bayes factors for the existing P98 data set.

The prior ranges for parameters play an important role in Bayesian model comparison. We will make a choice of the prior ranges of parameters in the different models by the resulting values of $`\mathrm{\Omega }_{m0}`$ and $`H_0t_0`$ that they give rise to. In particular, we will choose a conservative bound of $`0.5\mathrm{\Omega }_{\varphi 0}0.9`$. For the expansion age of the universe we choose $`0.72H_0t_01.15`$: for $`H_0=70\text{ km s}^1\text{Mpc}^1`$ it corresponds to a conservative bound of $`t_013\pm 3`$ Gyr consistent with recent estimates . The parameter space bounded by these two constraints is indicated by the contours of $`\mathrm{\Omega }_{m0}`$ and $`H_0t_0`$ in Figs. 1 and 2. In the large $`f`$ region in region I where the contours appear to be parallel to the $`f`$–axis, a cut-off at $`f10^{19}`$ GeV, the Planck scale, is chosen.

We calculate the average parameter likelihood by integrating $`(M,f)`$ over the prior ranges for the parameter and dividing the integral by the relevant area. For the full P98 data set the Bayes factor for the model without luminosity evolution versus the model with luminosity evolution, which is the ratio of the average parameter likelihoods, is $`B2.6`$. For the reduced “fit C” data set we obtain a Bayes factor of 2.5, which is the same as far as its interpretation is concerned. Thus the data alone cannot discriminate between the two hypotheses.

### B SNAP data sets

We will now simulate a data set which would be expected to be obtained by the SNAP satellite to investigate its potential for discriminating between the two possibilities. We will simulate two data sets from fiducial cosmologies chosen with parameters in each of regions I and II. In each case we will introduce a random error to the exact distance moduli to simulate a future supernova data set that has been converted to a table of effective magnitudes $`m_i`$ and redshifts $`z_i`$ of objects with a single fiducial absolute magnitude $`_0`$. We will consider both statistical and systematic uncertainties in the magnitudes. Typically the redshifts are known to sufficiently high precision that their uncertainties can be ignored. We assume the supernovae observed uniformly within four different redshift ranges with the following different sampling rates, which are the same as those chosen by Weller and Albrecht : In the first range from $`z=0`$$`0.2`$ we assume that there are $`50`$ observations, in the second and largest redshift range from $`z=0.2`$$`1.2`$ there are $`1800`$ SNe and in the two high redshift bins, $`z=1.2`$$`1.4`$ and $`z=1.4`$$`1.7`$, there are $`50`$ SNe and $`15`$ SNe observations respectively. The statistical error in magnitude is assumed to be $`\sigma _{\text{mag}}=0.15`$, including both measurement error and any residual intrinsic dispersion after calibration.

For data set A, we assume there is no luminosity evolution. We choose the fiducial parameters $`M=2.53\times 10^3h^{1/2}`$eV and $`f=0.58\times 10^{18}`$GeV; these parameters give $`(\mathrm{\Omega }_{m0},\mathrm{\Omega }_{\varphi 0})(0.23,0.77)`$ and $`H_0t_00.9`$. For data set B, we assume there is a luminosity evolution and we add a term $`\beta \mathrm{ln}(1+z)`$ to the distance moduli. We choose the fiducial parameters $`M=4\times 10^3h^{1/2}`$eV and $`f=0.676\times 10^{18}`$GeV. These parameters give $`(\mathrm{\Omega }_{m0},\mathrm{\Omega }_{\varphi 0})(0.25,0.75)`$ and $`H_0t_00.8`$. The best–fit value of the parameter $`\beta `$ for this particular choice of $`M`$ and $`f`$ is $`\widehat{\beta }=0.659`$ for the full data set, and $`\widehat{\beta }=0.727`$ for the “fit C” reduced data set. We will take $`\beta 0.659`$ as the fiducial parameter for data set B, as it is the somewhat more conservative value, representing less luminosity evolution.

The fiducial parameters for the two data sets have been deliberately chosen so as to give similar values of $`\mathrm{\Omega }_{m0},\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ on one hand, so as to be consistent with values favoured by other observational tests. On the other hand, the fiducial parameters for data set A are centred in region I, which corresponds to a universe in which the scalar field is rolling down the potential for the first time at the present epoch, whereas the fiducial parameters for data set B are centred in region II, which corresponds to a universe in which the scalar field is rolling down the potential for the second time at the present epoch. As was discussed in  region I is favoured if the SNe Ia luminosity distances are true cosmological distances, whereas region II becomes significantly favoured if the data hides a simple luminosity evolution.

We will use Bayesian analysis to estimate uncertainties. Rigorous calculation of the likelihood for the quintessential parameters $`𝒬`$ is very complicated, requiring the introduction and estimation of many additional parameters, including parameters from the light–curve model and parameters for characteristics of the individual supernovae. With several simplifying assumptions , the final result is relatively simple. One finds

$$(𝒬)\mathrm{e}^{\chi ^2/2},$$
(4)
where

$$\chi ^2(𝒬)=\underset{i}{}\frac{[\mu _{s,i}\mu (z_i;𝒬)]^2}{\sigma _{\mu ,i}^2}.$$
(5)
In the above equation,

$$\mu (z_i;𝒬)=5\mathrm{log}d_L(z_i;𝒬)+25$$
(6)
is the distance modulus predicted by each model with parameters $`𝒬`$, while $`\mu _{s,i}=m_i_0`$ is the simulated distance modulus, and its uncertainty is $`\sigma _{\mu ,i}=0.15`$. Note that we fix Hubble parameter to $`H_0=70\text{ km/s Mpc}^1`$ to simplify the computation.

The model with an unexpected luminosity evolution corresponds to replacing eq. (5) with

$$\chi ^2(𝒬,\beta )=\underset{i}{}\frac{[\mu _{s,i}\beta \mathrm{ln}(1+z_i)\mu (z_i;𝒬)]^2}{\sigma _{\mu ,i}^2}.$$
(7)
We will marginalize over $`\beta `$ with a flat prior, to obtain the marginal likelihood

$$(𝒬)=\frac{1}{\mathrm{\Delta }\beta }𝑑\beta \mathrm{e}^{\chi ^2/2}.$$
(8)
The above integration can be performed by an analytic marginalization . We separate $`q`$ from $`\chi ^2`$ where

$`q(𝒬)`$ $`=`$ $`{\displaystyle \frac{\widehat{\beta }^2(𝒬)}{\overline{\sigma }^2}}+{\displaystyle \underset{i}{}}{\displaystyle \frac{h^2(z_i;𝒬)}{\sigma _i^2}}`$ (9)
$`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{[h(z_i;𝒬)\widehat{\beta }(𝒬)\mathrm{ln}(1+z_i)]^2}{\sigma _{\mu ,i}^2}}`$ (10)
is independent of $`\beta `$. The integral is thus an integral over a Gaussian in $`\beta `$ located at $`\widehat{\beta }`$ with standard deviation $`\overline{\sigma }`$. $`\widehat{\beta }`$ is the best–fit value of $`\beta `$ given $`𝒬`$, and $`\overline{\sigma }`$ is its conditional uncertainty. They are given by

$`{\displaystyle \frac{1}{\overline{\sigma }^2}}`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{[\mathrm{ln}(1+z_i)]^2}{\sigma _{\mu ,i}^2}},`$ (11)
$`\widehat{\beta }(𝒬)`$ $`=`$ $`\overline{\sigma }^2\left[{\displaystyle \underset{i}{}}{\displaystyle \frac{h(z_i;𝒬)\mathrm{ln}(1+z_i)}{\sigma _{\mu ,i}^2}}\right],`$ (12)
where

$$h(z_i;𝒬)=\mu _{s,i}\mu (z_i;𝒬).$$
(13)
As long as $`\beta `$ is inside the prior range and $`\overline{\sigma }\mathrm{\Delta }\beta `$, the value of the integral is well approximated by $`\overline{\sigma }\sqrt{2\pi }`$, so that

$$(𝒬)=\frac{\overline{\sigma }\sqrt{2\pi }}{\mathrm{\Delta }\beta }\mathrm{e}^{q/2}.$$
(14)
The question of what bounds to place on the prior range for the parameter $`\beta `$ is an interesting one, since we are dealing with a purely empirical model of source evolution with no a priori restrictions on $`\beta `$, other than that very large values of $`\beta `$ would not be physically plausible. We will therefore seek the narrowest range of values of $`\beta `$ consistent with our numerical results. One condition on the prior range for $`\beta `$ is that $`\widehat{\beta }`$ also lies inside the prior range: setting too tight a bound on $`\mathrm{\Delta }\beta `$ would count against parameter values of $`M`$ and $`f`$ which might otherwise be included. By explicit numerical integration we have found that for the initial conditions and the range of parameter values of $`M`$ and $`f`$ we have considered $`\widehat{\beta }`$ lies in the range $`0.6<\widehat{\beta }<1.7`$. We will therefore take these bounds to be the prior range for $`\beta `$.

If we fit the data set obtained from fiducial model A assuming that there is no luminosity evolution, the parameters $`M`$ and $`f`$ are well constrained by the simulated data. The 95.4% confidence level contour bounds a very small region around the best–fit values of $`M2.55\times 10^3h^{1/2}`$eV and $`f0.592\times 10^{18}`$GeV, with $`\chi ^2=1910`$. These best–fit parameters are very close to the fiducial parameters that generate the data set. The average parameter likelihood is $`_{\text{ave}}7.1\times 10^4_0`$, where $`_0=\mathrm{exp}(1910/2)`$. If luminosity evolution is considered to exist (see Fig. 3), fitting the data set from fiducial model A we find that the 95.4% confidence level contour includes a larger region in region I of the parameter space around the fiducial parameters, and a region in region II. The average parameter likelihood is $`_{\text{ave}}1.9\times 10^3_0`$. Hence, the Bayes factor for the model with luminosity evolution versus the model without luminosity evolution is $`B2.7`$. We conclude that it is not possible to discriminate between the two hypotheses even though the underlying simulated data does not have any luminosity evolution.

If we fit the data set obtained from the fiducial model B by assuming that there is no luminosity evolution, then the 95.4% level bounds a very small region around the best–fit values of $`M2.27\times 10^3h^{1/2}`$eV and $`f0.764\times 10^{18}`$GeV. The average parameter likelihood is $`_{\text{ave}}9.4\times 10^{18}_0`$, where $`_0=\mathrm{exp}(1924/2)`$. If luminosity evolution is considered to exist (see Fig. 4), on the other hand, then the 95.4% level bounds a small region around the best–fit values of $`M4.00\times 10^3h^{1/2}`$eV and $`f0.695\times 10^{18}`$GeV, with $`\widehat{\beta }0.626`$ and normalized $`\chi ^2=1924`$. The average parameter likelihood is $`_{\text{ave}}8.3\times 10^4_0`$. This time, the Bayes factor for the model with luminosity evolution relative to the model without luminosity evolution is much greater than one, $`B9.1\times 10^{13}`$. Thus if the true data were derived from fiducial model B, it would show up clearly as giving a $`d_L(z)`$ which can be discriminated from the model without luminosity evolution.

This result in no way contradicts the inconclusive result obtained from data set A. If the underlying data were to truly follow a simple luminosity evolution, then the analysis above shows that this would stand out as a very strong positive signal in the Bayes factor test. If the underlying data has no luminosity evolution on the other hand, then it may be possible to obtain a good fit with an extra luminosity evolution parameter, simply because one has an extra parameter to fit. Thus a simple luminosity evolution is easy to rule in, but difficult to rule out. Combined together the two results demonstrate the efficacy of the Bayes factor approach. With data set B we would obtain a decisive result, but by Ockham’s razor we should reject the more complicated evolutionary hypothesis if we were to obtain an inconclusive result such as that pertaining to data set A. In fact, we have found that by narrowing the range of the prior $`\mathrm{\Delta }\beta `$ that the Bayes factor for data set A can even be made slightly greater than the cut-off value of 3 listed in Table 1. Clearly any weakly positive results should therefore also be treated with some caution.

## III PNGB versus other potentials

The constraints on cosmological parameters which arise from type Ia supernovae have been broadly studied for a wide range of quintessence models . However, all these models seem to have ranges of parameters that fit the existing data equally well, with no particular model standing out as being observationally preferred. The lack of data and the large statistical error bars are among the contributing reasons. In this section, we will investigate the potential of the SNAP supernovae data set to discriminate between competing quintessence models, as contrasted with the P98 data set, assuming no source evolution.

We consider some simple scalar field potential functions which have been widely studied in the literature. They are the simple exponential potential

$$V=V_A\mathrm{e}^{\lambda \kappa \varphi },$$
(15)
the inverse power–law potential

$$V=V_A\varphi ^\alpha ,$$
(16)
and the double–exponential potential

$$V=V_A\mathrm{exp}(A\mathrm{e}^{\sqrt{2}\kappa \varphi }),$$
(17)
where $`\kappa `$ is related to Newton’s constant by $`\kappa ^2=8\pi G`$. We evaluate the luminosity distance – redshift relation by numerically evolving the Einstein-scalar field evolution equations forward from the time of onset of matter domination, as what we did for the PNGB models. For the simple exponential and the double–exponential potentials, without lost of generality we choose the initial condition $`\varphi _i=0`$. For the power–law potential, we begin the integration by assuming that the scalar field has already reached the tracker solution. These leave us with two parameters in each case. In Figs. 5, 6, and 7, we display the confidence limits for the 54 supernovae Ia in the P98 “fit C” reduced data set, on the parameter spaces for these quintessence models. (The equivalent figures for the full data set have been omitted, as the 1$`\sigma `$ and 2$`\sigma `$ confidence regions are very similar to those shown in Figs. 57, but just slightly narrower.)

In this section we will only compare models (15)–(17) on the basis of assuming that the true data follows the PNGB model, with the fiducial parameters of data set A. Naturally, we could assume any of the other potentials (15)–(17) as our fiducial model, and then compare the other models on the basis of a new fiducial data set. This would allow us, for example, to compare the inverse power–law potential (16) with the simple– and double–exponential potentials (15), (17), a test which we have not performed here. However, we will restrict our attention to the fiducial data set A relevant to the PNGB potential (1), since the amount of computer time involved in computing the fiducial data sets is large, and an analysis based on a fiducial PNGB model will prove to be sufficient to illustrate the dramatically increased discriminatory power of a SNAP dataset as opposed to presently available SNe Ia data. Naturally, one could extend the discussion to other fiducial data sets, but we will leave that to future work.

### A P98 data set

Firstly we will calculate the Bayes factor of the P98 data set to compare the PNGB model with all the other quintessence models. In order to calculate Bayes factors which compare different theoretical models we must naturally make choices of the range of prior parameters integrated over, and we are dealing with different parameters in the different models. Different choices of priors for these parameters would lead to some variation in the resulting Bayes factors. We will make a choice of the range of prior parameters by using the resulting values of $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ that they give rise to. However, the parameter spaces are not completely bounded by the $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ constraints. Without a physical cut-off for the parameters, such as a cut-off for $`f`$ at Planck scale for the PNGB model, we will investigate in each case the dependence of the Bayes factor comparison on the prior ranges of the parameters. We will show that in some cases, the choice of the prior parameter ranges do not affect the conclusion from the Bayes factor analysis.

For the simple exponential potential model, as shown in Fig. 5, the $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ contours diverge at large $`V_A`$, putting a bound on the parameter space. In the small $`\lambda `$ region, where the potential is effectively a cosmological constant model whose value depends on $`V_A`$ only, the $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ contours remain near parallel. For the full P98 data set the average parameter likelihood over the region bounded by the $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ constraints, with a cut-off at $`\lambda =0.1`$, is $`_{\text{ave}}0.67_R`$, where $`_R`$ is the average parameter likelihood for the PNGB model. In the small $`\lambda `$ region the likelihood only changes slightly with decreasing $`\lambda `$ and tends to a constant, corresponding to an average parameter likelihood which we estimate to be $`_{\text{ave}}1.5_R`$. Since the small $`\lambda `$ region extends to an infinite range in terms of the parameterization shown in Fig. 5, if we enlarge the prior range of $`\lambda `$ to include smaller values, then it has the effect of bringing the average parameter likelihood closer to the average parameter likelihood in the small $`\lambda `$ region, which is effectively an upper bound. Thus the Bayes factor for the simple exponential potential model versus the PNGB model lies in the range $`0.67<B<1.5`$ for the full P98 data set. Similarly $`0.65<B<1.5`$ for the reduced “fit C” data set. In each case the larger value is the one appropriate to including the entire small parameter range for $`\lambda `$. Such a Bayes factor is too small to confidently discriminate between the two models.

For the double–exponential potential model, the average parameter likelihood over the region bounded by the $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ constraints, with a cut-off at $`A=0.01`$, is $`_{\text{ave}}0.92_R`$ for the full P98 data set. The average parameter likelihood in the small $`A`$ region, where the potential is effectively a cosmological constant model whose value depends on $`V_A`$ only, is $`_{\text{ave}}1.5_R`$. By a similar argument to above, the Bayes factor for the double–exponential potential model versus the PNGB model varies slightly over a range $`0.92<B<1.5`$ for the full P98 data set. Similarly, $`0.38<B<1.5`$ for the reduced “fit C” data set. Again, the Bayes factor is too small to confidently discriminate between the two models.

For the inverse power–law potential model, as shown in Fig. 6, the $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ contours remain nearly parallel towards the upper right-hand region where both $`V_A`$ and $`\alpha `$ are large. The average parameter likelihood over the region bounded by the $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ constraints with arbitrary cut-offs at $`\kappa ^{\alpha +2}V_A/H_0^2=10^{15}`$ and $`\alpha =15`$ is $`_{\text{ave}}0.14_R`$ for the full P98 data set, and $`_{\text{ave}}0.11_{RC}`$ for the reduced “fit C”data set, where $`_{RC}`$ is the average parameter likelihood for the PNGB model based on the “fit C” data set. This gives a Bayes factor for the PNGB model versus the inverse power–law potential model of $`B=7.2`$ for the full P98 data set, or $`B=9.2`$ for the “fit C” data set. Both values are large enough to provide slightly positive evidence that the P98 data set favours the PNGB model over the inverse power–law potential model. Since the Bayes factor is only weakly positive, however, we cannot not place strong confidence in this conclusion. Note that the $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ contours both extend towards the smaller likelihoods region. Therefore, increasing the prior ranges for $`V_A`$ and $`\alpha `$ will decrease the average parameter likelihood, and increase the slightly positive preference for the PNGB potential as compared to the inverse power–law potential in the Bayes factor test.

### B SNAP data set

We have shown that the P98 data set does not particularly favour any of the PNGB, simple exponential potential, and double–exponential potential models relative to the others. The P98 data set does disfavour the inverse power–law potential model as compared with the PNGB model. This might be considered to result from the fact that the $`2\sigma `$–confidence region only overlaps with a small region (the low $`\alpha `$ region) of the whole parameter space allowed by the priors set on $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$.

In this section we want to compare the typical results we should expect when a future SNAP data set is used. We will use data set A, simulated from a fiducial PNGB model assuming no luminosity evolution. We will study the potential of this data set to distinguish the PNGB model from other quintessence models.

For the simple exponential potential model, the average parameter likelihood over the region bounded by the $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ constraints, with a cut-off at $`\lambda =0.1`$, is $`_{\text{ave}}2.1\times 10^5_0`$, where $`_0=\mathrm{exp}(1910/2)`$. In the small $`\lambda `$ region, we estimate the average parameter likelihood to be $`_{\text{ave}}6.6\times 10^{31}_0`$. Comparing with the average parameter likelihood for the PNGB model, $`_{\text{ave}}7.1\times 10^4_0`$, the Bayes factor for the PNGB model versus the simple exponential potential model takes values in the range $`34`$$`10^{27}`$, depending on the prior range for $`\lambda `$. This would provide a very strong conclusion that data set A favours the PNGB model over the simple exponential potential model.

For the double–exponential potential model, the average parameter likelihood over the region bounded by the $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ constraints, with a cut-off at $`A=0.01`$, is $`_{\text{ave}}1.6\times 10^3_0`$. The average parameter likelihood in the small $`A`$ region, is $`_{\text{ave}}2.1\times 10^{29}_0`$. The Bayes factor for the PNGB model versus the double–exponential potential model lies in the range $`0.44`$$`10^{25}`$, depending on the the prior range for $`A`$. The lower bound of the Bayes factor, $`B=0.44`$, is too small to confidently discriminate between the two models. However, provided small values of the parameter $`A`$ are included in the prior range, the Bayes factor would provide a strong conclusion that data set A favours the PNGB model over the double–exponential potential model.

For the inverse power–law potential model, the average parameter likelihood over the region bounded by the $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ constraints with arbitrary cut-offs at $`\kappa ^{\alpha +2}V_A/H_0^2=10^{15}`$ and $`\alpha =15`$ is $`_{\text{ave}}1.9\times 10^{10}_0`$. This gives a Bayes factor for the PNGB model versus the inverse power–law potential model of $`B10^6`$. This provides very strong evidence that data set A would favour the PNGB model over the inverse power–law potential model.

## IV Discussion

Let us now consider our results in relation to previous work concerning the feasibility of determining the properties of the dark energy from future supernovae surveys . One approach in past studies of scalar field quintessence has been to assume a potential $`V(\varphi )`$ over which the scalar field, $`\varphi `$, would have slowly varied during cosmological time scales, and then test the efficacy of reconstructing the potential . Another approach has been to assume that the quintessence field can be described by a perfect fluid with slowly varying equation of state $`P=w(z)\rho `$, expand $`w(z)`$ in a power series, and then test the efficacy of determining the coefficients in the power series .

The conclusions of investigations to date are mixed . Weller and Albrecht  find that many models can be distinguished with a fit to a linear equation of state for the dark energy, $`P=w(z)\rho `$ with $`w(z)=w_0+w_1z`$, but only if the current mass density, $`\mathrm{\Omega }_m`$, is known to a high precision. Barger and Marfatia  find that, even by putting $`\mathrm{\Omega }_m=0.3`$ exactly, there is still a possibility of obtaining data sets which might not discriminate between quintessence and “$`k`$-essence”, namely an alternative form of dark energy with a scalar field characterized by non-linear kinetic terms . Wang and Garnavich  consider two classes of functions $`w(z)`$, corresponding respectively to a linear variation and to $`k`$-essence. Using somewhat different techniques to other authors they are more optimistic about prospects for determining $`w(z)`$ from future SNe Ia data.

In the present paper we have taken a different approach from those above, by considering a class of models – the PNGB models – which are very well motivated from the point of view of particle physics, but for which the above methods will not always be satisfactory in the case of all plausible parameter ranges, given the potentially oscillatory nature of $`w(z)`$ and the corresponding fact that the scalar field may have varied over a wide range of values of $`V(\varphi )`$ over observable time scales. We fitted the simulated SNAP data sets to the exact $`d_L(z)`$ of different models obtained by numerical integration, and compared them to other models.

One real drawback of all approaches is that as yet there is no preferred physical model for the dark energy. On one hand this means that any approximations made in potential reconstruction methods may be too restrictive, since many different potential energy functions $`V(\varphi )`$ are conceivable, and many of these may give results degenerate with each other. On the other hand, using a given Lagrangian for the quintessence field, as we have done, limits us to a model by model test.

Nevertheless, we find that data sets such as those that would be produced by SNAP promise to be very successful on some tests, even if they will probably be less successful on others. In particular, while existing data is not yet sufficiently large to discriminate between various quintessence models or models with evolution , we have shown that the much larger size and smaller error bars of the simulated SNAP type IA supernovae data sets provide much tighter constraints on the parameters for quintessence models such as those corresponding to pseudo Nambu–Goldstone bosons.

By evaluating Bayes factors in the context of the PNGB model, we have shown that future satellite SNe Ia data sets should have greater success in detecting whether the observed luminosity distance – redshift relation is purely cosmological in origin, or is significantly contaminated by evolutionary effects of the sources. The results of section IIB showed that although it may be difficult to completely rule out luminosity evolution, if the true data were from a population with luminosity evolution then this would provide a strong distinctive signal. We have only studied one simple illustrative supernova luminosity evolution function, but we expect that similar conclusions would apply to other simple luminosity evolution models.

We have further shown that with the future data it should be possible to discriminate PNGB model from some other particular types of quintessence; in particular, it gives a very different signature to that of simple inverse power–law potentials or simple exponential potentials. The case of a double–exponential potential gives a lower Bayes factor, and may therefore be more difficult to distinguish from the PNGB model. However, even in this case some distinction between the two models is possible if one allows a suitably large prior range of the parameter $`A`$, to include small values.

A number of obvious extensions of our analysis are possible. In the case of testing source evolution versus the case with no evolution, for example, it would be interesting to determine by how much we can reduce the parameter $`\beta `$ for the simulated SNAP data set B while still obtaining a very strong result for the discriminatory power of the relevant Bayes factor. Given the magnitude of the value obtained in section IIB, we suspect that the fiducial value of $`\beta `$ could be significantly reduced. Likewise many other tests could be performed with fiducial data sets based on other quintessence potentials.

In conclusion, we find that future supernovae measurements such as those that would be afforded by the SNAP satellite, will have the power to significantly increase our knowledge of the properties of the dark energy in the universe. To be completely confident, however, we will require a deeper theoretical understanding of the nature of the dark energy and hopefully new input from fundamental physics.

## Acknowledgements

DLW wishes to acknowledge the financial support of Australian Research Council grant F6960043.


format: text

text:
# Electrical transport through individual DNA molecules

## Acknowledgments

Support from the Research Grants Council of the Hong Kong Government and the National Natural Science Foundation of China is gratefully acknowledged.


text:
# Conductance scaling at the band center of wide wires with pure non–diagonal disorder

## I Introduction

The existence of peculiar properties that differentiate pure non–diagonal disorder models from disorder models including diagonal disorder has a long history. Probably, the first contribution along this line corresponds to Dyson’s work on a one–dimensional phonon model published in 1953dyson . The existence of a divergent density of states at the band center of such disordered systems implying a divergent localization length was pointed in later works1d\_a ; 1d\_b . Since then many works on this subject have been published (a representative list can be found in Ref.(mudry0, )). The existence of some kind of delocalization transition at the band center of such models certainly disturbs the widely accepted statement saying that the specific form of disorder does not matter in single parameter scaling theoryeconomou .

A bipartite lattice is a lattice that can be divided into two sublattices such that the Hamiltonian changes sign under a transformation that changes the sign of the wavefunction on one sublattice. When pure non–diagonal disorder is considered on a bipartite lattice, the electron-hole symmetry of the spectrum is not destroyed by disorder. This property has important consequences as, for example, the existence of an eigenstate at $`E=0`$ for any disorder realization of a system constituted by an odd number of sites (the spectrum shows $`\pm E_\alpha `$ pairs plus a state at $`E=0`$). Recently, several works pointing to the exotic behavior of transport properties of quantum wires showing chiral symmetry have been publishedmudry0 ; mudry . For example, scaling of the conductance strongly depends on the parity of the number of channels along the wire. Also, related activity in field theory has produced several models demonstrating a delocalization transition in the vicinity of the zero energy statecampos .

In this paper, the scaling properties of the simplest disorder model preserving chiral symmetry have been carefully analyzed at the band center of quantum wires. Numerical simulation has been used to get the whole distribution of conductances. While the mean value of the conductance decays algebraically for wires of odd width, an alternative measure of the central value of distributions -the median- shows an exponential decay. Actually, the larger part of the conductance measurements has an exponentially small value for large wires. Therefore, although the conductance distribution is certainly peculiar at the band center, I would still use the term exponential localization when referring to the scaling behavior of the conductance of long wires. Results change when width and length of the wire coincide. In this case, the numerical simulation presented in this paper uses a dot geometry point of view keeping the size of the contacts finite while dot area scaling proceeds. The conductance distribution converges to a well defined limit that is compatible with the predicted critical behavior of the state at the band center. It is shown that the analytic form used to fit conductance distributions comes from a new underlying wavefunction statistics describing the distribution of weight ratios of the $`E=0`$ wavefunction.

The format of the paper is as follows. Section II defines a quite simple chain model that allows some analytical results and an unbound numerical simulation. Section III gives a more general disorder model on finite rectangular clusters of the square lattice. The way in which the conductance is calculated is presented in Section IV. Numerical results are given in Section V, first for wide disordered wires and, second, for square clusters. The last Section of the paper compiles the main conclusions reached by this numerical study of conductance scaling.

## II Toy model

Let me begin with a detailed description of the scaling properties of the conductance of a chain showing pure non–diagonal disorder. In this case, sign changes of the hopping parameters do not matter and changes in their absolute value should be considered. One case that allows some analytic results consists of a chain with hopping parameters that randomly take values of $`1`$ and $`\sqrt{2}`$ with equal probability. Starting from the selfenergy of an ideal chain at its band center $`\mathrm{\Sigma }_0=i`$, successive selfenergies can be obtained through the disordered part of the chain by means of the usual iterative sequence (see, for example, Eqs. (7)-(10) in Ref.(1d\_b, )):

$$\mathrm{\Sigma }_{n+1}=\frac{t_{n}^{}{}_{}{}^{2}}{\mathrm{\Sigma }_n},$$
(1)
where $`t_{n}^{}{}_{}{}^{2}`$ takes values 1 or 2 with equal probability. Since selfenergy remains purely imaginary, a more convenient form is possible:

$$\mathrm{\Sigma }_{n+1}/i=\frac{t_{n}^{}{}_{}{}^{2}}{(\mathrm{\Sigma }_n/i)}.$$
(2)
The conductance $`G`$ of a sample is obtained from the selfenergy at the end of the disordered chain part $`\mathrm{\Sigma }_N/i`$:

$$G=\frac{4(\mathrm{\Sigma }_0/i)(\mathrm{\Sigma }_N/i)}{[1+(\mathrm{\Sigma }_0/i)(\mathrm{\Sigma }_N/i)]^2}=\frac{4(\mathrm{\Sigma }_N/i)}{[1+(\mathrm{\Sigma }_N/i)]^2}.$$
(3)
Eq. (3) shows that the conductance varies between 0 and 1 as it corresponds to a one–channel system. Actually, the conductance takes the value 1 only if $`\mathrm{\Sigma }_N/i=\mathrm{\Sigma }_0/i=1`$. Repeated use of Eq.(2) shows that the form in which the random hopping elements appear in the selfenergy expression is:

$$\frac{t_{0}^{}{}_{}{}^{2}t_{2}^{}{}_{}{}^{2}t_{4}^{}{}_{}{}^{2}\mathrm{}}{t_{1}^{}{}_{}{}^{2}t_{3}^{}{}_{}{}^{2}t_{5}^{}{}_{}{}^{2}\mathrm{}}.$$
(4)
Therefore, a perfect transmission through the chain is obtained when the number of normal (hopping equal to 1) and strong (hopping equal to $`\sqrt{2}`$) bonds of the numerator coincides with the corresponding numbers of the denominator. Solving this simple combinatorial problem, one gets a probability

$$p(G=1)=\frac{1}{2^N}\underset{n=0}{\overset{N/2}{}}\left[\frac{(N/2)!}{(N/2n)!n!}\right]^2$$
(5)
for the peak of the probability distribution of the conductance at $`G=1`$curiosidad . Notice that the distribution is a sum of delta functions since the hopping takes just two different values; for example, the peak below the one at $`G=1`$ appears at $`G=8/9`$ and corresponds to $`\mathrm{\Sigma }_N/i=2`$ or $`\mathrm{\Sigma }_N/i=1/2`$ (one extra strong bond either in the numerator or the denominator). An asymptotic expansion for large enough chain lengths can be found for the sum in (5):

$$p(G=1)\sqrt{\frac{2}{\pi N}}.$$
(6)
The scaling behavior of this peak is enough to explain an inverse square root law for the mean value of $`G`$ given that other peaks of the distribution decay in a similar or faster way with the size of the chain. This is an interesting behavior since it perfectly coincides with scaling predictions for wide wires of odd number of transversal modesmudry .

Numerical simulations can be done as precise as necessary for this simple model. This fact allows a detailed comparison with theoretical predictions. Fig. 1 shows the scaling of the conductance as described by its mean value. Root mean square deviation is about some tenths of the conductance unit. In any case, the average conductance at large sizes ($`N>1000`$) shows a power law scaling with the distance that is compatible with the inverse square root law given by Eq.(6). According with this result, a disordered chain with pure non–diagonal disorder shows non–standard scaling at the band center. Nevertheless, alternative measures of the central value of the distribution restore to some way the usual exponential decay of one–dimensional conductances, even in the presence of chiral symmetry. For example, the geometric mean111 The geometric mean is obtained from the average of the logarithm of the random variable by:
$$<x>_{\mathrm{geom}}=\mathrm{exp}<\mathrm{ln}(x)>$$
of the same conductance distributions shows a much more pronounced decrease with length (see Fig. 2) although the corresponding standard deviation is of the same order as the mean making doubtful its statistical relevance. But there are other alternatives which do give a good description of the overall scaling of the distribution. Both the median222 The median of a probability distribution function $`p(x)`$ is the value $`x_{\mathrm{med}}`$ for which larger and smaller values of $`x`$ are equally probable:
$$_{\mathrm{}}^{x_{\mathrm{med}}}p(x)𝑑x=\frac{1}{2}=_{x_{\mathrm{med}}}^{\mathrm{}}p(x)𝑑x$$
or any definition based on the value of the integral of the distribution between 0 and an arbitrary upper limit $`G_{\mathrm{max}}`$ flow to exponentially small values. The physical meaning is clear in this case: half or more of the measures are exponentially small at large chain lengths. Actually, the precise scaling law for these central value alternatives is:

$$G_0\mathrm{exp}\left(\sqrt{L/\xi }\right),$$
(7)
where $`\xi `$ gives a measure of the exponential localization length. Fig. 2 shows that fits according to this law are excellent over the whole length range. The ultimate reason for such statistically disappointing results is simply the unusual size scaling of the distribution. While the major part of the distribution below $`G=1`$ decays in an algebraic form, the weight of the distribution accumulates in an exponentially small region near $`G=0`$ (Fig. 4 in Section V illustrates graphically this behavior). In this way, the upper part of the distribution dominates the scaling behavior of mean, root mean square deviation, etc., while accumulation at $`G=0`$ gives the behavior of central value definitions based on the integral of the distribution. In my opinion, these last definitions are better suited for characterizing the whole distribution than standard averages. Ultimately, one should look at the precise experimental protocol followed to get a value of the conductance before making predictions about the result of the measurements.

## III Non–diagonal disorder model

The lattice Hamiltonian describing random hopping on a $`L\times M`$ cluster of the square lattice is :

$$\widehat{H}=\underset{<ll^{}>}{}t_{ll^{}}\widehat{c}_l^{}\widehat{c}_l^{},$$
(8)
where $`\widehat{c}_l^{}`$ creates an electron on site $`l`$, $`l`$ and $`l^{}`$ are nearest–neighbor sites, and $`t_{ll^{}}`$ is the hopping energy from site $`l`$ to $`l^{}`$. It takes values 1 and -1 with equal probability. Let me refer to this model as the Random Hopping Sign (RHS) model. Obviously, the square lattice can be divided into two sublattices such that atoms belonging to one of them hop only to sites belonging to the other sublattice when described by Hamiltonian (8). Therefore, this Hamiltonian changes sign under a transformation that changes the sign of the electron operators on one sublattice. Consequently, the spectrum is symmetric relative to the band center at $`E=0`$ for any disorder realization, i.e., for any values of the random variables of the model $`\{t_{ll^{}}\}`$.

The model given by Eq.(8) is probably the simplest two dimensional model showing chiral symmetry. Many changes can be done to this model preserving chirality. For example, absolute values of the hopping can fluctuate in addition to their signs or complex values of the hopping parameters can be considered if random magnetic fluxes are simulated. The results described below are not sensible to any of these modifications of the disorder model. Numerical values are somewhat changed but trends remain exactly the same.

## IV Conductance calculation

Randomly generated samples of $`L\times M`$ clusters are connected to ideal leads of width $`M`$. Typically, $`LM`$ in a wire geometry. The conductance of the whole system is obtained using Kubo formulakubo  within exact one-electron linear response theory. Computational details are given elsewherejav2 ; cpc . Let me just mention that the inversion of the Hamiltonian matrix needed to get the Green function of the wire cannot proceed slab by slab as it is usually done within an optimized code. Numerical divergences take place owing to the existence of a true eigenstate at exactly $`E=0`$ for any piece of the system showing an odd number of sites (the matrix $`E\widehat{H}`$ is locally singular). Nevertheless, numerical calculation proceeds straightforward when pivoting over the whole Hamiltonian matrix that includes the ideal leads is allowed.

## V Results

### V.1 Wide disordered wires

The first aim has been the recovery of some important results obtained for wide wires by Mudry, Brouwer, and Furusakimudry . In particular, the exotic dependence of the scaling law on the parity of the wire width is obtained for the present model (a simplification of their Random Flux Model). Open boundary conditions have been used to get conductances of stripes of fixed width and number of open channels (the wire width is equal to the number of channels at the band center). Sample lengths have been varied from 99 to 1980 in steps of 99. As many as $`10^4`$ samples are necessary to get good values of means and other central values of conductance distributions. Fig. 3 shows the scaling law for two typical odd widths (9 and 19 channels) in a log-log plot. Error bars are comparable to the symbols representing the conductance averages. A power law fit to the numerical data is compatible with a mean conductance scaling proportional to the inverse square root of the sample length $`L`$:

$$<G>\frac{1}{\sqrt{L}}.$$

This is precisely the scaling law obtained in Ref.(mudry, ) for quantum wires of an odd section. Notice that scaling proceeds smoothly without distinguishing odd and even wire lengths. As noted by these authors, the variance of the conductance distribution is as large as its mean, making the mean a poor characteristic of the whole distribution. Actually, the predicted relationship

$$\frac{<G^2>}{<G>}=\frac{2}{3}$$

is accurately reproduced by the numerical data (see filled symbols of Fig. 3 and the dashed lines that are simply $`2/3`$ of the previous fits). On the other hand, when an even number of channels is studied, exponential scaling of the average conductance together with an exponentially small typical deviation is obtained in complete agreement with previous workmudry .

A deeper understanding of scaling of the conductance mean value can be obtained from the analysis of the whole conductance distribution. Here, conductances for $`10^5`$ samples of width $`9`$ and lengths $`10^2`$, $`10^3`$, and $`10^4`$ have been compiled and the corresponding histograms plotted in Fig. 4. It is clear that the probability for large values of the conductance diminishes as the length of the wire increases. Actually, since the plot is semilogarithmic the roughly equal separation between filled and empty circles for one side and between empty circles and diamonds for the other side, implies a power law decrease of the probability. The weight that the probability distribution loses for large conductance values goes close to $`G=0`$. I have checked that the divergence of the probability at $`G=1`$ is proportional to $`1/\sqrt{1G}`$ while the divergence at the origin looks also algebraic but with an exponent starting close to $`\frac{1}{2}`$ for small disorder and decreasing towards $`1`$ as disorder increases.

Although the conductance corresponding to ideal (non disordered) wires equals their widths (9 and 19, in this case), Fig. 4 shows that conductance just reaches the value of 1 for disordered systems. It seems that just one channel is effective. Therefore, it is tempting to compare these conductance distributions with the one corresponding to one random channel within the orthogonal universality classcoe  arguing that before localization the effect of disorder is just randomizing transport coefficients. But the comparison is very bad since apart for small differences all theories show probabilities that continuously decrease from $`G=0`$ to $`G=1`$. In particular, the simplest result that applies to a random channel described by scattering matrices of the Circular Orthogonal Ensemble is:

$$p(G)=\frac{1}{2\sqrt{G}},$$
(9)
that is quite different from the one obtained by numerical simulation. There is a direct mathematical reason explaining this failure. Typically, several states contribute to the Green function calculated at an arbitrary energy within the spectrum of a disordered system. But this is not the case when the Green function of a chiral system is calculated at $`E=0`$ which is an eigenenergy of the isolated system with an odd number of sites. In this situation, both the Green function of the isolated finite system and the one corresponding to the extended system including the leads (and related to the previous one by a Dyson equation) are dominated by the pole at $`E=0`$, that is, are completely determined by the $`E=0`$ eigenfunction. In the next subsection, I will exploit this feature to analyze the conductance distribution as a consequence of a precise wavefunction statistics. Meanwhile, let us return to the theory of Ref. (mudry, ) to analyze the numerical results.

Mudry, Brouwer, and Furusaki give the following expression for the conductance distribution of wires with an odd number of channels in the localized regime (see Eq. (4.10) of the second paper of Ref. (mudry, )):

$$p(G)=\sqrt{\frac{C}{\pi }}\frac{\mathrm{exp}\{C[\mathrm{arcosh}(G^{1/2})]^2\}}{G\sqrt{1G}},$$
(10)
where

$$C=\frac{l}{4L}\frac{M^2+M2}{M1},$$
(11)
with $`l`$ the mean free path, $`L`$ the wire length, and $`M`$ the wire width (which coincides with the number of channels at the bandcenter). A second scaling parameter $`\eta `$ that characterizes the disorder on a microscopic scale does not appear because it vanishes for the random flux model (RFM) and the model studied here (see Section III) is just a special case of RFM modeljav2 . When $`C`$ is small enough ($`C<0.01`$) the distribution is very well described by a much simpler expression (error smaller than one percent for $`G0`$):

$$p(G)\sqrt{\frac{C}{\pi }}\frac{1}{G\sqrt{1G}}$$
(12)
which allows the evaluation of the mean and the variance of the distribution:

$$<G>2\sqrt{\frac{C}{\pi }}$$

$$<G^2>=\frac{2}{3}<G>.$$

Using the explicit form of $`C`$ (Eq.(11)) an expression is obtained that can be used to fit the remaining parameter $`l`$ once an enough number of widths and/or lengths have been studied:

$$<G>^2L\frac{l}{\pi }\frac{M^2+M2}{M1}.$$
(13)
Fig. (5) shows a fit to the mean and variance values obtained for sets of $`10^4`$ randomly generated samples of length $`L=2000`$ and width $`M`$ from three to 37. While Eq.(13) works reasonably good for large widths it fails near $`M=3`$. Actually, an alternative fit by a linear law works sensibly better:

$$<G>^2L=aM+b.$$
(14)
Consequently, the distribution given by Eq. (10) will be used with

$$C=\frac{0.3515M+2.3050}{L}$$
(15)
to describe the results obtained by numeric simulation.

Probability $`p(G)`$ is shown in Fig. 4 as continuous lines. Although the three numerical distributions are nicely reproduced, accuracy is better for wider wires as expected from the theory ($`ML`$ is assumed in the theory leading to Eq. (10)). The divergence at $`G=1`$ is of the inverse square root form while the apparent $`G^1`$ non–integrable divergence at the origin is regularized by the complex numerator which vanishes at $`G=0`$. In conclusion, Eq.(10) is a very good description of the numerical data once the constant $`C`$ is properly estimated.

Let us try an alternative way of characterizing the central value of distributions of the form shown in Fig. 4. Previous experience with disorder models of this kind proves that the median is more robust than the mean for some distributions showing large $`G`$ variances but smaller $`\mathrm{ln}(G)`$ variances, i.e., when geometric mean does it better that the usual arithmetic meanjav2 . Unfortunately, the distribution of $`\mathrm{ln}(G)`$ is also very broad in the present case. Scaling results for the toy model of Section II suggest the use of the median. Fig. 6 gives the median scaling obtained for the previously collected conductance distributions. In can be seen that median values are exponentially smaller than mean values at large wire lengths. The fact that the median scales towards 0 can be inferred from the scaling of the histograms given in Fig. 4 although, the accumulation near exponentially small values of $`G`$ is not visible in the Figure. Even the scaling law of medians is the same that describes the chain model (see, Eq. (7)). Alternatively, the exponential decay of the median as a function of the square root of the wire length can be inferred from Eq. (10) once the distribution is written as a logarithmic normal distribution using the variable change described in the next subsection.

The practical meaning of the result is quite clear, more than half of the conductance observations are exponentially decreasing as the wire length increases. Actually, the fraction of the samples showing exponentially small conductances increases with wire length because probability below $`G=1`$ decreases monotonically. Just the difference with a conventional one–dimensional scaling (as, for example, the one obtained for wires of even width) comes from the power law decrease of the upper part of the probability distribution that should be compared with the exponential decrease characteristic of standard scaling.

### V.2 Two dimensional system

While a standard study of the scaling of the conductance in 2D would imply the calculation of conductances of increasing $`L\times L`$ samples connected to ideal leads of width equal to the square side $`L`$, in this work I have used a dot setup specially designed for the study of just one conducting statenote1 . Certainly, when the number of incoming channels is fixed by a point contact geometry, the only factor that affects the value of the conductance is the size (area) of the dot. In this way, the presumably increase of the conductance due to wider contacts does not obscure the underlying scaling law strictly due to the increased size of the system. Two different limits are well known for large values of $`L`$. First, ballistic transport through the sample can occur as it happens in chaotic cavities. This limit is described in a first approximation by scattering matrices of the Circular Orthogonal Ensemble (COE)coe . Roughly speaking, a conductance about $`\frac{1}{2}`$ per channel can be expected. Second, Anderson exponential localization would imply an exponentially small value of conductances for large enough dot sizes. This limit applies to diagonal disorder, for example. At the band center of a system with chiral symmetry, numerical simulation shows a behavior similar to the ballistic one (see Fig. 7). Mean conductance converges to a well–defined finite limit while the whole conductance distribution is perfectly described by Eq. (10). Since now the system is not quasi–onedimensional as wide wires are, one is forced to conclude that there should be deep general reasons for the validity of Eq. (10) in this context.

Let me briefly describe the numerical procedure. The transmission between any couple of points within the dot is calculated and the corresponding conductance distribution obtained. To this end, two clean infinite chains (the leads) are attached through two arbitrarily chosen lattice sites within the disordered square sample (the dot). For this geometry, the transmission from site $`𝐫`$ to site $`𝐫^{}`$ is given by the following expression:

$$T=\frac{4\psi (𝐫)^2\psi (𝐫^{})^2}{\left(\psi (𝐫)^2+\psi (𝐫^{})^2\right)^2},$$
(16)
where $`\psi (𝐫)`$ is the wavefunction at the band center. From this equation, $`0T1`$, and the conductance is between zero and one quantum unit for this numerical simulation. While the large number of transmission evaluations ($`\frac{1}{2}L^4`$) allows for a very precise calculation of the whole conductance distribution, the dependence of the conductance on the separation between point contacts is not given by the procedure. The eigenstate at the band center is obtained by direct inversion of the Schrödinger equation ($`\widehat{H}\psi =0`$) where $`\widehat{H}`$ is given by (8). Once the $`E=0`$ wavefunction statistics is known, Eq. (16) can be used to get a conductance distribution. For example, if the Porter-Thomas formpoto  were valid:

$$f(t)=\frac{1}{\sqrt{2\pi t}}\mathrm{exp}(t/2),$$
(17)
where $`t=N\psi (𝐫)^2`$ being $`N`$ the number of sites, the conductance distribution would be given by:

$$p(G)=_0^{\mathrm{}}_0^{\mathrm{}}𝑑t𝑑t^{}f(t)f(t^{})\delta \left(G\frac{4tt^{}}{(t+t^{})^2}\right).$$
(18)
which can be integrated to give the final result:

$$p_{\mathrm{PT}}(G)=\frac{1}{\pi }\frac{1}{\sqrt{G(1G)}}.$$
(19)
Although probability $`p_{\mathrm{PT}}(G)`$ reproduces some features of the conductance distributions shown in Figs. 4 and 7 for the smaller cluster sizes (for example, the square root divergence at $`G=1`$), it is clearly non comparable to the accurate result given by Eq.(10). What comes as a surprise is the fact that a logarithmic normal distribution of the wavefunction ratios $`t/t^{}`$ exactly gives the conductance distribution proposed by Mudry et al.mudry . That is, assuming

$$g(x)=\sqrt{\frac{C}{4\pi }}\mathrm{exp}(\frac{C}{4}x^2)$$
(20)
being $`x=\mathrm{ln}(t/t^{})`$, the conductance distribution is given by

$$_{\mathrm{}}^{\mathrm{}}𝑑xg(x)\delta \left(G\frac{4}{e^x+e^x+2}\right).$$
(21)
This integral is easily solved giving $`p(G)`$ of Eq.(10) with

$$G=\frac{4}{e^x+e^x+2}$$

or equivalently,

$$\frac{x}{2}=\mathrm{arcosh}(G^{1/2}).$$

This result gives some clue over the complicate $`G`$ dependence that happens to appear in its distribution function. Numerical simulation (see Fig. 8) shows that Eq.(20) accurately describes the wavefunction squared ratios of large disordered two-dimensional systems and, consequently, conductance distributions of the form given by Eq.(10) are valid in this case.

Let us now discuss the scaling properties of the mean conductance at the band center. Conductance has been averaged over a significant number of randomly generated samples of increasing linear sizes ($`L`$) keeping almost constant the total number of measurements. Fig. 9 shows the results obtained by this numerical procedure. The scaling of the mean conductance can be fitted to a model of the form:

$$<G>=a+bL^{\frac{1}{4}}$$
(22)
with a very good precision. The asymptotic value corresponding to the infinite limit is 0.245. Finite values of the conductance of two–dimensional systems showing chiral symmetry have been predicted in a number of papersgade ; furusaki ; mudry0 ; campos ; mudry2 . Although, the present numerical simulation nicely supports these theories some caution must be used for two main reasons. First, a somewhat practical reason. Results heavily depend on the fact that the Green function of this problem is just given by only one particular state. While the energy of this state is well defined theoretically, it could be difficult to make an experiment at just a particular energy. Previous authors on the subject have clearly shown that chirality is lost as soon as $`E=0`$ is leftmudry . Second reason is a bit more technical. Scaling properties have been obtained for clusters of an odd number of sites and, therefore, a state at $`E=0`$. Present computational facilities do not allow to prove that scaling of clusters of an even number of states proceeds in the same way. (Note that in this case the Green function should be recalculated for any new position of the point contacts since storing of the whole Green function matrix of the isolated cluster is not possible). Nevertheless, I have checked that odd/even differences are minimal for small systems.

The universal conductance distribution at the band center of a two–dimensional chiral system shown by the thick line in Fig.(7) is probably the most original result in the paper. It corresponds to the $`L\mathrm{}`$ limit of Eq.(22) which is obtained for $`C=0.0471`$ in Eq.(10) but does not differ very much from the distributions corresponding to the simulated sizes ($`L=599,699,799`$) plotted as three thin indistinguishable lines. All three are accurately described by Eq.(10) with $`C=0.0948`$. Although obtained for point contacts, previous experience with wide wires proves that minor differences should be expected for wider but finite contactsnote2 . Let me insist in the idea that perfect agreement with the theoretical conductance distribution is also explained by a particular distribution of wavefunction weight ratios (see Eq.(20)) that is confirmed by numerical simulation. Fig. 8 give the distributions obtained for a small cluster size ($`L=99`$) -higher maximum value- and for a larger size ($`L=499`$) together with its corresponding Gaussian (dashed line). The thick line gives the universal prediction for this statistical measure of the wavefunction. Consequently, I think that unusual phenomena as they have been described in the literature are quite compatible with present results for the conductance distributionfurusaki ; mudry2 ; campos .

## VI Discussion

The scaling properties of the conductance distribution of wide wires showing chiral symmetry has been analyzed by numerical solution of a simple model Hamiltonian. Predictions of a power law decrease of the mean conductance have been confirmed. Besides, the strange scaling properties of the distribution give rise to unusual statistical properties that have been described in the main text. For example, while the slow algebraic decay of the mean conductance is due to few but very large fluctuations, it has been shown that almost all conductance measurements are exponentially small for large wires. Actually, an exponential law is obtained for the median of the distribution and for any generalization of the median based on the integral of the distribution. Nevertheless, exponential decay is not proportional to the sample length but to its square root. That is a significant deviation from standard localization theory.

Surprisingly, statistical conductance properties of two–dimensional systems connected to point contacts are described at the band center by the same distribution that describes wide wires. It has been shown that this statistics appears due to the existence of an underlying Gaussian distribution of the logarithm of the ratio of the wavefunction weights at $`E=0`$. The main difference with quasi one–dimensional systems is that now both the mean and the median conductance seem to reach a finite limit. The universal statistical distributions of conductance and wavefunction ratios have been obtained. Although they share some properties with the corresponding statistics of the Circular Orthogonal Ensemble, it has been shown that they differ. This numerical result calls for some theory able to explain the origin of such a simple statistics of the wavefunction at the band center of a chiral system.

###### Acknowledgements.

This work has been partially supported by Spanish Comisión Interministerial de Ciencia y Tecnología (grant PB96-0085) and Dirección General de Enseñanza Superior e Investigación y Ciencia (grant 1FD97-1358).


text:
# Accelerated Particles from Shocks Formed in Merging Clusters of Galaxies

## 1 Introduction

Rich clusters contain thousands of galaxies and are the largest gravitationally bound systems in nature. Masses for rich clusters are $`10^{15}M_{}`$, with $`5`$–10% of the mass found in a hot intergalactic gas at temperatures of 2–12$`\mathrm{keV}`$. Rich clusters emit thermal bremsstrahlung with luminosities $`L_x10^{44}10^{45}`$$`\mathrm{ergs}\mathrm{s}^1`$(Allen & Fabian, 1994). Poor clusters contain hundreds of galaxies, have total masses $`10^{14}M_{}`$, have a hot intergalactic gas of temperatures 1–5$`\mathrm{keV}`$ and X-ray luminosities $`L_x10^{42}10^{43}`$$`\mathrm{ergs}\mathrm{s}^1`$(Dahlem & Thiering, 2000). Approximately 90% of the total mass of clusters is in the form of nonluminous dark matter.

In the hierarchical merging cluster scenario, poor clusters merge together to form richer clusters. Approximately 30–40% of galaxy clusters show evidence of substructure in both the optical (Beers, Geller & Huchra, 1982) and X-ray wavelengths (Forman et al., 1981). Velocity differences between the observed structures is $`1000`$$`2000`$$`\mathrm{km}\mathrm{s}^1`$. With gravitational forces driving the interaction between the two systems, cluster mergers are consistent with highly-parabolic orbits. Typical sound speeds within the intergalactic medium (IGM) are $`800`$$`\mathrm{km}\mathrm{s}^1`$, so shocks will form at the interaction boundary of the two systems. Computer models of merging clusters support the development of shocks in the IGM (Roettiger, Burns, & Loken, 1993; Ricker, 1998; Takizawa, 1999). Dimensional arguments show that a cluster merger releases $`10^{63}10^{64}`$$`\mathrm{ergs}`$ of gravitational potential energy when initial separations are of the order $``$ Mpc.

Only the most massive and X-ray luminous galaxy clusters have extended diffuse radio sources. With projected linear sizes $`1`$$`\mathrm{Mpc}`$, these diffuse sources have no known optical counterparts. The diffuse radio emissions have two distinct characteristics. The extended diffuse emission found in the central region of a galaxy cluster with a regular, azimuthally symmetric shape is known as a radio halo. The diffuse emission found on the cluster periphery are the cluster radio relics. These features often have irregular shapes with signs of filamentary structure. Radio relics are associated only with clusters that show evidence of a recent or ongoing merger event.

Shock fronts that form in the IGM as a result of a cluster merger event are thought to be associated with the cluster radio relics. The shock compression will orient any existing cluster magnetic field into the plane of the shock. Radio relics are characterized by highly organized magnetic fields with field strengths in the $`1`$$`\mu `$G range with linearly polarized field lines in the vicinity of the shock (Enßlin, 1998). The shock front will accelerate a fraction of the thermal particles within the IGM by first-order Fermi acceleration.

Recent work has highlighted the importance of nonthermal radiation from particles accelerated by shocks formed in merging clusters. Loeb and Waxman (Loeb & Waxman, 2000) have proposed that cluster mergers are the dominant contributor to the diffuse $`\gamma `$-ray background, provided the efficiency to convert the available gravitational energy into nonthermal electron energy is $`5`$%. Some unidentified EGRET sources are claimed to be associated with $`\gamma `$-ray emission from galaxy clusters (Colafrancesco & Blasi, 1998). Excess EUV emission from Coma can be explained by nonthermal electrons accelerated at merger shocks (Atoyan & Völk, 2000). Variations in radio surface brightness will result from superposition of cluster emissions (Waxman & Loeb, 2000).

In order to examine this question in more detail, we have modified a supernova remnant code (Sturner, et al., 1997) to calculate nonthermal radiation spectra from primary and secondary particles in merger shocks. We go beyond previous treatments by considering both primary electron and proton acceleration, a time-dependent treatment of radiation losses, and radiation signatures of secondaries from proton-nuclear interactions.

## 2 Model

We have adapted a supernova remnant code (Sturner, et al., 1997) to treat the cluster merger scenario. The code is designed to calculate time-dependent particle distribution functions evolving through adiabatic and radiative losses for electrons and protons accelerated by the first-order Fermi process at the cluster merger shock.

The electron and proton distribution functions originate from a momentum power-law injection spectrum. In terms of total particle energy $`E=m\gamma c^2`$, the injection function is

$$Q_{e,p}(E,t)=Q_{e,p}^0\left[\frac{(pc)^s}{\beta }\right]\mathrm{exp}\left[\frac{E}{E_{max}(t)}\right],$$
(1)
where $`p=\beta \gamma `$ is the dimensionless momentum and $`s`$ is the injection index. The maximum particle energy $`E_{max}`$ is determined by the maximum energy associated with the available time since the beginning of the merger event, by a comparison of the Larmor radius with the size scale of the system, and by a comparison of the energy-gain rate through first-order Fermi acceleration with the energy-loss rate due to adiabatic, synchrotron and Compton processes. Particle injection ceases after the age of the shock front exceeds $`t_{acc}=10^9`$ yrs. The constant $`Q_{e,p}^0`$ normalizes the injected particle spectrum over the entire volume $`V(t)`$ swept out by the shock front, and is determined by

$$E_{e,p}^{tot}=_0^{t_{acc}}𝑑t_0^{E_{max}}𝑑EEQ_{e,p}(E,t)V(t).$$
(2)
We assume a total available energy $`E_{e,p}^{tot}=\eta _{e,p}10^{63}`$$`\mathrm{ergs}`$, and an efficiency factor $`\eta _{e,p}=5\%`$ for both protons and electrons. Although the injection index $`s`$ depends upon the Mach number of the shock, here we present calculations for a fixed index $`s=2`$.

The time evolving particle spectrum is determined by solving the Fokker-Planck equation in energy space for a spatially homogeneous IGM, given by

$$\frac{n(E,t)}{t}=\frac{}{E}[\dot{E}_{tot}(E,t)n(E,t)]+\frac{1}{2}\frac{^2}{E^2}[D(E,t)n(E,t)]+Q(E,t)\frac{n(E,t)}{\tau _{pion}(E,t)}.$$
(3)
The quantity $`\dot{E}_{tot}(E,t)`$ represents the total synchrotron, Compton, Coulomb, and adiabatic energy-loss rate for electrons, and the sum of the Coulomb and adiabatic energy-loss rates for protons. Both protons and electrons are subject to diffusion in energy space by Coulomb interactions. The protons experience catastrophic losses due to proton-proton collisions on the time scale $`\tau _{pion}`$. The spectra of secondary electrons and positrons are calculated from pion-decay products, and are subject to the same physical processes as the primary electrons.

The synchrotron, Compton, bremsstrahlung, and pion-decay $`\gamma `$-ray spectral components are calculated from the particle spectra following the methods described by Sturner, et al. (1997). We use a standard parameter set with a mean IGM number density $`n_{IGM}=10^3\mathrm{cm}^3`$, a uniform cluster magnetic field of $`B=0.1\mu \text{G}`$, a constant shock speed $`v_s=1000`$$`\mathrm{km}\mathrm{s}^1`$, and an acceleration period of 1$`\mathrm{Gyr}`$. Thus, $`V(t)=A(t)v_st`$. We assume a constant surface area with a 1 Mpc radius. Particles loose energy due to adiabatic expansion according to the relation $`\dot{\gamma }/\gamma =\dot{V}(t)/V(t)=1/t`$.

## 3 Results and Discussion

Figure 1 shows nonthermal photon spectra calculated at 1$`\mathrm{Gyr}`$ and 5$`\mathrm{Gyr}`$ for our standard parameters. The system is very luminous at radio frequencies during the particle acceleration phase because of intense Compton losses of nonthermal electrons on the CMB. After the acceleration period ends, the radiation from primary electrons is dominated by emission from secondaries. The $`\pi ^0`$ bump at 70 MeV is hidden beneath the Compton-scattered CMB radiation from primary electrons during the acceleration period, but dominates the $`\gamma `$-ray spectrum after the acceleration period is over. The intensity of the $`\pi ^0`$ bump depends sensitively upon the relative efficiencies $`\eta _{p,e}`$ for proton and electron acceleration. If no pion signature is found in nonthermal $`\gamma `$ rays from merging clusters, then $`\eta _p\stackrel{<}{}\eta _e`$ for our standard parameters.

From the total photon spectra in Figure 1, we have calculated light curves at various photon energies in Figure 2 for a merging cluster at a distance of 100$`\mathrm{Mpc}`$. The maximum peak 1.4 GHz radio flux density of $``$ 12$`\mathrm{Jy}`$ occurs at the end of the acceleration period. The radio emission from the merger event is, however, distributed over an angular region $`0.5^{}`$ of the merger shock, making it difficult to detect.

Figure 2 presents calculations of the energy fluxes from a merger shock at 100 Mpc, for photon energies $`>100`$$`\mathrm{MeV}`$, $`>1`$$`\mathrm{GeV}`$, and $`>300`$$`\mathrm{GeV}`$. The maximum luminosity occurs at the end of the acceleration period due to the accumulation of nonthermal protons. Cooling times due to synchrotron and Compton energy losses are very short for high-energy electrons; electrons that Compton-scatter microwave photons to energies $`E=100E_{100}`$ MeV will cool on a time scale of $`t_{cool}6\times 10^6/\sqrt{E_{100}}`$ yrs. Consequently, the $`\gamma `$-ray emission declines sharply after particle acceleration ceases, and approaches a level where all the $`\gamma `$-ray production originates from proton interactions. The maximum flux for secondary electrons is $`2`$ orders of magnitude less than the maximum flux for the primary electrons, as can be seen by comparing the catastrophic proton loss time scale of $`30`$$`\mathrm{Gyr}`$ with the 1 Gyr injection time scale. Because of adiabatic losses and the long catastrophic pion production cooling times, protons inject a slowly declining rate of secondary electrons in the 1-4 Gyrs after the initial acceleration event.

## 4 Summary

We have studied the acceleration of protons and electrons in shocks formed by merging clusters of galaxies, and calculated the expected synchrotron, Compton, bremsstrahlung, and pion emission from the accelerated particles. Both the radio flux density at 1.4$`\mathrm{GHz}`$ and the $`\gamma `$-ray flux at photon energies $`>100`$$`\mathrm{MeV}`$, $`>1`$$`\mathrm{GeV}`$ and $`>300`$$`\mathrm{GeV}`$ for a source at 100$`\mathrm{Mpc}`$ as a function of time are presented. The flux is greatest at the end of the particle acceleration phase, and merging clusters of galaxies at a characteristic distance of 100$`\mathrm{Mpc}`$ should be detectable with radio telescopes and GLAST. Radio and $`\gamma `$-ray detectability of merger events must contend with the angular extent of the emission.

Because the merger event lasts for only $`1`$ Gyr, a large fraction of clusters will no longer be experiencing ongoing nonthermal particle acceleration from cluster mergers. GLAST, with a limiting sensitivity of $`4\times 10^{13}`$ ergs cm-2 s-1 for a one year all-sky survey, will still be able to detect such systems if they are at optimal distances. Depending on the shock speed and magnetic field, cluster merger events can also be detectable sources of TeV radiation (see Fig. 2). We therefore expect detectable $`\gamma `$-ray emission from clusters at distances $`100`$$`\mathrm{Mpc}`$. Within this distance a number of radio relics have been detected (Giovannini & Feretti, 2000; Kempner & Sarazin, 2001).


text:
# Interacting Stellar Wind and Photoionization Models of the SN 1987A remnant

## 1. Introduction

Observations of the SN 1987A nebula give us the clues that we will need to model its structure. Figure 1 shows an HST image of the nebula (Burrows et al. 1995) and a schematic diagram of the nebular structure. The HST image shows two of the salient features of the nebula: the inner and outer rings. The diagram shows a slice containing the symmetry axis. The inner ring lies in the equatorial plane at a radius of about $`0.6`$ lt yr and is expanding with a velocity of 10.3 km s<sup>-1</sup> (Crotts & Heathcote 1991), while the outer rings are at an angle of about 55 from the equatorial plane and an axial radius of about $`1.4`$ lt yr. A tenuous shell connects the inner and outer rings (Crotts, Kunkel, & Heathcote. 1995), but the shell does not extend past the outer rings; that is, the nebula is open along the symmetry axis.

In an interacting winds model for the nebula, a fast, hot wind ejected during a blue supergiant phase overtakes a slow, dense wind ejected from a previous red supergiant phase. In order to produce the equatorial ring the RSG wind must be asymmetrical, with higher density at the equator and lower density at the symmetry axis. Blondin & Lundqvist (1993) performed numerical simulations of this interacting winds scenario for adiabatic and radiatively cooling gas. These simulations showed that the interacting winds model could produce the inner ring and the overall shape of the nebula. However, the Blondin and Lundqvist models did not produce the outer rings, and the shell in these models is complete all the way to the symmetry axis. Therefore, the interacting wind models need some refinement if they are to describe the SN 1987A nebula.

Radio and x-ray observations of the nebula since 1990 provide additional information. These observations indicate that the supernova blast has encountered a region of density intermediate to the inner and outer winds (Gorenstein, Hughes, & Tucker 1994; Beuermann, Brandt, & Pietsch 1994; Staveley-Smith et al. 1993; Gaensler et al. 1997). The slow expansion of the radio source indicates that this material is too dense to be shocked BSG wind. Chevalier & Dwarkadas (1995) proposed that the intermediate density material is RSG wind that was photoionized by radiation from the progenitor, a B3 Ia star.

In this paper we explore the interacting winds scenario, including photoionization from the progenitor star, as a model for the SN 1987A nebula. Section 2. discusses the features of our numerical fluid dynamics code. Section 3. gives some preliminary results from our numerical simulations, and section 4. compares the results to SN 1987A observations and discusses our plans for future refinements to the models.

## 2. Numerical methods

For our calculations we use a two-dimensional Eulerian finite difference code with a van Leer monotonic transport algorithm. We use an ideal gas equation of state, $`p=(\gamma 1)ϵ`$, where $`ϵ`$ is the internal energy density, and the adiabatic constant, $`\gamma =5/3`$. The coordinate system is spherical-polar with reflecting boundary conditions at the equator and pole. The inner boundary condition is an inflow condition, set to match the characteristics of the BSG wind, and similarly the outer boundary condition is an outflow condition set to the characteristics of the RSG wind.

In addition to the standard fluid dynamics equations, the code has several features that are necessary for this particular problem. These features include grid expansion, ionization and recombination, radiative cooling, and separate advection for neutral and ionized gas. They are described fully in Rosenberg, Link, & Chevalier (2000).

The RSG wind is characterized by its the mass loss rate, velocity, and temperature, $`\dot{M}_R`$, $`v_R`$, and $`T_R`$. Since the RSG wind is asymmetric, $`\dot{M}_R`$ is actually the mass loss rate that would be experienced if the rate were everywhere equal to its value at the equatorial plane. The angular dependence of the mass loss rate is described by

$$\dot{M}_R(\theta )=\dot{M}_R\left(\frac{\pi }{2}\right)\left[1A\frac{\mathrm{exp}(2\beta \mathrm{cos}^2(\theta ))1}{\mathrm{exp}(2\beta )1}\right],$$
(1)
where $`A=11/R`$, $`R`$ is the ratio of equatorial to polar density in the wind, and $`\beta `$ controls the steepness of the falloff from the equatorial value. Larger values of $`\beta `$ produce distributions more strongly peaked at the equator; while smaller values produce a more gradual falloff. The BSG wind is similarly described by its mass loss rate, velocity, and temperature, $`\dot{M}_B`$, $`v_B`$, and $`T_B`$. Finally, the ionizing flux is described by the total flux of photons above the ionization threshold, $`S_{}`$. The temperatures $`T_A`$ and $`T_B`$ are sufficiently low that they do not affect the results. Compared to the computations of Blondin & Lundqvist (1993), $`S_{}`$ is the only additional parameter.

Table 1 gives the parameters for the models used in our calculations. Model ‘A’ is a fiducial model that was chosen for comparison with previous results (*e.g.* Blondin & Lundqvist 1993).

## 3. Results

Figure 2 shows a greyscale image of the final state of the fiducial calculation. The criterion for stopping the calculation is that the radius of the inner ring in the calculation reaches $`6\times 10^{17}`$ cm, which is approximately the measured radius of the inner ring in the nebula. The figure shows clearly the inner ring and HII region. However, there is no clear outer ring. Moreover, although the ionization front has broken out at the middle latitudes, the neutral shell is still intact near the pole, in contrast to observations. Further evolution of the model would result in photoionization of the polar shell.

Model ‘XA’ (fig. 2) has identical parameters to model ‘A’, except that the photoionizing flux has been turned off. Since there is no photoionizing flux, there is no HII region. The properties of the inner ring, including the expansion velocity, are mostly unchanged from the photoionization case; however, the structure in the middle latitudes is quite different. Particularly striking is the fragmentation in the neutral shell, which leads to dense clumps that look like they might form outer rings; however, these clumps appear to be formed by a Rayleigh-Taylor or similar instability, and therefore they should not be expected to exhibit the azimuthal symmetry required to form the outer rings.

The left panel of figure 3 shows the final state of Model ‘B’. Model ‘B’ has a stronger BSG wind than the fiducial model, and it also has a more strongly peaked (higher $`\beta `$) RSG wind. We see in figure 3 the beginning of a high-latitude ionization front breakout that could eventually destroy the neutral shell near the pole. There is also a neutral clump that might form an outer ring; however, a mid-latitude ionization front breakout has disconnected the clump from the rest of the neutral shell. The HII region is thick throughout most of the model; however, right at the equatorial plane it becomes very thin. Apparently the ionized gas is advected away from the equator, possibly due to the strong peak in the density of the RSG wind.

Finally, Model ‘C’ is shown in the right panel of figure 3. The density of the RSG wind is lower than in the other simulations we have discussed; consequently, most of the RSG wind above about $`30^{}`$ is ionized. A dense blob of neutral gas remains at mid-high latitude. This neutral gas appears to survive longer than the rest of the shell because it is compressed by converging shocks preceding the ionization fronts. It is possible (although not certain) that such a blob could form a ring in three dimensions; however, it does not seem possible to form both a ring and a connecting shell through this mechanism. Moreover, the ‘ring’ is short-lived. The photoionizing flux must last long enough to blow away the surrounding shell and then turn off before it blows away the ring as well. Such a coincidence seems unlikely.

Table 2 summarizes the results of our calculations. We have chosen to focus on three diagnostic properties: the thickness of the HII region at the equator (expressed as a ratio of its inner and outer radii), the presence or absence of a density enhancement that (if actually azimuthally symmetric) could form an outer ring, and the expansion velocity of the inner ring.

## 4. Discussion

For the most part the interacting winds scenario seems to be on the right track for explaining the SN 1987A nebula. The expansion velocity of the inner ring is a bit low in these simulations, but it is close enough that it is reasonable to expect that some combination of parameters could be found to produce an acceptable match to the observed value.

The HII region is the most encouraging development in the simulations. Although the HII region is somewhat thin in the equatorial plane, the overall thickness seems to indicate that this problem could be fixed by judicious choice of wind parameters. Since the equatorial thickness of the HII region seems to go down with increasing $`\beta `$, further calculations will concentrate on calculations with smaller $`\beta `$.

Although we do see mid-latitude clumps in several of our simulations, we are skeptical that this means a mechanism has been found for producing the outer rings. We have seen that it is possible for these clumps to form even in the nonphotoionizing case, as the neutral shell breaks up due to Rayleigh-Taylor instabilities. Consequently, it seems likely that the clumps would not form axisymmetric rings in three dimensions. However, the edge of the nebular shell created by the ionization front could have an axisymmetric structure because its position is related to the density gradient in the RSG wind. Whether a dense ring could form there depends on the detailed cooling properties of the gas, which are only approximately accounted for in our models.

### Acknowledgments.

This work was supported in part by NASA grant NAG5-8232.

## References

Beuermann, K., Brandt, S., & Pietsch, W., A&A, 281, L45

Blondin, J. M., & Lundqvist, P. 1993, ApJ, 405, 337

Burrows, C. J., et al. 1995, ApJ, 452, 680

Chevalier, R. A. & Dwarkadas, V. V. 1995, ApJ, 452, L45

Crotts, A. P. S. & Heathcote, S. R. 1991, Nature, 350, 683

Crotts, A. P. S., Kunkel, W. E., & Heathcote, S. R. 1995, ApJ, 438, 724

Gaensler, B. M., Manchester, R. N., Staveley-Smith, L., Tzioumis, A. K., Reynolds, J. E., & Kesteven, M. J. 1997, ApJ, 479, 845

Gorenstein, P., Hughes, J. P., & Tucker, W. H. 1994, ApJ, 420, L25

Rosenberg, D. L., Link, R., & Chevalier, R. A. 2000, in preparation

Staveley-Smith, et al. 1993, Nature, 366, 136