text:
# Cosmic microwave background and large scale structure limits on the interaction between dark matter and baryons

## I Introduction

It is now generally accepted that the energy density of the universe includes a substantial fraction of cold dark matter (CDM), constituting roughly $`0.30.4`$ of the closure density. It is normally assumed that the CDM particle does not interact with ordinary baryonic matter or with photons; hence the “dark” in cold dark matter. However, a small scattering cross section cannot be entirely ruled out, only constrained by observations.

Several different models are possible for interactions between ordinary matter and CDM. One possibility is strongly interacting dark matter (also known as strongly interacting massive particles, or SIMPs), in which the CDM particle couples to baryons, but not to photons or electrons. Recently, to solve the small scale problem in CDM models, strongly self-interacting particles were suggested Spergel ; Wandelt , with

$$8\times 10^{25}<\frac{\sigma /\mathrm{cm}^2}{m/\mathrm{GeV}}<10^{23}.$$
(1)
The range quoted here (from Ref. Wandelt ) is narrower than originally proposed in Spergel  because of additional constraints. While the original Spergel-Steinhardt cross section applies only to dark matter self-interaction, if such an interaction is mediated by the strong force, then the dark matter might also interact strongly with baryonsWandelt . Examples of such dark matter includes, e.g., quark-gluino bound states Farrar , strangelets Jaffe , gauge singlet mesons Bento , and Q-balls Kusenko . Limits on such interactions were investigated by Starkman et al. Starkman , and recently reexamined in references Wandelt ; QW ; Moh ; Cyburt . A second possibility is to couple the dark matter particle electromagnetically, e.g., by giving it a tiny electromagnetic charge (see, e.g., reference millicharge  and references therein); such a CDM particle will scatter primarily off of photons and free electrons.

Any interactions of this sort will affect the spectrum of the cosmic microwave background (CMB) fluctuations, since they act to transfer momentum from the dark matter to the baryon-photon fluid, although for WIMPs this effect is negligible Chen ; Hofmann . Current CMB observations can place constraints on such models. The case of direct interactions between the CDM particle and photons has already been examined Boehm , so we do not consider it here. While the model discussed in reference Boehm  is not identical to the case of electromagnetically charged dark matter, since it includes photon-CDM interactions but not electron-CDM interactions, we expect it to be qualitatively similar to the charged dark matter model. Therefore, we consider only the case of strongly-interacting dark matter, in which the dark matter particle couples to baryons (only).

In the next section, we give the perturbation equations for CDM which couples to baryons, and we show the matter transfer function and CMB power spectra for some representative cases. In Sec. 3, we compare our calculations with the CMB observations to derive constraints on the cross section for dark matter-baryon scattering, and we compare to previously-derived limits. Our conclusions are summarized in Sect. 4.

## II Fluctuation growth with strongly-interacting dark matter

### II.1 The Basic Modified Equations

We begin by assuming a scattering cross section $`\sigma `$ between the CDM particle of mass $`m_c`$ and protons with mass $`m_b`$, and we derive the changes to the perturbation equations produced by this scattering interaction. Our starting point is the set of perturbation equations in reference MB95 , which form the basis for CMBFAST CMBFAST . We will not reproduce here all of the perturbation equations, but simply give those equations which are modified.

Following reference MB95 , we work in the synchronous gauge. However, because the CDM is interacting it does not provide a natural way of defining the synchronous coordinates, as opposed to standard CDM MB95 . Let $`\delta _c`$ and $`\delta _b`$ be the density fluctuation ($`\delta \rho /\rho `$) for the cold dark matter and the baryons, respectively, and let $`\theta _c`$ and $`\theta _b`$ be the corresponding velocity divergences. Then for the cold dark matter, we have

$`\dot{\delta }_c`$ $`=`$ $`\theta _c{\displaystyle \frac{1}{2}}\dot{h},`$ (2)
$`\dot{\theta }_c`$ $`=`$ $`{\displaystyle \frac{\dot{a}}{a}}\theta _c+c_c^2k^2\delta _b+R_c`$ (3)
For baryons the corresponding terms are

$`\dot{\delta }_b`$ $`=`$ $`\theta _b{\displaystyle \frac{1}{2}}\dot{h},`$ (4)
$`\dot{\theta }_b`$ $`=`$ $`{\displaystyle \frac{\dot{a}}{a}}\theta _b+c_b^2k^2\delta _b+R_\gamma +R_b.`$ (5)
with

$`R_b`$ $``$ $`K_c(\theta _b\theta _c)={\displaystyle \frac{m_cn_c\sigma v_{\mathrm{rel}}a}{m_b+m_c}}(\theta _b\theta _c),`$ (6)
$`R_c`$ $``$ $`K_b(\theta _c\theta _b)={\displaystyle \frac{m_bn_b\sigma v_{\mathrm{rel}}a}{m_b+m_c}}(\theta _c\theta _b),`$ (7)
$`R_\gamma `$ $``$ $`K_\gamma (\theta _\gamma \theta _b)={\displaystyle \frac{4\rho _\gamma }{3\rho _b}}an_e\sigma _T(\theta _\gamma \theta _b)`$ (8)
In these equations $`c_c`$ and $`c_b`$ are the adiabiatic sound speed of the CDM and baryons, respectively, $`v_{\mathrm{rel}}`$ is the average relative velocity between the CDM and the baryons, and $`R_c,R_b,R_\gamma `$ are the momentum transfer terms. The sound speeds are given by

$`c_b^2`$ $`=`$ $`(k_BT_b/\mu )\left(1{\displaystyle \frac{1}{3}}{\displaystyle \frac{d\mathrm{ln}T_b}{d\mathrm{ln}a}}\right),`$ (9)
$`c_c^2`$ $`=`$ $`(k_BT_c/m_c)\left(1{\displaystyle \frac{1}{3}}{\displaystyle \frac{d\mathrm{ln}T_c}{d\mathrm{ln}a}}\right),`$ (10)
the average relative velocity by

$$v_{\mathrm{rel}}=\sqrt{\frac{3(m_cT_b+m_bT_c)}{m_bm_c}},$$
(11)
In calculating average velocities we have assumed thermal Maxwell-Boltzmann distributions for both CDM and baryons. The temperature evolution of the dark matter and baryons is given by

$`\dot{T}_b+2{\displaystyle \frac{\dot{a}}{a}}T_b`$ $`=`$ $`K_\gamma (T_bT_\gamma )K_c(T_bT_c),`$ (12)
$`\dot{T}_c+2{\displaystyle \frac{\dot{a}}{a}}T_c`$ $`=`$ $`K_b(T_cT_b).`$ (13)
It is possible to obtain analytical solutions for the temperature evolution equations, but they are complicated and cumbersome to use. Instead, we use a semi-implicit scheme to integrate these differential equations. We tested this method in a few cases, and found that it produces solutions which agree with the analytical solutions within 0.2%, and the effect on the CMB spectrum is negligible.

### II.2 The Effect of Primordial Helium

This set of equations ignores the fact that the baryonic fluid consists of two species: hydrogen and helium-4, in a roughly 3 to 1 ratio by mass. To account for the presence of helium as well as hydrogen, we make the following replacements:

$`{\displaystyle \frac{m_bn_b\sigma v_{\mathrm{rel}}}{m_b+m_c}}`$ $``$ $`{\displaystyle \frac{m_\mathrm{H}n_\mathrm{H}\sigma _\mathrm{H}v_\mathrm{H}}{m_\mathrm{H}+m_c}}+{\displaystyle \frac{m_{\mathrm{He}}n_{\mathrm{He}}\sigma _{\mathrm{He}}v_{\mathrm{He}}}{m_{\mathrm{He}}+m_c}},`$ (14)
$`{\displaystyle \frac{m_cn_c\sigma v_{\mathrm{rel}}}{m_b+m_c}}`$ $``$ $`{\displaystyle \frac{m_cn_c\sigma _\mathrm{H}v_\mathrm{H}}{m_\mathrm{H}+m_c}}(1Y)+{\displaystyle \frac{m_cn_c\sigma _{\mathrm{He}}v_{\mathrm{He}}}{m_{\mathrm{He}}+m_c}}Y,`$ (15)
where $`Y`$ is the primordial mass fraction of helium-4 (we take $`Y=0.24`$ throughout), and we define separate velocities for the hydrogen and helium relative to the dark matter:

$`v_\mathrm{H}`$ $`=`$ $`\sqrt{{\displaystyle \frac{3(m_cT_b+m_\mathrm{H}T_c)}{m_\mathrm{H}m_c}}}`$ (16)
$`v_{\mathrm{He}}`$ $`=`$ $`\sqrt{{\displaystyle \frac{3(m_cT_b+m_{\mathrm{He}}T_c)}{m_{\mathrm{He}}m_c}}}.`$ (17)
In our calculations, $`\sigma `$, the scattering cross section between the CDM particle and the proton, is a free parameter which we will constrain from the observations. However, we need to make some sort of assumption regarding $`\sigma _{He}`$, the scattering cross section between CDM and helium nuclei. For coherent scattering, we expect $`\sigma _{He}=16\sigma `$, while a spin-dependent cross section would yield $`\sigma _{He}=0`$. For now, we will examine the intermediate case of incoherent scattering, for which $`\sigma _{He}=4\sigma `$. Later, in Sec. 3, we will extend our calculation to consider both coherent scattering and spin-dependent scattering.

Taking $`m_{\mathrm{He}}4m_\mathrm{H}`$, and $`n_b=\rho _b/m_\mathrm{H}`$, the above substitution can be achieved by simply replacing everywhere

$`{\displaystyle \frac{m_b\sigma v_{\mathrm{rel}}}{m_b+m_c}}`$ $``$ $`{\displaystyle \frac{m_\mathrm{H}\sigma }{m_\mathrm{H}+m_c}}\left[(1Y)v_\mathrm{H}+Y\left({\displaystyle \frac{\sigma _{\mathrm{He}}}{\sigma _\mathrm{H}}}\right){\displaystyle \frac{m_\mathrm{H}+m_c}{4m_\mathrm{H}+m_c}}v_{\mathrm{He}}\right]`$ (18)
$`{\displaystyle \frac{m_c\sigma v_{\mathrm{rel}}}{m_b+m_c}}`$ $``$ $`{\displaystyle \frac{m_c\sigma }{m_\mathrm{H}+m_c}}\left[(1Y)v_\mathrm{H}+Y\left({\displaystyle \frac{\sigma _{\mathrm{He}}}{\sigma _\mathrm{H}}}\right){\displaystyle \frac{m_\mathrm{H}+m_c}{4m_\mathrm{H}+m_c}}v_{\mathrm{He}}\right].`$ (19)
### II.3 Tight-Coupling Approximation

Following reference MB95 , we use a separate set of equations in the limit of tight coupling between baryons and photons. In the tight-coupling regime, the temperature evolution is given by

$`\dot{T}_b`$ $`=`$ $`2{\displaystyle \frac{\dot{a}}{a}}T_b+{\displaystyle \frac{8}{3}}{\displaystyle \frac{\mu }{m_e}}{\displaystyle \frac{\rho _\gamma }{\rho _b}}an_e\sigma _T(T_\gamma T_b)+K_c(T_cT_b)`$ (20)
$`\dot{T}_c`$ $`=`$ $`2{\displaystyle \frac{\dot{a}}{a}}T_c+K_b(T_bT_c)`$ (21)
where $`\mu `$ is the mean molecular weight. During tight coupling, $`T_cT_bT_\gamma `$, $`\dot{T_b}\frac{\dot{a}}{a}T_b`$, and $`\dot{v}_{\mathrm{rel}}=\frac{1}{2}\frac{\dot{a}}{a}v_{\mathrm{rel}}`$. Setting $`R=4\rho _\gamma /3\rho _b`$, and $`\tau _\gamma =(an_e\sigma _T)^1`$, we have

$$\dot{R}_b+\frac{\dot{a}}{a}R_b=\frac{m_cn_c\sigma av_{\mathrm{rel}}}{m_b+m_c}\left[(\dot{\theta }_c\dot{\theta }_b)\frac{3}{2}\frac{\dot{a}}{a}(\theta _c\theta _b)\right]$$
(22)
$$(1+R)\dot{\theta _b}+\frac{\dot{a}}{a}\theta _bc_b^2k^2\delta _bk^2R\left(\frac{1}{4}\delta _\gamma \sigma _\gamma \right)+R(\dot{\theta }_\gamma \dot{\theta }_b)R_b=0,$$
(23)
$$\theta _b\theta _\gamma =\frac{\tau _\gamma }{1+R}\left[\frac{\dot{a}}{a}\theta _b+k^2(c_b^2\delta _b\frac{1}{4}\delta _\gamma +\sigma _\gamma )+\dot{\theta }_\gamma \dot{\theta }_b+R_b\right],$$
(24)
$$\dot{\theta }_b\dot{\theta }_\gamma =\frac{\tau _\gamma }{1+R}\left[\frac{\ddot{a}}{a}\theta _b\frac{\dot{a}}{a}k^2\frac{\delta _\gamma }{2}+k^2\left(c_b^2\dot{\delta }_b\frac{1}{4}\dot{\delta }_\gamma \right)+\frac{\dot{a}}{a}R_b+\dot{R}_b\right]+\frac{2R}{1+R}\frac{\dot{a}}{a}\left(\theta _b\theta _\gamma \right).$$
(25)
### II.4 Applicable range of our calculation

In the calculation described above, we have assumed that the dark matter is made of a single type of particle, with density parameter $`\mathrm{\Omega }_m`$. We have also assumed that it is a “cold dark matter” particle, i.e. during the whole range of our calculation it is non-relativistic, which translates to $`m>1`$ MeV. We have also assumed that during each scattering only one baryon interacts with the dark matter particle. If the interaction length of the particle is $`a`$ ($`\sigma 4\pi a^2`$), then we require that this interaction length be less than the inter-particle spacing:

$$a<n_b^{1/3},$$
(26)
where $`n_b10^5(1+z)^3h_0^2\mathrm{\Omega }_b\mathrm{cm}^3`$. Processes which affect CMB anisotropy happen at redshift less than $`10^5`$, so we can put a limit on the applicable cross section, neglecting factors of order unity:

$$\sigma _{\mathrm{max}}10^6\mathrm{cm}^2.$$
(27)
Our calculation must be modified for $`\sigma >\sigma _{\mathrm{max}}`$ to take into account multi-particle scattering.

### II.5 The matter and CMB power spectra

We have modified the CMBFAST code CMBFAST  to include the effects discussed in the previous sections. For initial conditions, we assume that the dark matter and baryons started tightly coupled with the same temperature. Tensor modes have been ignored, and we consider only flat geometries ($`\mathrm{\Omega }=1`$).

In Fig. 1 we show the CMB fluctuation spectra and the corresponding square of the matter transfer function for some representative cases. The coupling of the dark matter to the baryons damps the matter power spectrum on small scales, as can be seen in Fig. 1(b).

One might imagine that the effect on the CMB of coupling the baryons to the dark matter would be equivalent to a standard model with a larger value of $`\mathrm{\Omega }_b`$. This is not the case, as can be seen in Fig. 1(a). Although both models produce an increase in the amplitude of the first acoustic peak, the large $`\mathrm{\Omega }_b`$ model produces an increase in the amplitude of the third peak, while the strongly-interacting dark matter model yields a decrease in the third peak amplitude (and in all of the other peaks as well).

One can understand the difference between the CMB spectra in the following way: The baryon-photon oscillations are damped oscillations, with a damping rate (viscosity) inversely proportional to the baryon-photon interaction rate. Interacting CDM with intermediate cross section amounts to adding “baryons” but simultaneously also increasing the viscosity of the plasma. Therefore the acoustic oscillations damp quickly and all peaks but the first are suppressed. Scales around the first acoustic peak have only undergone part of an oscillation and therefore have not been damped significantly. On such large scales the increase in oscillation amplitude due to added “baryons” is more important and the overall power spectrum increases.

On the other hand, increasing the baryon density has the well-known effect of increasing the height of the compression peaks (1st, 3rd,…) and decreasing the height of the decompression peaks (2nd, 4th,…). Of course, for very large cross sections, the effective viscosity of the CDM-baryon-photon plasma approaches that of a pure baryon-photon plasma and the CDM acts exactly like baryons. An effect similar to the increased viscosity due to CDM-baryon interactions can be seen in models with increased width of the last scattering surface due to non-standard recombination, where the effective viscosity of the baryon-photon plasma also increases around the time of recombination HS01 .

The effect on the matter power spectrum seen in Fig. 1(b) can be understood in the same way. The increased viscosity in the interacting model leads to a much stronger damping on small scales than in the model with increased baryon density.

## III Comparison with observations

In order to constrain models with CDM-baryon interactions we compare the CMB and matter power spectra of such models with recent observational data.

CMB data set — Several data sets of high precision are now publicly available. In addition to the COBE Bennett:1996ce  data for small $`l`$ there are data from BOOMERANG boom , MAXIMA max , DASI dasi  and several other experiments WTZ ; qmask . Wang, Tegmark and Zaldarriaga WTZ  have compiled a combined data set from all these available data, including calibration errors. In the present work we use this compiled data set, which is both easy to use and includes all relevant present information.

LSS data set — At present, by far the largest survey available is the 2dF 2dF  of which about 147,000 galaxies have so far been analysed. Tegmark, Hamilton and Xu THX  have calculated a power spectrum, $`P(k)`$, from this data, which we use in the present work. The 2dF data extends to very small scales where there are large effects of non-linearity. Since we only calculate linear power spectra, we use (in accordance with standard procedure) only data on scales larger than $`k=0.2h\mathrm{Mpc}^1`$, where effects of non-linearity should be minimal.

The CMB fluctuations are usually described in terms of the power spectrum, which is again expressed in terms of $`C_l`$ coefficients as $`l(l+1)C_l`$, where

$$C_l|a_{lm}|^2.$$
(28)
The $`a_{lm}`$ coefficients are given in terms of the actual temperature fluctuations as

$$T(\theta ,\varphi )=\underset{lm}{}a_{lm}Y_{lm}(\theta ,\varphi ).$$
(29)
Given a set of experimental measurements, the likelihood function is

$$(\mathrm{\Theta })\mathrm{exp}\left(\frac{1}{2}x^{}[C(\mathrm{\Theta })^1]x\right),$$
(30)
where $`\mathrm{\Theta }=(\mathrm{\Omega },\mathrm{\Omega }_b,H_0,n,\tau ,\mathrm{})`$ is a vector describing the given point in parameter space, $`x`$ is a vector containing all the data points, and $`C(\mathrm{\Theta })`$ is the data covariance matrix. This applies when the errors are Gaussian. If we also assume that the errors are uncorrelated, it can be reduced to the simple expression, $`e^{\chi ^2/2}`$, where

$$\chi ^2=\underset{i=1}{\overset{N_{\mathrm{max}}}{}}\frac{(C_{l,\mathrm{obs}}C_{l,\mathrm{theory}})_i^2}{\sigma (C_l)_i^2},$$
(31)
is a $`\chi ^2`$-statistic and $`N_{\mathrm{max}}`$ is the number of power spectrum data points oh . In the present letter we use equation (31) for calculating $`\chi ^2`$. In the case where we also use LSS data $`\chi ^2`$ is instead given by

$$\chi ^2=\underset{i=1}{\overset{N_{\mathrm{max},\mathrm{CMB}}}{}}\frac{(C_{l,\mathrm{obs}}C_{l,\mathrm{theory}})_i^2}{\sigma (C_l)_i^2}+\underset{j=1}{\overset{N_{\mathrm{max},\mathrm{LSS}}}{}}\frac{(P(k)_{\mathrm{obs}}P(k)_{\mathrm{theory}})_j^2}{\sigma (P(k))_j^2}.$$
(32)
The procedure is then to calculate the likelihood function over the space of cosmological parameters. The 2D likelihood function for $`(m_c,\sigma )`$ is obtained by keeping $`(m_c,\sigma )`$ fixed and marginalising over all other parameters.

As free parameters in the likelihood analysis we use $`\mathrm{\Omega }_m`$, the matter density, $`\mathrm{\Omega }_b`$, the baryon density, $`H_0`$, the Hubble parameter, $`n`$, the scalar spectral index, $`\tau `$, the optical depth to reionization, and $`Q`$, the overall normalization of the data. When large scale structure constraints are included we also use $`b`$, the normalization of the matter power spectrum, as a free parameter. This means that we treat $`Q`$ and $`b`$ as free and uncorrelated parameters. This is very conservative and eliminates any possible systematics involved in determining the bias parameter. We constrain the analysis to flat ($`\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }=1`$) models, and we assume that the tensor mode contribution is negligible ($`T/S=0`$). These assumptions are compatible with analyses of the present data WTZ , and relaxing them does not have a big effect on the final results. For maximizing the likelihood function we use a simulated annealing method, as described in reference Hannestad:2000wx .

Table I shows the different priors used. In the “CMB” prior the only important constraint is that $`0.4h0.9`$ ($`hH_0/(100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1)`$). For the CMB+$`H_0`$+BBN prior we use the constraint $`H_0=72\pm 8\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ from the HST Hubble key project freedman  (the constraint is added assuming a Gaussian distribution) and the constraint $`\mathrm{\Omega }_bh^2=0.020\pm 0.002`$ from BBN Burles:2000zk . Finally, in the $`H_0`$+BBN+LSS case, we add data from the 2dF survey THX .

Following Ref. Wandelt , we plot the constraints from various experiments in the $`(m_c,\sigma )`$ plane in Fig. 2; limits obtained with different priors are plotted with different line styles. The shaded region is the “Spergel-Steinhardt” region as given in equation (1); note that it is actually narrower than plotted in Fig. 3 of Ref. Wandelt . This figure clearly shows that constraints from the CMB data alone for large CDM mass ($`m_c>m_b1`$ GeV) are comparable to limits from galactic disk heating arguments Starkman . On the other hand, for small CDM mass ($`m_c<m_b1`$ GeV), the limit from the CMB alone is much stronger. This difference arises because the limits from disk heating are related to energy transfer, which decreases as $`(m_c/m_b)^2`$ for $`m_cm_b`$, while the limits from the CMB and large scale structure are based on momentum transfer, which decreases only as $`m_c/m_b`$ for $`m_cm_b`$. At the same time, the relative velocity also increases for low CDM mass (as $`m_c^{1/2}`$). This is the reason that, for $`m_cm_b`$, our bound goes as $`m_c^{1/2}`$, but for the disk heating argument it goes as $`m_c^1`$. Our bound is therefore much stronger at low mass, even if only CMB data is used.

The addition of priors on $`\mathrm{\Omega }_b`$ and $`H_0`$ does not increase the exclusion region by much, because there is very little degeneracy between $`(m_c,\sigma )`$ and these two parameters. Only when LSS data is added do we find a significant improvement. The reason is that CDM-baryon interactions lead to a significant suppression of small scale power, as discussed in the previous section. Even for relatively small cross sections this effect leads to a disagreement with the 2dF data. With the inclusion of LSS data, our bound becomes much stronger than the disk heating argument of Starkman et al. Starkman , even at large masses. A good analytical fit to our 95% C.L. limit is

$$\frac{\sigma }{10^{24}\mathrm{cm}^2}<\frac{0.63x^{1/2}+0.22x^{3/2}}{1+x^{1/2}}$$
(33)
where $`x=m_c`$/GeV. Equation (33) is accurate to within 10%. For this case, we also consider the effect of coherent ($`\sigma _{He}=16\sigma `$) and spin-dependent ($`\sigma _{He}=0`$) scattering between the dark matter and helium. The limits for these two cases are shown in Fig. 2 as thinner dashed curves above (for $`\sigma _{He}=0`$) and below (for $`\sigma _{He}=16\sigma `$) the limit for $`\sigma _{He}=4\sigma `$. Altering $`\sigma _{He}`$ in this way changes the upper bounds by a factor of roughly $`23`$, which is barely noticeable on the scale of this graph.

## IV Conclusions

CMB observations constrain any scattering interaction between dark matter and baryons. Our results indicate that the limit from current CMB observations is comparable to previous limits from disk heating Starkman  for masses greater than 1 GeV, and the CMB limit is much stronger at smaller masses. If we also include large-scale structure data, then our limit is more stringent that the disk-heating limit at all masses. An analytical fit of our combined CMB+LSS limit is given in Eq. (33).

Our CBM+LSS limits exclude the region discussed in Ref. Wandelt , in which self-interacting dark matter interacts with baryons with roughly the same cross section with which it interacts with itself (Fig. 2). This, by itself, does not exclude the self-interacting dark matter scenario, since the dark matter is not required to interact with baryons at all, but it does exclude models with the indicated dark matter-baryon scattering cross section.

While space and balloon-based experiments can provide tighter constraints for specific regions in the parameter space of CDM mass and cross section, (Fig. 2 and reference Wandelt ), the combination of CMB and large-scale structure appears to provide the best general upper limit on the CDM-baryon scattering cross section for arbitrary CDM masses (and it is much stronger than all other limits at low masses). Our constraints are somewhat stronger for coherent scattering from helium nuclei and weaker for a spin-dependent interaction. Perhaps more importantly, these limits will only get better with new CMB data.

###### Acknowledgements.

X.C. is supported by the NSF under grant PHY99-07949, R.J.S. is supported by the DOE under grant DE-FG02-91ER40690.


format: text

text:
# Limits to the Cas A 44Ti Line Flux and Constraints on the Ejecta Energy and the Compact Source

## 1 Introduction

Optical observations suggest that Cas A is the remnant of a Type Ib/c supernova whose progenitor was a WN Wolf-Rayet star with initial mass greater than 25 M (Fesen, 2001; Fesen & Becker, 1991; Jansen et al., 1988). The discovery of a point source at the center of the remnant by *Chandra* (Tananbaum, 1999) is strongly indicative of the creation of a neutron star or black hole in the event, and would be consistent with such a progenitor. Subsequent observations by *Chandra* (Chakrabarty et al., 2001; Murray et al., 2001) and *XMM* (Mereghetti, Tiengo, & Israel, 2001) have confirmed the point source nature of the object, but the existence of pulsations (and therefore a neutron star) is still questionable. Comparison of supernova explosion models (Woosley & Weaver, 1995; Nakamura et al., 2001) with measurements of the amount of 44Ti created, can be used to estimate the remnant core mass and, therefore, the nature of the compact object.

44Ti is produced by explosive Si burning and the freeze out from nuclear statistical equilibrium in supernovae, and is believed to be the primary source of 44Ca (Timmes et al., 1996). 44Ti decays with a 59.2$`\pm `$0.6 year half-life (Ahmad et al., 1998; Görres et al., 1998; Norman et al., 1998) to 44Sc producing two nuclear lines at 67.9 and 78.4 keV of essentially equal intensity. The 44Sc decays to the first excited state of 44Ca which then emits a gamma-ray of 1157 keV with a 3.93 hour half-life to reach the ground state. All three nuclear lines are expected to have essentially the same line strength. Given the half-life of 44Ti, the distance to Cas A (3.4 kpc; Reed et al. (1995)), the time since the supernova (317 yr for a mean observation date of 1997; Ashworth (1980)), and the observed line flux, one can estimate the amount of 44Ti created.

Attempts to detect the 44Ti lines began with the analysis of galactic scanning observations by the germanium spectrometer on *HEAO-3* (Mahoney et al., 1992). A 1$`\sigma `$ limit of 8.3$`\times `$10-5photons cm-2s-1 from a point source anywhere in the galaxy was determined for the 44Ti line flux. Initial results from *CGRO*/COMPTEL observations of Cas A announced the discovery of 1.157 MeV line flux from the 44Sc decay (Iyudin et al., 1994), and this was revised after further observations to be 3.3$`\pm `$0.6$`\times `$10-5photons cm-2s-1 (Iyudin, 1999). *CGRO*/OSSE observations of Cas A provided for a simultaneous fit to the three nuclear lines, and yielded a 99% confidence upper limit of 5.1$`\times `$10-5photons cm-2s-1 (The et al., 1995), with detection of the continuum radiation to 40 keV (The et al., 1996). With the launches of *RXTE* and *BeppoSAX* in 1995 and 1996, two powerful hard X-ray instruments became available for 44Ti line searches. The initial *RXTE*/HEXTE result was a 2.4$`\sigma `$ detection (4.3$`\pm `$1.8$`\times `$10-5photons cm-2s-1) during one long observation and no detection (-1.4$`\pm `$1.7$`\times `$10-5 photons cm-2s-1) in another (Rothschild et al., 1998a). Similarly, initial *BeppoSAX*/PDS measurements could only claim an upper limit ($``$5$`\times `$10-5photons cm-2s-1; Vink (1998)). After additional observations, Vink et al. (2001) were finally able to claim a good detection with a flux of 2.1$`\pm `$0.7$`\times `$10-5photons cm-2s-1 for the 68 and 78 keV lines from 44Ti. In this article we present a reanalysis of the *RXTE*/HEXTE observations of Cas A that were begun before the Vink et al. (2001) announcement of the detection.

The $`>`$10 keV continuum upon which the 44Ti lines lie has been fit historically to a power law with photon index $``$3 or a thermal bremsstrahlung model whose temperature is dependent upon the extent of the high energy limit to the data. Allen et al. (1997) used *ASCA*, *RXTE*, and *CGRO* data to show that the thermal models successful below 10 keV (2 Raymond-Smith components and a 6.4 keV Gaussian) fall well below the data above 10 keV. The addition of a broken power law ($`\mathrm{\Gamma }`$=1.8$`{}_{0.6}{}^{}{}_{}{}^{+0.5}`$ below 15.9$`{}_{0.4}{}^{}{}_{}{}^{+0.3}`$ keV and $`\mathrm{\Gamma }`$=3.04$`{}_{0.13}{}^{}{}_{}{}^{+0.15}`$ above) best described the data to 120 keV. Several authors (Laming, 2001; Atoyan et al., 2000; Baring et al., 1999; Gaisser, Protheroe & Stanev, 1998) have described the processes that may be involved in the photon continuum production in supernova remnants such as Cas A. A non-thermal bremsstrahlung component may become the dominant contributor at $``$100 keV, and therefore, its presence may affect the inferred photon flux from the 44Ti lines for observations with high sensitivity above the line complex.

## 2 Observations

*RXTE* observed Cas A several times for calibration purposes (ObsIds 00022, 10418, 30804, and 40806) and twice as part of an investigation into the 44Ti emission lines (ObsIds 10271 and 20253). The broad-band (2-60 keV) spectral results have been published by Allen et al. (1997), and this analysis concentrates on the 15-240 keV data from the High Energy X-ray Timing Experiment (HEXTE; Rothschild et al. (1998b)). The HEXTE is a set of eight NaI(Tl)/CsI(Na) scintillation detectors grouped into two clusters of four detectors each. The HEXTE detectors are mechanically collimated to a 1 FWHM field of view, and cover the 15-250 keV range. The full collecting area for HEXTE is 1400 cm2 (spectral capability was lost from one HEXTE detector early in the mission).

The observation dates and accumulated livetimes are given in Table 1. Since the livetime associated with the calibration observations is quite small compared to that of the two long observations, the analysis presented here will be restricted to the latter set of data.

## 3 Data Reduction and Analysis

The HEXTE data were accumulated for each cluster using a script developed at UCSD and the University of Tübingen. The cluster data were separated into on-source, off-source-plus and off-source-minus, since each HEXTE cluster collects real-time background data from two independent positions $`\pm `$1.5 to either side of the on-source position. The two clusters’ rotation axes are orthogonal to each other, and in this way four independent background regions are sampled. By comparing the two off-source positions for a given cluster, one can determine if a confusing source is contaminating one of the background regions, and if so, eliminate it from further analysis. The background observations are 75% as long as the on-source observations, since the time to move the pointing position on- and off-source comes out of the background observations. This ensures HEXTE continual on-source coverage.

The two background histograms from each cluster were tested for a confusing source by specifying the off-minus data as the source and the off-plus data as the background in the HEASARC spectral fitting program XSPEC (Arnaud, 1996). The net count rate was compared to the ideal value of 0.0, and the spectra were examined for a continuous deviation from zero at low energies, as a sign of another source. Table 2 gives the net 15-240 keV rates for the two clusters for ObsIds 10271 and 20253. Both clusters in ObsId 10271 show a negative net flux, with that for cluster A being 2$`\sigma `$ below 0.0. Investigation of the net spectra of each revealed a negative residual at 30 and 65 keV — the location of the prominent background lines. This is an indication of an under- or over-estimated livetime. The HEXTE livetime estimation is based upon a model using the Upper Level Discriminator (250 keV) and the eXtreme Upper Level Discriminator (20 MeV) rates, and on a daily timescale is accurate to $``$1%. Changing the livetime of cluster A off-plus by 2.1% and cluster B off-plus by 1.3% brought the net rates to 0.0 within uncertainties. ObsId 20253 data required no such change in livetime.

Since the HEXTE background has emission line features at 67 keV due to activation of the NaI(Tl) and at the lead K-lines at 74 and 85 keV due to the collimator, one must demonstrate that any claimed line features are not due to imperfect background subtraction. The brightest background line complex is at 30 keV, with a flux of 3.85$`\times `$10-2 photons cm-2s-1. The 90% upper limit to its flux in the net Cas A spectrum is 2.18$`\times `$10-5 photons cm-2s-1, or 0.07% of background. The 67 and 74 keV background lines have fluxes of 1.9$`\times `$10-2 photons cm-2s-1 and 1.2$`\times `$10-2 photons cm-2s-1, respectively. Using the percentage of background upper limit at 30 keV, we estimate the systematic sensitivity limit for the 44Ti lines in this observation is 1$`\times `$10-5 photons cm-2s-1.

Separate good time intervals were calculated for the three pointing positions of each cluster. The good time intervals required the pointing direction to be within 0.01 of the center of the Cas A remnant, its elevation above the Earth’s horizon to be greater than 10, and the time since the start of the most recent South Atlantic Anomaly passage to be greater than 20 minutes.

For each individual observation (17 each within ObsIds 10271 and 20253), the plus and minus off-source data sets were combined to form the background data set for each cluster, and the individual source/background data for each ObsId were summed for each cluster. Finally, the ObsId 10271 and 20253 cluster A and B data were summed to form a single on-source and single off-source spectral accumulation. The resulting pulse height histogram from the two 200 ks observations contains 226 ks livetime on-source and 170 ks of real-time background. Subsequent spectral fitting was performed using these two files, as well as for the combined data in each ObsId.

The pulse height data were binned into single, 1 keV bins to 30 keV, double width bins to 89, and quadruple width bins to 100 keV. Above that energy they were combined as 100-140, 140-180, 180-220, and 220-240 keV. Three sets of data were fitted: ObsId 10271, 20253, and their sum. The data were then fitted with a power law plus two gaussian lines at the expected energies for 44Ti$`^{44}`$Sc decay (67.9 and 78.4 keV). The intensities of the two lines were linked to be the same value, and line widths were considered: 1) narrow lines and 2) lines broadened by 2.5% to 1.7 and 2.0 keV respectively to account for the Cas A expansion velocity of 7500 km/s (Fesen, Becker & Goodrich, 1988). The results are given in Table 3. The fit to the continuum for the two independent observations yields consistent indices and fluxes, and the best fit to the summed data is a photon index $`\mathrm{\Gamma }`$=3.125$`\pm `$0.050 and 20-100 keV flux of 4.60$`\pm `$0.18$`\times `$10-11 ergs cm-2s-1keV-1. Fitting the 44Ti line flux in the two independent observations yielded fluxes at the -0.87$`\sigma `$ and 1.69$`\sigma `$ level, respectively. The best fit line flux to the combined data was 1.57$`\pm `$2.81$`\times `$10-5 photons cm-2s-1, which yielded a 90% confidence upper limit on the flux from each of the 44Ti lines of 3.6$`\times `$10-5 photons cm-2s-1. These values are consistent with those measured by BeppoSAX (Vink et al., 2001), and are unaffected by line widths from zero to 2 keV. The continuum was detected to 100 keV, thereby extending the maximum detected energy reported in Allen et al. (1997). The data plus best fit model for the combined data are shown in Figure Limits to the Cas A 44Ti Line Flux and Constraints on the Ejecta Energy and the Compact Source.

Residuals to the best fit model in the case of the ObsId 10271 observation contain a large fluctuation ($`\mathrm{\Delta }\chi ^2`$= 8.4) in the single bin at 70 keV. Since the 2 keV width of this bin is less than the detector resolution at that energy (10 keV), the deviation is considered a statistical fluctuation. This conclusion is supported by the residuals in the ObsId 20253 data, where the residual at that energy is significantly less ($`\mathrm{\Delta }\chi ^2`$= 3.3).

In order to investigate the effect of a non-power law continuum, we have tried two other continuum forms: 1) a broken power law with photon index $``$3 up to a break energy, after which the continuum flattens (to mimic the onset of the non-thermal bremsstrahlung component as shown in Figure 1 of Ellison et al. (1999)), and 2) a high energy thermal instead of the power law to represent a steepening of the spectrum at high energies. Fitting the 20253 data, which does not have the over-subtraction of background, shows that a $`\mathrm{\Gamma }`$=3 power law that breaks at 43 keV (best fit) to a $`\mathrm{\Gamma }`$=2.4 power law (best fit), does not present a better fit ($`\mathrm{\Delta }\chi ^2`$= 0.8 for loss of 2 degrees of freedom). The 90% upper limit on the 44Ti lines is reduced from 9.0 $`\times `$10-5 photons cm-2s-1 to 5.8 $`\times `$10-5 photons cm-2s-1. Using the thermal bremsstrahlung model in place of the power law results in a worse fit ($`\mathrm{\Delta }\chi ^2`$= 10.1 for the same number of degrees of freedom). This model under-estimates the flux before and after the 44Ti complex. The lines do become more prominent due to the dropping continuum with the upper limit rising to 1.07 $`\times `$10-4 photons cm-2s-1. A steepening of the continuum, as opposed to a straight power law, is not supported by the data.

## 4 Discussion

The observed 44Ti flux (Iyudin, 1999; Vink et al., 2001) has very interesting implications for the mass of the progenitor, the supernova ejecta energy, and the question of whether the compact object left after the explosion is a neutron star or a black hole. Taking the best estimate of the mean COMPTEL and BeppoSAX flux in each of the 44Ti lines to be 2.7$`\pm 0.5\times 10^5`$ photons/cm2 s, we find that the 44Ti yield from the Cas A supernova should be 1.5$`\pm 0.3\times 10^4`$ M, assuming an age of 317 yr for a 1997 mean date of the observations, a 44Ti half-life of 59.2 $`\pm `$ 0.6 yr, and a distance of 3.4 kpc. This can be compared with theoretical calculations of 44Ti production in core collapse supernovae.

First, we find that the 44Ti yield is roughly a factor of 2 or more greater than the expected (Woosley & Weaver, 1995; Nakamura et al., 2001) yields in any core collapse supernovae with ejecta energies close to $``$1$`\times `$1051 ergs, which can explain the observed properties of most (8 to 25 M) Type II supernovae. However, the observed yield is consistent with that expected (Woosley, Langer & Weaver, 1993; Nakamura et al., 2001) from Type Ib/c supernovae of Wolf-Rayet stars with initial masses greater than 25 M — implied by the optical observations (Fesen, 2001; Fesen & Becker, 1991; Jansen et al., 1988). In particular for a pre-collapse mass of $``$ 4 to 6 M, that would be expected (Woosley, Langer & Weaver, 1993) for a WN Wolf-Rayet star that explodes before reaching the WC phase, the observed yield implies an ejecta energy of $``$2$`\times `$1051 ergs. These calculations for WN Wolf-Rayet pre-collapse stars also all predict Type Ib/c supernova explosions that leave a neutron star remnant, not a black hole.

Secondly, the predicted (Nakamura et al., 2001) 44Ti/56Ni ratio for the Type Ib/c supernovae of such stars is $`1.7\times 10^3`$, which is similar to more general predictions (Woosley & Hoffman, 1991). Therefore the discrepancy between the measured 44Ti yield and the estimated brightness of the historical supernova still seems to require the added effects (Nagataki et al., 1998) of an asymmmetric explosion, which has also been supported by the recent observations of (Fesen, 2001).

Finally, the consistency between the observed 44Ti yield and that calculated for a Type Ib/c supernova of a WN Wolf-Rayet progenitor, thus would appear to rule a black hole as the point source discovered by *Chandra* (Tananbaum, 1999).

Although no pulse period has yet been detected from the Cas A compact object, comparisons of the spectral index and luminosity of the *Chandra* x-ray source with those of other x-ray pulsars suggest (Chakrabarty et al., 2001; Mereghetti, Tiengo, & Israel, 2001) that it looks more like an Anomalous X-ray Pulsar/Soft Gamma-ray Repeater (AXP/SGR) than a radio pulsar. We point out that another indicator for such an association is the fact that the Cas A remnant is expanding in the warm denser phase of the interstellar medium (e.g. Higdon & Lingenfelter (1980); Hatsukade & Tsunemi (1992)), which is also the site of the bulk of the AXP/SGRs (Marsden et al., 2001). The majority of the radio pulsars, on the other hand, are observed in the hot tenuous medium where most of the core-collapse supernovae occur. Fallback disk (Alpar, 2001) and push-back disk (Marsden et al., 2001) models for AXP/SGRs utilize ejecta material to form a small accretion disk around a conventional neutron star (i.e., with magnetic field of 1010-13 Gauss), which provides the additional torque to spin-down the neutron star rapidly and the accreted matter to generate the greater X-ray luminosity seen in AXP/SGRs. Fallback disks form within a few days from material that cannot escape the newly formed neutron star, whereas push-back disks form many years later from material decelerated by the Sedov-phase reverse shock and subsequently captured by the neutron star. We suggest Cas A is possibly in the early stages of the AXP/SGR scenario in which the push-back disk has yet to form, since the reverse shock has not reached the material in the vicinity of the neutron star. Thus, one might expect the present emission to be indicative of a conventional cooling neutron star, as per Chakrabarty et al. (2001), and when the disk does form, the accretion will increase the luminosity to that of present-day AXP/SGRs and pulsed emission will commence.

The authors would like to thank Jacco Vink, Shigehiro Nagataki, and Anatoli Iyudin for fruitful discussions, and Harish Khandrika for data processing assistance. The authors also acknowledge the useful comments by the referee. This work was supported by NASA contract NAS5-30720.


text:
# Resistance and Resistance Fluctuations in Random Resistor Networks Under Biased Percolation

## I Introduction

The study of electrical and mechanical stability of disordered systems is attracting a considerable interest in the recent literature - because of its implications on both material technology -, - and fundamental aspects related to the response of these systems to high external stresses -,-. Indeed, the application of a finite stress (electrical or mechanical) to a disordered material generally implies a nonlinear response, which ultimately leads to an irreversible breakdown (catastrophic behavior) in the high stress limit -,. Such catastrophic phenomena have been successfully studied by using percolation theories -, -,- . Critical phenomena near the percolation threshold have been widely investigated in the electrical breakdown of granular metals or conductor-insulator composites -,-, -. The associated critical exponents have been measured -,  and theoretically studied using continuum or lattice percolation models -,-. In particular, large attention has been devoted to the critical exponents describing the resistance and its relative noise in terms of the medium properties (e.g. conducting particle fraction, defect concentration, etc.) -,, -. However, very few attempts -, have been made so far to describe the behavior of a disordered medium over the full range of the applied stress, i.e. by studying the response of the system to an external bias when the bias strength covers the full range of linear and nonlinear regimes. Therefore, a satisfactory understanding of breakdown phenomena over the full dynamical regime is still missing ,. On the other hand, relevant information can be obtained from such a study, like: precursor phenomena, role of the disorder, existence of scaling laws, predictability of breakdown, etc. -,.

Here we present a percolative model of sufficient generality to address the above issues. Our aim is to provide a theoretical framework to study response and fluctuation phenomena under linear and nonlinear regimes in a wide class of disordered systems. To this purpose, we study the evolution of a random resistor network (RRN) in which two competing processes are present, defect generation and defect recovery, which determine the values of the elementary network resistances. Both processes are driven by the joint effect of an electrical bias and of the heat exchange between the network and the thermal bath. The bias is applied through a constant voltage or, alternatively, a constant current. Monte Carlo (MC) simulations are performed to investigate the network evolution in the full range of bias values. Depending on the bias strength, an irreversible failure or a stationary state of the RRN can be achieved. By focusing on the steady-state, we analyze the behavior of the average network resistance, $`<R>`$, and the properties of the resistance fluctuations as a function of the bias. In constant voltage conditions, a scaling relation is found between $`<R>/<R>_0`$ and $`V/V_0`$, where $`<R>_0`$ is the linear regime resistance and $`V_0`$ the threshold value for the onset of nonlinearity. A similar relation is found in constant current conditions. The relative variance of resistance fluctuations also exhibits a strong nonlinearity in the pre-breakdown regime whose properties are discussed. The fluctuation analysis is completed by a further investigation of the noise resistance spectrum and of the Gaussian features of the fluctuation amplitudes. Theoretical results agree with electrical breakdown measurements in thin films of composites  and of other conducting  or insulating materials .

The paper is organized as follows. In Sect. 2 we briefly describe the model used. Section 3 presents the results of the MC simulations for the resistance and its fluctuations. The main conclusions are drawn in Sect. 4.

## II Model

We study a two-dimensional random resistor network of total resistance $`R`$, made of $`N_{tot}`$ resistors, each of resistance $`r_n`$, disposed on a square lattice. We take a square geometry, $`N\times N`$, where $`N`$ determines the linear size of the lattice, with the total number of resistors being $`N_{tot}=2N^2`$. For the comparison with resistivity measurements in thin films, the value of $`N`$ can be related to the ratio between the size of the sample and the average size of the grains composing the sample. An external bias, represented by a constant voltage $`V`$ or by a constant current $`I`$, is applied to the RRN through electrical contacts realized by perfectly conducting bars at the left and right hand sides of the network. A current $`i_n`$ is then flowing through each resistor. The RRN interacts with a thermal bath at temperature $`T_0`$ and the resistances $`r_n`$ are taken to depend linearly on the local temperatures, $`T_n`$, as:

$$r_n(T_n)=r_0[1+\alpha (T_nT_0)]$$
(1)
In this expression $`\alpha `$ is the temperature coefficient of the resistance and $`T_n`$ is calculated by adopting a biased percolation model  as:

$$T_n=T_0+A\left[r_ni_n^2+\frac{B}{N_{neig}}\underset{l=1}{\overset{N_{neig}}{}}\left(r_li_l^2r_ni_n^2\right)\right]$$
(2)
Here, $`N_{neig}`$ is the number of first neighbours around the nth resistor, the parameter $`A`$, measured in (K/W), describes the heat coupling of each resistor with the thermal bath and it determines the importance of Joule heating effects. The parameter $`B`$ is taken to be equal to $`3/4`$ to provide a uniform heating in the perfect network configuration. We note that Eq. (2) implies an instantaneous thermalization of each resistor at the value $`T_n`$, therefore, by adopting Eq. (2), for simplicity we are neglecting time dependent effects which are discussed in Ref. \[\].

In the initial state of the network (corresponding to the perfect network configuration with no heating) all the resistors are identical: $`r_nr_0`$. Now, we assume that two competing percolative processes act to determine the RRN evolution. The first process consists of generating fully insulating defects (resistors with very high resistance, i.e. broken resistors) with probability  $`W_D=\mathrm{exp}(E_D/K_BT_n)`$, where $`E_D`$ is an activation energy characteristic of the defect and $`K_B`$ the Boltzmann constant. The second process consists of recovering the insulating defects with probability $`W_R=\mathrm{exp}(E_R/K_BT_n)`$, where $`E_R`$ is an activation energy characteristic of this second process. For $`A0`$, Eq. (2) implies that both the processes (defect generation and defect recovery) are correlated percolations. Indeed, the probability of breaking (recovering) a resistor is higher in the so called ”hot spots” of the RRN . On the other hand, for $`A=0`$ Eq. (2) yields $`T_nT_0`$, which corresponds to random percolations . The same is true for vanishing small bias values, when Joule heating effects are negligible.

As a result of the competition between these two percolations and depending on the parameters related to the particular physical system ($`E_D`$, $`E_R`$, $`A`$, $`\alpha `$, $`r_0`$, $`N`$) and on the external conditions (specified by the bias conditions and the bath temperature), the RRN reaches a steady state or exhibits an irreversible breakdown. In the first case, the network resistance fluctuates around an average value $`<R>`$. In the second case, a critical fraction of defects $`p_c`$, corresponding to the percolation threshold, is reached, i.e. $`R`$ diverges due to the existence of at least one continuous path of defects between the upper and lower sides of the network . We note that in the limit of a vanishing bias (random percolation) and infinite lattices ($`N\mathrm{}`$), the expression: $`E_R<E_D+K_BT_0\mathrm{ln}[1+\mathrm{exp}(E_D/K_BT_0)]`$ provides a sufficient condition for the existence of a steady state .

The evolution of the RRN is obtained by MC simulations carried out according to the following iterative procedure. (i) Starting from the perfect lattice with given local currents, the local temperatures $`T_n`$ are calculated according to Eq. (2); (ii) the defects are generated with probability $`W_D`$ and the resistances of the unbroken resistors are changed as specified by Eq. (1); (iii) the currents $`i_n`$ are calculated by solving Kirchhoff’s loop equations by the Gauss elimination method and the local temperatures are updated; (iv) the defects are recovered with probability $`W_R`$ and the temperature dependence of unbroken resistors is again accounted for; (v) $`R`$, $`i_n`$ and $`T_n`$ are finally calculated and the procedure is iterated from (ii) until one of the two following possibilities is achieved. In the first, the percolation threshold is reached. In the second, the RRN attains a steady state; in this case the iteration runs long enough to allow a fluctuation analysis to be carried out. Each iteration step can be associated with an elementary time step on an appropriate time scale (to be calibrated with experiments). In this manner it is possible to represent the simulation of the resistance evolution over either a time or a frequency domain according to convenience.

In the simulations, as reasonable values of the parameters, we have taken: $`N=75`$, $`r_0=1(\mathrm{\Omega })`$, $`\alpha =10^3`$ (K-1), $`A=5\times 10^5`$ (K/W), $`E_D=0.17`$ (eV), $`E_R=0.043`$ (eV), and $`T_0=300`$ (K) if not stated otherwise. The values of the external bias range from $`0.05I2.8`$ (A) under constant current conditions, and from $`0.05V3.5`$ (V) under constant voltage conditions.

## III Results

Figures 1 (a) and 1 (b) report a sampling of the resistance evolutions coming out from the simulations under steady state conditions for the cases of constant voltage (Fig. 1 (a)) and constant current (Fig. 1 (b)), respectively. At increasing the external bias, a remarkable increase of the average resistance and of the amplitude of fluctuations is evident in both cases while approaching the breakdown. This is reached by applying an external bias just above the highest values shown in these figures. Precisely, the RRN becomes unstable under constant voltage already by applying a voltage just above $`3.0`$ (V) and, under constant current, by applying a current just above $`2.1`$ (A).

In the following, Figs. 2 to 4 will detail the behavior of the average resistance while Figs. 5 to 10 will focus on the results of resistance fluctuations.

### A Resistance

Figure 2 reports the average value of the RRN resistance as a function of the applied bias (current or voltage). Each value is calculated by considering the time average on a single steady-state realization and then averaging over 20 independent realizations (RRNs subjected to the same bias conditions). The numerical uncertainty is found to be within 0.01% at worst. At the lowest biases the resistance takes a value $`<R>_0`$ which represents the intrinsic linear response property of the network (Ohmic regime) . We note that for the activation energy values here chosen, the average fraction of defects at the lowest bias $`<p>_0`$ is very small ($`<p>_0=0.545\times 10^2`$) so that $`<R_0>=0.9978`$ ($`\mathrm{\Omega }`$) is about 1% greater than the perfect network resistance. Under constant current conditions, when the current overcomes a certain value, $`I_0`$, the average resistance starts to become dependent on the bias. Thus, $`I_0`$ sets the current scale value for the onset of nonlinearity. By adopting as criterium for the onset of nonlinearity a resistance increase of $`0.05\%`$ over the linear response value, we determine $`I_0=0.090\pm 0.005`$ (A). For currents above $`I_0`$ the resistance increase is smooth at moderate bias and it exhibits a sharp ramp (typical of a catastrophic behavior) at high bias, until a threshold current value $`I_b=2.10`$ (A) above which the RRN undergoes an irreversible breakdown. A step of current values $`\delta I=0.05`$ (A) has been used to determine $`I_b`$, thus this value of $`\delta I`$ gives the maximum uncertainty on $`I_b`$. Accordingly, we have found $`I_b/I_0=23\pm 1`$ and $`<R>_b/<R>_0=1.6\pm 0.1`$, where $`<R>_b`$ is the last stable value of the resistance calculated before the breakdown. We remark that the uncertainty on the ratio $`<R>_b/<R>_0`$ is mainly due to the uncertainty on $`<R>_b`$ which reflects, amplified by the nonlinearity, the uncertainty on $`I_b`$. Similarly, under constant voltage conditions we have found $`V_b/V_0=32\pm 2`$ and $`<R>_b/<R>_0=1.5\pm 0.1`$, where $`V_0=0.095\pm 0.005`$ (V) and $`V_b=3.00`$ (V) are the voltage values corresponding respectively to the nonlinearity onset and to the electrical breakdown. The maximum uncertainty on $`V_b`$ is $`\delta V=0.05`$ (V). The significantly higher value of the ratio $`V_b/V_0`$, when compared with that of the ratio $`I_b/I_0`$, is a quantitative indication that the system is more robust when biased under constant voltage than under constant current conditions. This property is further emphasized by the fact that the increase of the resistance in the pre-breakdown region exhibits a lower slope under constant voltage than under constant current conditions. It must be noticed that in spite of the significant difference of the ratios $`V_b/V_0`$ and $`I_b/I_0`$, the ratio $`<R>_b/<R>_0`$ remains the same, within the error, under the different bias conditions. This result agrees with measurements of the ratio $`<R>_b/<R>_0`$ performed in composites under the Joule regime .

To better analyze the dependence on the bias of the average resistance, Figs. 3 (a) and 3 (b) report the log-log plot of the relative variation of the resistance as a function of $`V/V_0`$ and $`I/I_0`$, respectively. Figure 3 (a) shows that the relative variation of the average resistance scales with the ratio $`V/V_0`$ in the whole region of applied voltages up to breakdown as:

$$\frac{<R>_V}{<R>_0}=1+a\left(\frac{V}{V_0}\right)^\theta $$
(3)
where $`a=(1.1\pm 0.1)\times 10^4`$ is a dimensionless coefficient and $`\theta =2.1\pm 0.1`$. Figure 3 (b) shows that also under constant current the relative variation of the average resistance scales with the ratio $`I/I_0`$ in the moderate bias region, as:

$$\frac{<R>_I}{<R>_0}=1+a^{}\left(\frac{I}{I_0}\right)^\theta $$
(4)
with $`a^{}=(3.8\pm 0.1)\times 10^4`$ and the same exponent $`\theta `$ of Eq. (3), within numerical uncertainty. However, in this case we notice that in the pre-breakdown region the relative variation of the average resistance exhibits a superquadratic behavior characterized by a power law $`(I/I_0)^{\theta _I}`$ with $`\theta _I=3.7\pm 0.1`$.

This behavior of the average resistance with the applied bias can be understood as follows. From a macroscopic point of view, in a degradation process associated with a resistance increase, constant current conditions lead, at increasing bias, to a superquadratic increase of the power dissipated through Joule heating and thus to a major efficiency in the defect generation. By contrast, constant voltage conditions lead, at increasing bias, to a growth of the dissipated power which is subquadratic in the applied voltage, thus implying a minor efficiency in the generation of defects. In this respect, it must be underlined that the resistance increase driven by both kind of bias is further enhanced by the positive value of the temperature coefficient of the resistance.

From a microscopic point of view, the RRN resistance depends on the average fraction of defects, $`<p>`$, according to the expression :

$$<R>|
p_c|^\mu $$
(5)
where, for biased percolation, $`p_c`$ and $`\mu `$ are functions of the bias strength . Of course, for vanishing bias, $`p_c`$ and $`\mu `$ take the well known values corresponding to random percolation: $`p_c=0.5`$ (for a square lattice) and $`\mu =1.303`$ (universal value) . On the other hand, the dependence of $`p_c`$ and $`\mu `$ on the bias makes the analysis of $`<R>`$ in terms of Eq. (5) quite problematic. Nevertheless, it is interesting to analyze the behavior of $`<p>`$ on the bias strength and, to this purpose, Fig. 4 reports on a log-log plot the relative variation of $`<p>`$ as a function of the normalized bias. As evidenced by the figure, for the case of constant voltage conditions the average fraction of defects scales with the same exponent $`\theta `$ over the full range of bias values, as:

$$\frac{
_V}{
_0}=1+b\left(\frac{V}{V_0}\right)^\theta $$
(6)
with $`b=(5.8\pm 0.3)\times 10^3`$ a dimensionless coefficient. A similar expression holds under constant current conditions at moderate bias:

$$\frac{
_I}{
_0}=1+b^{}\left(\frac{I}{I_0}\right)^\theta $$
(7)
with $`b^{}=(5.3\pm 0.3)\times 10^3`$. However, under constant current, there is a significant increase of the slope in the pre-breakdown region up to a nearly above cubic power law.

The behavior of the average fraction of defects in the moderate bias region can be understood by generalising the results obtained for random stationary percolations  to the presence of a moderate external bias responsible of Joule heating. In fact, for a RNN in a stationary state resulting from the competition of two random percolations, it has been found  that $`<p>=x/(1+x)`$, with $`x=W_{0D}(1W_{0R})/W_{0R}`$, where $`W_{0D}`$ and $`W_{0R}`$ are the probabilities defined in section II with $`T_nT_0`$. At moderate bias, when the defect distribution is rather homogeneous and the variation of $`<R>`$ small, we can take $`x=x(T)`$ with $`T=T_0+\mathrm{\Delta }T`$, where in the spirit of a mean-field theory, the average temperature increase can be expressed as $`\mathrm{\Delta }T=A<R>_0I^2/(2N^2)`$. Therefore, by differentiating:

$$\mathrm{\Delta }
=\frac{1}{(1+x)_0^2}\left(\frac{dx}{dT}\right)_0\mathrm{\Delta }T$$
(8)
which by simple manipulation gives:

$$\frac{\mathrm{\Delta }
}{
_0}=\left[\frac{A<R>_0I_0^2}{2N^2
_0[(1+x)^2]_0}\left(\frac{dx}{dT}\right)_0\right]\left(\frac{I}{I_0}\right)^2$$
(9)
The term in the squared brackets of Eq. (9) can thus be identified with the $`b^{}`$ coefficient in Eq. (7). In fact, by introducing the parameters used in the simulations we find for the squared brackets in Eq. (9) the value $`5.5\times 10^3`$ which well agrees with the previous reported value of $`b^{}`$. For constant voltage the factor $`A<R>_0I_0^2`$ inside Eq. (9) must be replaced with $`AV_0^2/<R>_0`$. This gives for the square bracket in Eq. (9) the value of $`6.1\times 10^3`$, which again well agrees with the value of $`b`$ in Eq. (6). Therefore, at moderate bias, the quadratic behavior of $`<p>`$ is well explained by Eq. (9). In the pre-breakdown region, under constant current conditions, the significant increase of the resistance and the defect filamentation characteristic of biased percolation  add to give a superquadratic increase of $`<p>`$ with the bias. By contrast, under constant voltage, these two effects act in an opposed manner, thus their compensation results in a quadratic dependence of $`<p>`$ over the full range of bias values.

From the above behavior of the average defect fraction we can now understand the bias dependence of $`<R>`$. At moderate bias ($`I/I_0`$ or $`V/V_010`$) the increase of $`<R>`$ is rather small. Accordingly, super and subquadratic effects on the dissipated power are negligible and the average fraction of defects over the intrinsic value grows quadratically with the bias strength (Eq. (9)); at these bias values $`<p>1`$ thus $`<R><p>`$ and also the average resistance grows quadratically with bias. At high bias (i.e. $`I/I_0`$ or $`V/V_0>10`$), when entering the pre-breakdown region, the increase of $`<R>`$ is significant. Accordingly, super and subquadratic effects on the dissipated power are relevant as well as the increase of the average fraction of defects which now takes a nonuniform distribution (biased percolation). As a consequence, the percolation threshold is lowered with respect to the random percolation value . Depending on constant current or constant voltage conditions, this biased percolation effect adds to further enhancing the increase of the resistance well above the quadratic law or to keep the quadratic behavior, respectively.

This interpretation is confirmed by MC simulations performed with the same parameters but setting the temperature coefficient of the resistance $`\alpha =0`$. In this case, we have found that the onset of nonlinearity as well as the electrical breakdown occur at higher values of the bias. Precisely we have found $`I_0^{}=0.25\pm 0.03`$ (A), $`I_b^{}=2.65\pm 0.05`$ (A) (thus $`I_b^{}/I_0^{}=11\pm 1`$), $`<R>_b/<R>_0=1.2\pm 0.1`$ under constant current and $`V_0^{}=0.25\pm 0.03`$ (V), $`V_b^{}=3.25\pm 0.05`$ (V) (thus $`V_b/V_0=13\pm 1`$) and $`<R>_b/<R>_0=1.3\pm 0.1`$ under constant voltage conditions. The quadratic regime is always present at moderate bias values. However, in the pre-breakdown regime we have found $`\theta _I^{}=3.1\pm 0.1`$ for constant current and $`\theta _V^{}=2.8\pm 0.1`$ for constant voltage conditions, respectively.

We finally conclude this section by remarking that the behavior described by Eqs. (3) and (4) well agrees with measurements of the resistance of composites made in the Joule regime up to breakdown .

### B Fluctuations

It is well known that important information about the properties and the stability of different systems (physical, biological, social, etc.) can be obtained by studying the fluctuations around the average value of some characteristic quantities of the system -,-, -,-. Several features of the fluctuations are usually considered, and among these, the most important are the variance (or the relative variance), the power spectrum and the Gaussianity property. Here, we analyze the fluctuations of the RNN resistance in the full range of bias values, devoting particular attention to the features that can be identified as failure precursors. Figure 5 reports the relative variance of resistance fluctuations, $`\mathrm{\Sigma }<\delta R^2>/<R>^2`$, as a function of the applied bias. The data correspond to the same simulations presented in Fig. 2. The same procedure of time averaging over a single simulation and then making ensemble averages over $`20`$ realizations provides an uncertainty of $`3\%`$ at worst. At low bias, the relative variance is found to achieve a constant value $`\mathrm{\Sigma }_0=(4.0\pm 0.1)\times 10^6`$. This value represents an intrinsic property of the system determined by the competition of two random percolations, as already discussed for the case of the resistance. A detailed analysis of the scaling properties of $`\mathrm{\Sigma }_0`$ can be found in Ref. \[\]. At increasing bias, when $`I>I_0`$ or $`V>V_0`$, the systematic increase of $`\mathrm{\Sigma }`$ reveals the existence of a nonlinear regime of the response. Similarly to the behavior of the resistance, by approaching the electrical breakdown $`\mathrm{\Sigma }`$ exhibits a significant increase which is steeper under constant current than under constant voltage conditions. This feature is ascribed to the better stability of the RRN under constant voltage, since resistance fluctuations in excess over the mean value are damped in this condition. By contrast, the same kind of fluctuations are enhanced under constant current conditions. To detail the dependence of $`\mathrm{\Sigma }`$ on the bias, we have reported in Fig. 6 (a) its relative variation, $`[\mathrm{\Sigma }_V\mathrm{\Sigma }_0]/\mathrm{\Sigma }_0`$, as a function of $`V/V_0`$ under constant voltage conditions. Two regions can be identified in the nonlinear regime: a moderate bias region ($`V/V_0<10`$), characterized by a quadratic dependence on the applied voltage, and a pre-breakdown region at the highest voltages ($`V/V_0>10`$), where a superquadratic dependence up to a cubic power law is evidenced. Accordingly, in the first region it is:

$$\frac{\mathrm{\Sigma }_V}{\mathrm{\Sigma }_0}=1+c\left(\frac{V}{V_0}\right)^2$$
(10)
with $`c=(1.14\pm 0.06)\times 10^2`$ a dimensionless coefficient, while in the second region we have found:

$$\frac{\mathrm{\Sigma }_V}{\mathrm{\Sigma }_0}\left(\frac{V}{V_0}\right)^{\eta _V}$$
(11)
where $`\eta _V=3.0\pm 0.1`$. This behavior of $`\mathrm{\Sigma }`$ under constant voltage conditions should be compared with that of the average resistance which remains always a quadratic function up to the breakdown, as shown in Fig. 3 (a). The emergence of a superquadratic behavior of $`\mathrm{\Sigma }`$ in the pre-breakdown region reflects the higher sensitivity of the resistance fluctuations with respect to the average value of the resistance to the instability of the network. The same analysis performed for the case of constant current is shown in Fig. 6 (b). Again, at moderate bias $`[\mathrm{\Sigma }_I\mathrm{\Sigma }_0]/\mathrm{\Sigma }_0`$ increases quadratically with the current as:

$$\frac{\mathrm{\Sigma }_I}{\mathrm{\Sigma }_0}=1+c^{}\left(\frac{I}{I_0}\right)^2$$
(12)
with $`c^{}=(6.6\pm 0.4)\times 10^3`$ a dimensionless coefficient. On the other hand, the superquadratic dependence characterizes the pre-breakdown region as:

$$\frac{\mathrm{\Sigma }_I}{\mathrm{\Sigma }_0}\left(\frac{I}{I_0}\right)^{\eta _I}$$
(13)
with $`\eta _I=5.4\pm 0.1`$. We note, that $`\eta _I`$ is significantly greater than $`\eta _V`$, according to the behaviors shown in Fig. 5. Moreover, $`\eta _I`$ is greater than the resistance exponent in the pre-breakdown region $`\theta _I`$. It is noteworthy, that by neglecting the effect of the temperature coefficient of the resistance, i.e. by taking $`\alpha =0`$, the moderate bias region remains characterized by a a quadratic increase which is common to both constant current and constant voltage conditions. However, in the pre-breakdown region we have found $`\eta _V=\eta _I=4.0\pm 0.1`$. Thus, the above results prove that the quadratic dependence of $`\mathrm{\Sigma }`$ on the bias in the moderate bias region is a feature independent of the conditions on which the bias is applied and also independent of the value of the thermal coefficient of the resistance. This reflects the fact that at moderate bias the two following conditions are satisfied: i) the variation of the average fraction of defects follows Eq. (9); ii) the average fraction of defects is much smaller than the percolation threshold value, i.e. $`<p>p_c`$. Therefore, in this nearly perfect regime $`\mathrm{\Sigma }<p><R>`$ . On the other hand, the superquadratic behavior of $`\mathrm{\Sigma }`$ in the pre-breakdown region is a consequence of the fact that at high biases the local correlations in the defect generation and recovery processes become strong. Defect filamentation, characteristic of biased percolation , implies that a small variation in the number of defects bring strong resistance fluctuations.

Figure 7 reports the relative variance of resistance fluctuations normalized to the linear regime value as a function of the normalized resistance. We can see that the former quantity exhibits an increase over nearly two orders of magnitude in the same range of bias where the latter quantity increases for about $`60\%`$ only. The log-log plot of $`[\mathrm{\Sigma }\mathrm{\Sigma }_0]/\mathrm{\Sigma }`$ as a function of $`[<R><R>_0]/<R>_0`$, shown in Fig. 8, confirms the proportionality of $`\mathrm{\Sigma }`$ with $`<R>`$ in the moderate bias region and evidences the following power-law in the pre-breakdown region.

$$\frac{[\mathrm{\Sigma }\mathrm{\Sigma }_0]}{\mathrm{\Sigma }_0}\frac{[<R><R>_0]}{<R>_0}^\zeta $$
(14)
where the value of the exponent, $`\zeta =1.5\pm 0.1`$, within the statistical uncertainty is the same for constant voltage or constant current conditions. This result is consistent with Eqs. (3) and (4).

To complete the study of resistance fluctuations, we have investigated the Gaussian properties of the fluctuation amplitudes and the spectra in the frequency domain. Figures 9 (a) and 9 (b) report the distribution function of the resistance fluctuations $`p(R)`$ for different currents and different bath temperatures, respectively. The dotted and the long-dashed curves in Fig. 9 (a) correspond to fit with a Gaussian distribution the data corresponding to different current values in the pre-breakdown region at $`T_0=300`$ (K), while in Fig. 9 (b) to fit the data corresponding to different bath temperatures at $`I=1.5`$ (A). When approaching the breakdown conditions, at increasing current in Fig. 9 (a) and at increasing temperature in Fig. 9 (b), the simulations show the onset of a non-Gaussian behavior characterized by the enhancement of the probability for the positive fluctuations with respect to the Gaussian distribution. The emergence of a non-Gaussian behavior near the breakdown, in agreement with experiments ,, can be considered a relevant precursor of failure.

Figures 10 (a) and 10 (b) show the spectral densities of resistance fluctuations for the same conditions of Fig. 9 (a) and 9 (b). The spectral densities have been calculated by Fourier transforming the corresponding correlation functions $`C_{\delta R}(t)`$. We have found Lorentzian spectra in all the cases. For a given temperature, within the numerical uncertainty, the corner frequency is found to be independent of the applied current, while the value of the plateau increases more than quadratically with the current in the pre-breakdown region. For a given current in the pre-breakdown region, both the corner frequency and the value of the plateau are found to increase at increasing temperatures. This fact indicates that the characteristic times of fluctuations while depending on temperature they are independent of the applied bias. We finally remark that the power spectral density in Figs. 10 (a) and 10 (b) are in good agreement with experiments ,.

## IV Conclusions

We have studied the stationary state of a two-dimensional RRN resulting from the competition between two biased percolations. The two percolative processes consist of the breaking and recovering of elementary resistors and they are driven by the joint effects of an electrical bias and of the heat exchange with a thermal bath. The electrical bias is set up by applying a constant voltage or a constant current. MC simulations have been performed to analyze the network resistance and its fluctuation properties in the full range of bias values, covering linear and nonlinear regimes up to the breakdown limit. The nonlinear regime starts for biases greater than the threshold values, $`V_0`$ or $`I_0`$ (nonlinearity onset values), and it extends until the values $`V_b`$ or $`I_b`$ (breakdown values). The ratios $`V_b/V_0`$, $`I_b/I_0`$, $`<R>_b/<R>_0`$ can be considered as relevant indicators characterizing the breakdown properties of the system. We have found that the ratio $`<R>_b/<R>_0`$ is independent of the bias conditions but it depends on the thermal coefficient of the resistance. This result agrees with measurements of this ratio performed in composites . Moreover, we have found that under constant voltage conditions the relative variation of the average resistance scales quadratically with the ratio $`V/V_0`$ over the full nonlinear regime. A similar relation has been found under constant current conditions, but in this case a superquadratic dependence emerges in the pre-breakdown region ($`I>10I_0`$). For what concerns the relative variance of resistance fluctuations we have found that at moderate bias ($`I_0<I<10I_0`$) it grows quadratically with the bias, independently of constant-current or constant-voltage conditions. On the other hand a superquadratic dependence appears in the pre-breakdown region. The presence of two distinct regions in the nonlinear regime has been explained in the following terms. At moderate bias there is a small Joule heating giving rise to the competitions of two nearly random percolations and to a small variation of the RRN resistance. At high bias, local correlations typical of biased percolation become important. Moreover, also the resistance variation becomes strong. Therefore, as a consequence of both these effects the resistance noise depends superquadratically on the bias. The power spectral density of resistance fluctuations is found to exhibit a Lorentzian spectrum and in the pre-breakdown region the amplitude of fluctuations shows the onset of non-Gaussian behaviors. Theoretical results agree qualitatively, and in some cases quantitatively, with breakdown experiments performed in composites , and in other conducting , or nonconducting materials . We conclude that the study of the stationary state of a RRN, resulting from the competition of biased percolations, provides an unifying framework for the interpretation of several non-linear transport phenomena in a variety of disordered systems.

ACKNOWLEDGEMENTS

This research is performed within the STATE project of INFM. Partial support is also provided by the MADESS II project of the Italian National Research Council and the ASI project, contract I/R/056/01. We thank Prof. L.B. Kish and Dr. Z. Gingl who introduced us to the percolation theory.

FIGURE CAPTIONS

Fig. 1 (a) - Resistance evolutions of a steady state $`75\times 75`$ network under constant voltage conditions for increasing values of the bias. Going from the lower to the upper curves, the voltage is $`V=0.2`$, $`0.5`$, $`1.0`$, $`1.5`$, $`2.0`$, $`3.0`$ (V). The values of all the other parameters are specified in the text.

Fig. 1 (b) - Resistance evolutions of a steady state $`75\times 75`$ network under constant current conditions for increasing values of the bias. Going from the lower to the upper curves, the current is: $`I=0.2`$, $`0.5`$, $`1.0`$, $`1.5`$, $`2.1`$ (A). The values of all the other parameters are specified in the text.

Fig. 2 - Average network resistance as a function of the external bias. The averages are calculated over $`20`$ networks of sizes $`75\times 75`$ subjected to the same bias conditions. The circles correspond to biasing the network by an external constant current, the triangles correspond to a constant voltage bias.

Fig. 3 (a) - Log-Log plot of the relative variation of the network resistance calculated under constant voltage as a function of the ratio $`V/V_0`$. The data are the same of those in Fig. 2 and only the nonlinear regime is shown. The dotted line represents the power-law fit with a value for the exponent $`\theta =2.1\pm 0.1`$.

Fig. 3 (b) - Log-Log plot of the relative variation of the network resistance calculated under constant current as a function of the ratio $`I/I_0`$. The data are the same of those in Fig. 2 and only the nonlinear regime is shown. The long-dashed line represents the power-law fit of the data in the moderate bias region, the value of the exponent is $`\theta =2.1\pm 0.1`$. The solid line fits the data in the pre-breakdown region with a power-law, the value of the exponent is $`\theta _I=3.7\pm 0.1`$.

Fig. 4 - Log-Log plot of the relative variation of the average number of defects under constant voltage (triangles) and under constant current conditions (open circles) as a function of the normalized bias ($`V/V_0`$ or $`I/I_0`$). The solid line represents the power-law fit with a value for the exponent $`\theta =2.1\pm 0.1`$. The long-dashed line fits constant-current data in the pre-breakdown region with a power-law of exponent $`3.3\pm 0.1`$.

Fig. 5 - Relative variance of resistance fluctuations, $`\mathrm{\Sigma }`$, under constant voltage conditions (triangles) and under constant current conditions (open circles) as a function of the external bias. Each point has been obtained by averaging the variances calculated for $`20`$ networks of sizes $`75\times 75`$ subjected to the same bias conditions.

Fig. 6 (a) - Log-Log plot of the relative variation of $`\mathrm{\Sigma }`$ under constant voltage conditions as a function of $`V/V_0`$. The long-dashed line represents the power-law fit to the data in the moderate bias region, the value of the exponent is $`2.0\pm 0.1`$. The solid line fits data in the pre-breakdown region with a power-law of exponent $`3.0\pm 0.1`$.

Fig. 6 (b) - Log-Log plot of the relative variation of $`\mathrm{\Sigma }`$ under constant current conditions as a function of $`I/I_0`$. The long-dashed line represents the power-law fit to the data in the moderate bias region, the value of the exponent is $`2.1\pm 0.1`$. The solid line fits data in the pre-breakdown region with a power-law of exponent $`5.4\pm 0.1`$.

Fig. 7 - Relative variance of resistance fluctuations normalized to the linear regime value as a function of the normalized average resistance. Open circles refer to constant current condition and triangles to constant voltage condition, respectively.

Fig. 8 - Log-Log plot of the relative variation of $`\mathrm{\Sigma }`$ as a function of the relative variation of the average resistance under constant voltage condition (open circles) and constant voltage condition (triangles). The long-dashed line gives the best fit with a power-law with exponent $`1.0\pm 0.1`$. The solid and the dotted lines fit with a power-law the data in the pre-breakdown region, the exponent are $`1.5\pm 0.1`$ and $`1.6\pm 0.1`$, respectively.

Fig. 9 (a) - Distribution function of the resistance fluctuations for two values of the applied current. Open circles refer to $`I=1.5`$ (A) and full diamonds to $`I=1.8`$ (A). The scale is a linear-log, therefore the dashed and the long-dashed curves, coming from parabolic fits of the two sets of data, correspond to Gaussian distributions.

Fig. 9 (b) - Distribution function of the resistance fluctuations under constant current conditions for two values of the substrate temperature. Open circles refer to $`T_0=300`$ (K) (the same data reported in Fig. 9 (a) ) and full triangles to $`T_0=450`$ (K). The scale is linear-log, therefore the dashed and the long-dashed curves, coming from parabolic fits of the two sets of data, correspond to Gaussian distributions.

Fig. 10 (a) - Power spectral density of resistance fluctuations under constant current conditions for $`I=1.5`$ (A) (dotted curve) and for $`I=1.8`$ (A) (solid curve). In both cases the substrate temperature is $`300`$ (K).

Fig. 10 (b) - Power spectral density of resistance fluctuations under constant current condition for $`T_0=300K`$ (dotted curve) and for $`T_0=450`$ (K) (solid curve). In both cases the bias current is $`1.5`$ (A).