text:
# IR Detection of Low-Mass Secondaries in Spectroscopic Binaries

## 1 Introduction

This paper describes a technique to measure the radial velocities of faint secondaries residing in known single-lined spectroscopic binaries (SB1s) using high-resolution infrared (IR) spectroscopy. The paper outlines the approach and reports on the detection of low-mass secondaries in three main-sequence binaries. From the radial velocities of the faint secondaries, combined with the already known velocities of the primaries, we derive the mass ratios of these binaries, and estimate the masses of the secondaries.

In most of the known spectroscopic binaries only the primary spectrum is detected and measured, rendering such systems as SB1s. Measurements of SB1s yield only the mass function, in which the mass of the primary, secondary, and inclination ($`M_1`$, $`M_2`$, and $`i`$) are inseparable. To derive the mass ratios of these systems it is necessary to turn them into double-lined binaries (SB2s). To obtain the masses we further need to obtain information about the inclination, either from the light curves in eclipsing binaries, or from astrometric orbits. When the inclination of an SB2 is not known, it is possible to derive the secondary mass if the primary mass can be estimated from the system’s photometry and spectrum.

Accurate mass ratios and masses are important for many aspects of the astrophysics of binaries. We are applying the technique we describe here in two areas of particular interest to us. Accurate masses for pre$``$main-sequence (PMS) stars with low masses (M$`<`$1 M) are known for only very few objects. As a result, calculations of PMS evolution are incompletely tested. Comparison of evolutionary tracks published in the past 5 years show that, for a given location of a star in the HR diagram, the different tracks can produce a scatter of 2$``$3 in the inferred mass and age (Simon et al., 2000). One of our goals is, therefore, to improve the calibration of the theoretical tracks by deriving the mass ratios, and, eventually, the masses of the low-mass components of PMS binaries.

Another of our goals is to derive the distribution of mass ratios and secondary masses in main-sequence binaries. The secondary-mass distribution touches upon one of the unanswered basic questions of star formation, namely how binaries, and close binaries in particular, are formed. The secondary distribution provides one way to distinguish (Boss, 1993) between the star formation models (e.g. Clark 1995; Bate 2000). In addition, the lower end of the secondary mass distribution is of crucial importance to an issue that has attracted recent attention: the distinction between low-mass secondaries and extrasolar planets (Han et al., 2001; Pourbaix, 2001; Pourbaix & Arenou, 2001; Zucker & Mazeh, 2001a). Some recent studies (Basri & Marcy, 1997; Mazeh, 1999; Halbwachs et al., 2000; Jorissen et al., 2001; Zucker & Mazeh, 2001b) suggested that the two populations could be separated by comparing their mass distributions.

Spectroscopic binaries have been studied mostly in the visible, where sensitive detectors have been available for a few decades. Unfortunately, detection of secondaries in visual spectra favors systems with light ratio, and hence mass ratio, close to unity. Most of the known spectroscopic binaries are thus SB1s and the distribution of mass ratios in the known double-lined systems (SB2s) is biased toward unity.

Our approach is based on the premise that in each binary the secondary has a lower mass and is therefore cooler than the primary, favoring its detection in the IR. Consider, for example, a binary with 0.2 and 1.0 M stars at the age of 5 Gyr. The results of Baraffe et al. (1998) and Henry & McCarthy (1990) indicate that the flux ratio between the primary and the secondary is about $`8\times 10^4`$ at 5500 Å, but only 0.01 at 1.6 $`\mu `$m. Thus, the advent of high-resolution IR spectrometers with efficient large format detectors offers the means to complement the visible radial-velocity studies with IR observations that are considerably more sensitive to low-mass companions. This enables the measurement of the secondary radial velocities which convert the SB1s into SB2s. The derivation of mass ratios by this approach does not require extensive observations since the parameters of the primary orbit have already been derived from previous work so only a small number of secondary velocities are required.

To analyze the composite target spectra we use TODCOR, a two-dimensional correlation algorithm developed by Zucker & Mazeh (1994). The algorithm assumes that the observed spectrum is a combination of two known spectra, shifted by the radial velocities of the two components, and calculates the correlation of the observed spectrum against a combination of two templates with different shifts. TODCOR generates a two-dimensional correlation function, whose peak simultaneously identifies the radial velocities of both the primary and the secondary (e.g. Zucker & Mazeh 1994). One of the advantages of TODCOR is its ability to use different templates for the primary and the secondary. When the primary and the secondary are of different spectral types, the use of two different templates enables us to utilize all the spectral information contained in the observed spectrum of the combined system. This feature of TODCOR is most important when deriving the radial velocity of a faint secondary.

This paper outlines the new approach and reports on the detection of low-mass secondaries in three main-sequence binaries: G 147-36, G 164-67, and HD 144284. The former two systems were identified and studied as SB1s by Latham et al. (2001), who searched for spectroscopic binaries within a sample of high-proper-motion stars (Carney & Latham, 1987). The latter was found by Duquennoy & Mayor (1991), who studied the nearby G stars. The three stars presented here were the three most significant detections in our main-sequence program to date. Preliminary reports on the application of this technique to main-sequence and PMS SBs appear in Mazeh et al. (2000), Prato et al. (2001a), Steffen et al. (2001), and Prato et al. (2001b). Section 2 presents four approaches to the derivation of mass ratios from the primary and secondary radial velocities. Our observations are described in Section 3 and the spectral type templates in Section 4. We present our derivation of the radial velocities in Section 5, while Section 6 shows the calculation of the mass ratios for the three binaries. A short discussion of our results appears in Section 7.

## 2 The Derivation of the Mass Ratios

In this section we present four different approaches to derive the mass ratio $`q=M_2/M_1`$ of a spectroscopic binary from a few radial velocities of both the primary and the secondary. Choosing between these approaches depends upon the extent of the information available from the previous SB1 solution.

1. When no previous orbital information is available, any two pairs of velocities can yield the mass ratio of a binary (Wilson 1941). Suppose we have observed the system at two epochs, $`t_i`$ and $`t_j`$, with velocities $`\{V_{1,k},V_{2,k};k=i,j\}`$ where $`V_{1,k},V_{2,k}`$ are the velocities of the primary and the secondary, respectively. Then the mass ratio is

$$q_{i,j}=\frac{V_{1,i}V_{1,j}}{V_{2,i}V_{2,j}}.$$

The ratio of the velocity difference between the i-th and the j-th observations is exactly $`q`$. This approach will be most effective with measurements made at such phases that the difference of the velocities is large compared to the uncertainty of their measurements. When more than two pairs of velocities are available, any two pairs of velocities can yield another value of the mass ratio. Note, however, that not all the $`q`$ values are independent (see also Wilson 1941).
2. When the systemic velocity, $`\gamma `$, is available from the SB1 observations, each pair of velocities, obtained at time $`t_i`$, can yield a value of the mass ratio

$$q_i=\frac{V_{1,i}\gamma }{V_{2,i}\gamma }.$$
3. When the elements of the single-lined orbit are available, the expected primary velocity $`V_{1,i}^{exp}`$ at the time of the observation $`t_i`$ can be derived from the elements. In such a case, each secondary velocity $`V_{2,i}`$ can yield a value of the mass ratio

$$q_i=\frac{V_{1,i}^{exp}\gamma }{V_{2,i}\gamma }.$$

In the case in which the zero point of the new velocities is uncertain, and the primary velocity was also measured, $`V_{2,i}`$ can be replaced by $`V_{2,i}+(V_{1,i}^{exp}V_{1,i})`$. This assumes that the new primary and secondary velocities share the same zero point.
4. In the case in which the original single-lined radial velocity measurements are available, one can combine the old velocities with the new double-lined measurements and obtain a new solution for the orbital elements. The new set of elements includes $`K_2`$, the semi-amplitude of the secondary, from which one can derive the mass ratio.

The three binaries considered here have known elements together with the individual measurements (Duquennoy & Mayor, 1991; Latham et al., 2001). We therefore opted to apply the second and the last approaches listed above. The advantage of the second approach is that the spread of the mass ratio values can be examined and evaluated directly, while in the last option this is all hidden in the orbital solution.

## 3 Observations

To search for the spectroscopic signature of the secondaries in our SB1 sample, we observed the target systems and standard stars in 1997 December and 1999 February, May, and November with CSHELL on the IRTF 3-m and in 1999 December with Phoenix at the KPNO 4-m. CSHELL is a high resolution IR spectrometer equipped with a $`256\times 256`$ InSb array that provided a free spectral range of $`c\mathrm{\Delta }\lambda /\lambda 730`$ km/s. We observed with a $`0.5^{\prime \prime }`$ slit, yielding a measured resolution of $`30,000`$. Phoenix, NOAO’s high resolution IR spectrometer, has a $`512\times 1024`$ InSb array which provides a $`1440`$ km/s free spectral range. With Phoenix, we used a $`0.8^{\prime \prime }`$ slit which gave a measured resolution of $`35,000`$.

The spectrometers were centered at $`1.5548\mu `$m for all of our observations. Inspection of atlases of solar photospheric and sunspot spectra showed that this region contains lines suitable for the characterization of stars in the mass range of interest. The region includes atomic lines for identifying the G and early K spectral type stars and molecular lines, most notably OH and H2O, suitable for late K and M spectral type stars (Livingston & Wallace, 1991; Wallace & Livingston, 1992, see also Meyer et al. 1998 for an overview of stellar spectra in the $`H`$ band). This wavelength region is also advantageously free of strong terrestrial absorption lines.

We calculated the dispersion and wavelength calibration of CSHELL and Phoenix at the $`1.5548\mu `$m setting using spectral lamps internal to the instruments. Both spectrometers showed small flexure, at the level of $`3`$ km/s during a night of observations. We corrected for the flexure by reobserving the lamps during the night and also by monitoring the apparent positions of the OH lines detected from the night sky. The spectra were extracted by procedures written by one of us (Prato, 1998).

## 4 Templates

In addition to the target spectra, we also observed a sample of main-sequence standard stars, from spectral type F8 through M6.5, in order to provide templates for the two-dimensional cross-correlation analysis using TODCOR. When possible, we chose a standard whose metallicity would be close to solar, on the expectation that its lines would be deeper than those of older stars. The templates are usually slow rotators, and if needed we generated new calculated templates by broadening the template lines with profiles representing different rotational velocities (see below). We also expected that these standards would provide better templates for our PMS star program.

Table 1 lists the observed standards. The third column provides an estimate of their population type. A numerical entry in this column represents the value of \[Fe/H\] from the catalog of Cayrel de Strobel, et al. (1992). In cases where the catalog listed several determinations, we used the most recent one. A dagger ($``$) in this column is Keenan & McNeil (1989)’s designator for stars they judged to have the best determination of spectral type. OD and YD are Leggett’s (1992) “old disk” and “young disk” designations, based on kinematic and photometric criteria. Column 4 lists the assigned radial velocities of the observed spectra (see below). Column 5 indicates the instrument used, CSHELL (C) or Phoenix (P), on the date indicated in Column 6.

Figure 1 shows all the template spectra measured with the Phoenix spectrometer. The spectra are shifted to the laboratory frame. Some of the prominent lines are indicated following the identifications given by Livingston & Wallace (1991) and Wallace & Livingston (1992). Since the TODCOR analysis correlates spectra of the templates and targets, only the presence of lines matters; their specific identification is unimportant. The sawtooth pattern evident in the spectrum of HR 1543 is an artifact of the illumination of the array for that particular observation. As Table 1 indicates, we observed a similar range of spectral types, and often the same stars, with CSHELL and Phoenix. In Figure 2, we show the CSHELL and Phoenix observations for the M stars only. Comparison shows that most of the features in the complex spectra of the late M types, which appear to be noise, are in fact real.

The radial velocities of the observed spectra of the standards are crucial for our analysis. As a starting point, we assigned to each spectrum the radial velocities of these template stars measured either by the Geneva (Duquennoy & Mayor, 1991) or by the CfA groups (Latham, 1985). However, not all our templates were observed by these two groups. Moreover, some of the observed spectra of the templates were shifted by about 1–2 km/s, apparently the result of random jitter of the spectrographs. To overcome these two difficulties we performed a few iterations to rederive new assigned velocities for each of the observed template spectra. This was done by correlating each observed spectrum with templates of similar classes, using the assigned velocities of the previous iteration. A few iterations led us to assign a velocity for each observed spectrum as given in Table 1.

The assigned velocities are not meant to measure the actual radial velocities of these stars. They are our best estimates for the radial velocity shifts of the observed spectra, and therefore these are the values we used in our analysis with TODCOR. In an ideal world, without any outliers and spectrograph jitters, these velocities would be our best estimate for the actual stellar radial velocities. Actually, we find that for most of the stars for which we had a prior radial velocity measurement (Duquennoy & Mayor, 1991; Latham et al., 2001), our assigned radial velocity was within 1$``$2 km/s of the previous measurement. We therefore estimate that our zero-point is good to within 1.5 km/s. Only one star, GJ 475, displayed a discrepancy of 6 km/s between the CSHELL (8.0 km/s) and Phoenix (1.9 km/s) measurements. Duquennoy & Mayor (1991) derived a velocity of $`6.39\pm 0.08`$ km/s. We, nevertheless, used the two spectra with the assigned velocities, assuming some non-random mechanical or optical shift of 4.5 km/s of the Phoenix instrument had occurred during that exposure.

## 5 The Radial Velocities

To derive the radial velocities of the two stars in each binary, TODCOR requires two specific template spectra, one for the primary and the other for the secondary. We ran TODCOR for a variety of different pairs of templates from our accumulated library to optimize the match between the templates and the observed composite spectra. To each pair of templates we assigned, as a measure of its fit to the binary spectra, the correlation peak obtained by TODCOR for that pair of templates, averaged over all the observed target spectra. We then chose the pair of templates that gave the highest average value for the peak.

To improve the fit between the observed binary spectrum and the two templates, we matched the rotational broadening of the primary and secondary lines with those of the templates. This was done by broadening the template lines with profiles representing different rotational velocities, forming an extended template library. We then chose the broadened templates that best match the observed spectrum of the binary. To extract the radial velocities for the two stars, the brightness ratio of the two templates, $`\alpha `$, was also required. Assuming that each observation of a given binary system had the same brightness ratio, we determined a value to be used for all of its spectra. To choose the best brightness ratio, we ran TODCOR for a grid of selected $`\alpha `$ values, searching again for that value which gave the highest averaged correlation peak. We do not consider the templates and $`\alpha `$ values chosen by this procedure as representative of the true spectral types and light ratios of the binary targets because the template and target stars may differ in abundances and, for PMS stars, surface gravity. We regard the template spectral types and $`\alpha `$ as a set of parameters that yields the best radial velocities. Choosing the best templates and $`\alpha `$ comprised the major part of the first step of our procedure. Only after the best templates and the brightness ratio were determined did we proceed to derive the radial velocities, given in Table 2. Since the essential inputs to TODCOR are the existence and locations of lines in the spectra of the templates and targets, the fundamental outputs are the velocities of the primary and secondary, and hence their mass ratio. The templates and light ratios used are necessary to derive the velocities but are not intended as models of the target binary.

While all observed spectra yielded reliable primary velocities, not all spectra gave reliable secondary velocities. This is because the secondary peak was not always sufficiently prominent, depending on several factors including the SNR of the observed spectrum, the velocity difference between the primary and the secondary, the small differences between the actual spectra of the primary and the secondary on one hand and their templates on the other, and the light ratios between the primary and the secondary. The first two can change between one exposure and another, and can explain why we could derive a reliable secondary velocity from one spectrum and not from the other. For our analysis we chose only spectra that yielded cross correlation functions for which we judged the secondary peak to be high enough above the other spurious peaks of the correlation. Our somewhat stringent criterion led us to accept only 10 secondary velocities for the three stars, all listed in Table 2. Figures 3$``$5 depicts the correlation functions derived by TODCOR for G147-36, G164-47, and HD 144842. The velocity uncertainty for the former two objects is 1.5 km/s and for the latter is 2 km/s, because of its high rotational velocity.

## 6 The Mass Ratios

### 6.1 Derivation from Systemic Velocities

Table 2 shows the results of deriving the mass ratios using only the $`\gamma `$ velocities, drawn from the previous single-lined orbital solutions. The $`\gamma `$ velocity appears as the first entry for each binary. The Julian date, instrument used, component velocities, and derived mass ratio are then listed for each individual observation. The last entry in each section is the average value and uncertainty of the individual mass ratios.

### 6.2 Global Solution for Binary Parameters

In this subsection we present the full double-lined orbital solutions for the three orbits, derived with our ORB code (Mazeh et al., 1993, 2001). To derive the double-lined solution, we assign different uncertainties to the primary and secondary velocities of the different measurements. We estimated that the uncertainty of the previous velocities of G147-36 and G164-67 (Latham et al., 2001) are 0.5 km/s, and those of HD 144284 (Duquennoy & Mayor, 1991) are 1 km/s. The latter somewhat bigger uncertainty is attributable to the large rotational width of the stellar lines. The scatter we find for these velocities (see Table 3) is very close to these estimates.

The uncertainty of the secondary velocities depends strongly on the free spectral range, the number and depth of the absorption lines, and, mainly, on the light ratio between the primary and the secondary. Using two different sets of data, obtained with CSHELL and Phoenix, complicates the analysis. Therefore, our uncertainty estimates are poorly constrained. We estimated that the uncertainty of our secondary velocities for G147-36 and G164-67 are 1.5 km/s and those of HD 144284 are 2 km/s, again because of the large rotational velocity of the latter.

Table 3 presents our results, giving the derived values of the period, P, in days, eccentricity, e, longitude of the periastron, $`\omega `$, in degrees, time of periastron passage T, in J.D.$``$2440000, primary velocity amplitude, $`K_1`$, in km/s, and the center-of-mass velocity, $`\gamma `$, in km/s. The only new parameter coming from the secondary radial velocities is the velocity amplitude of the secondary, $`K_2`$, in km/s, from which we derived the mass ratio $`q`$. Finally, Table 3 lists $`\sigma _{1,old}`$, $`\sigma _{1,new}`$ and $`\sigma _{2,new}`$, in km/s — the r.m.s. of the residuals of the old and new velocities of the primary, and the new secondary velocities. The orbital elements derived here are entirely consistent with the original published values. The orbits are plotted in Figures 6$``$8.

The scatter for the primary and the secondary velocities, listed in Table 3, gives an external estimate of our precision. As can be seen from Table 3, these estimates are based on extremely small numbers, only two measurements in the case of G164-67. It is therefore not surprising that the secondary radial-velocity scatter ranges from 4.0 to 0.23 km/s. Our estimates for the uncertainty of the secondary velocities are between these two values.

Table 4 is a summary of the mass ratios and masses estimated for the three binaries. The primary masses for G147-36 and G164-67 are from Carney et al. (1994) and for HD 144284 from Mazeh et al. (1997). The latter is based only on the spectral type of the system. Arbitrarily, we assigned an uncertainty of 0.1 $`M_{\mathrm{}}`$ to these primary mass estimates. The table also lists the minimum mass ratios derived from the SB1 orbits. Statistically, the derived mass ratios are expected to be close to the minimum values, but never smaller. This is indeed the case, within the $`1\sigma `$ range, for all three binaries.

## 7 Discussion

The work presented here demonstrates the potential of IR spectroscopy for detecting faint low-mass secondaries in single-lined spectroscopic binaries and deriving their mass ratios. We have shown that with the CSHELL and the Phoenix spectrographs we can reach a mass ratio of $`0.38\pm 0.01`$. To date, this is among the smallest known mass ratios for a main-sequence spectroscopic binary. Mazeh et al. (1997), observing in the red part of the visual spectrum, derived mass ratios of $`0.57\pm 0.02`$ and $`0.48\pm 0.03`$ for two binaries.

The detection of a secondary spectrum and the derivation of its radial velocity is hindered by the presence of the bright primary in the combined spectrum. For our purposes, the bright spectrum of the primary is a strong background source that has to be removed from the combined spectrum before we can measure the secondary radial velocity. The subtraction of the primary leaves behind noise, which ideally would be determined by Poisson statistics, but in reality is determined by the imperfect match between the primary spectrum and its template. To reach the faint stellar secondaries, we have to exploit every possible advantage available to improve the SNR of the secondary spectrum.

In this work we have used two such advantages. We observed the binaries in the IR regime, where the light ratio of the primary/secondary is much smaller than in the visual, and we used the two-dimensional cross-correlation analysis of TODCOR to extract the maximum amount of spectral information from the observed spectrum. To improve on this further, we need high sensitivity IR spectra covering as large a spectral range as possible. Fortunately, IR echelle spectrometers providing several orders simultaneously are now coming into operation at the large aperture telescopes. We have begun systematic surveys to detect low-mass secondaries in both main-sequence and PMS binaries (Prato et al., 2001b).

We thank W. Vacca (IRTF), R. Joyce (NOAO) and K. Hinkle (NOAO) for thorough support of CSHELL and Phoenix and for advice and help in their use. We thank J. Carr for advice about spectral type standards and for sharing with us his list of standards. We thank the referee for an exceptionally informative and thorough report. Our research made extensive use of the SIMBAD database, operated at CDS, Strasbourg, France. This research was supported in part by US-Israel Binational Science Foundation Grant No. 97-00460 and by the Israel Science Foundation, grant no. 40/00 (to T.M.), and by the US NSF, Grant 98-19694 (to M.S.).


format: text

text:
# 1 Introduction

## 1 Introduction

It is widely believed that the role of spin is decreasing with the rise of beam energy and transverse momentum ($`p_T`$) of secondary particles. In particular, the perturbative QCD predicts the vanishing of spin effects with the rise of $`p_T`$ . Experiments performed during last years have shown much more complicated picture. It was found that polarization of $`\mathrm{\Lambda }`$-hyperons produced in hadron-hadron collisions is almost energy independent and is rising with $`p_T`$ and $`x_F=2p_Z^{}/\sqrt{s}`$, where $`p_Z^{}`$ and $`\sqrt{s}`$ are, respectively, longitudinal momentum of secondary particle and total reaction energy in C.M. . Single-spin asymmetry ($`A_N`$) is rising with $`p_T`$ and $`x_F`$ at least for $`\pi `$\- mesons and it differs from zero up to the highest available beam energies (200 GeV) .

It seems, that the most promising among theoretical approaches to the solution of large single-spin asymmetry problem, observed in experiments, is to take into account higher twists contributions (see, for example ).

The existing data on single-spin asymmetries have in most cases, limited statistical accuracy, $`x_T=2p_T/\sqrt{s}`$ and $`x_F`$ ranges as well as incomplete particle identification. That does not allow one to make definite conclusions about asymmetry dependence on kinematical variables and quantum numbers of hadrons.

The goal of the experiment is to study single-spin asymmetry ($`A_N`$) in reactions:

$$p+p(A)h+X,$$
(1)
where $`h`$ denotes $`\pi ^\pm `$, $`K^\pm `$, p, p̄ or hadron pair .

It includes, in particular, the following points:

* $`A_N`$ vs $`x_T`$ and $`x_F`$ ($`x_T`$ and $`x_F`$ scaling check);
* $`A_N`$ vs quark flavor (u,d,s);
* $`A_N`$ vs A - atomic number (nature of A-dependence);
* $`A_N`$ vs $`x_T`$ for hadron pairs.

The experiment has been carried out with FODS-2 double arm spectrometer designed to study high $`p_T`$ single and pair hadron production in pp-, $`\pi `$p-, pA- and $`\pi `$A- collisions . A new IHEP polarized 40 GeV/c proton beam gives an additional tool to study high $`p_T`$ processes in reactions (1). The beam energy and available intensity of polarized beam allow one to measure $`A_N`$ in a wide $`x_F`$ and $`x_T`$ range ($`x_T<0.9`$) for not so high $`p_T`$ (up to 3.4 GeV/c) of charged hadrons.

Preliminary results of this experiment, based on half the statistics, were reported at SPIN-95 workshop in Protvino .

## 2 Beam

Polarized protons  are produced via decays of $`\mathrm{\Lambda }`$-hyperons. A slow extracted 70 GeV/c proton beam with intensity up to $`10^{13}`$ ppp strikes the primary target. Charged particles, produced in the primary target, are deflected and absorbed in a beam dump. Protons from $`\mathrm{\Lambda }`$-decays are collimated in vertical direction to select polarization sign and beam intensity. Average polarization of 40 GeV/c beam (for $`\mathrm{\Delta }p/p=\pm 4.5\%`$) is $`39_3^{+1}\%`$ and it was estimated using Monte Carlo simulation of the polarized beam production. The vertical beam position on the secondary target is corrected by two magnets. The beam polarization sign changes automatically each 15 minutes during 3-4 spills. Beam monitoring is performed by ionization chambers and scintillator monitors (arm triggers) with relative precision better than 1%. Calculated and measured beam parameters are in a good agreement, which gives us confidence in the knowledge of its polarization calculation.

## 3 Spectrometer

The experimental setup FODS-2 is a rotating double arm spectrometer  shown schematically in Fig.1. It includes a two gap magnet, drift chambers to measure momentum and angles of particles production, two Cerenkov ring imaging detectors (SCOCH) to identify particles, scintillation counters, hodoscopes and hadronic calorimeters to make trigger decision. The SCOCH detector is able to identify hadrons in the following momentum ($`p`$) ranges: $`\pi ^\pm `$ ($`2<p<24`$ GeV/c), $`K^\pm `$ ($`4<p<24`$ GeV/c), p and p̄ from 8 up to 24 GeV/c. Protons and antiprotons were also identified in the momentum range from 5 to 8 GeV/c using the absence of a signal in the SCOCH detector. The transverse position ($`X,Y`$) of an interaction point in the hydrogen target ($`D`$=7 cm, $`L`$=40 cm, $`L/\lambda _{abs}`$=5%) is measured with beam hodoscopes. Longitudinal coordinate ($`Z`$) is determined as a result of matching $`X,Y`$-coordinates, measured with beam hodoscopes, with a particle trajectory, reconstructed using drift chambers downstream of the magnet. Hodoscopes resolution is $`\pm `$1 mm, both for $`X`$ and $`Y`$ coordinates. The total number of channels in each hodoscope plane is 32, and average multiplicity was around 2 per plane. When the multiplicity is not equal to unity a hodoscope channel closest to the center of gravity for the beam during the spill is used. Ionization chambers were used to determine the center of gravity for the beam.

## 4 Measurements

The measurements were carried out during 6 days of the November Run in 1994. The mean intensity of the polarized beam was $`910^6`$ ppp, limited by radiation environment. Now the shielding has been improved and nominal intensity up to $`2.610^7`$ can be achieved. The trigger for each arm required signals from the scintillation counters and energy deposition above threshold in the hadronic calorimeter. The total number of polarized protons that striked the target was about $`1.810^{11}`$ for each polarization sign. The corresponding number of recorded events is $`6.310^5`$/polarization/arm. 24% of them have been reconstructed and identified for the final analysis.

## 5 Data reduction

The reconstructed trajectory of a particle downstream the spectrometer magnet is used to determine its momentum and production angles. Since the magnet deflects charged particles in vertical direction it does not change essentially polar angle ($`\theta `$), but changes azimuthal angle ($`\varphi `$) depending on momentum and sign of the particle. The mean value of $`cos\varphi `$ varies from 0.80 to 0.89.

Particle momentum was determined by matching incident particle transverse coordinates in the target ($`X,Y`$), measured with the beam hodoscopes having the trajectory downstream the magnet. The momentum ($`p`$) of hadron is used to determine its mass ($`m`$):

$$m^2=p^2(1\beta ^2)/\beta ^2,$$
(2)
where $`\beta `$ is hadron velocity, measured with the SCOCH detector . Mass spectrum ($`m^2`$) of hadrons measured with the SCOCH is shown in Fig.2. For heavy hadrons (protons) the main error in $`m^2`$ is due to the uncertainty in the momentum. For light particles ($`\pi `$-meson) the error in $`m^2`$ is mainly due to the uncertainty in $`\beta `$.

## 6 Results

Single-spin asymmetry has been calculated for two arms ($`L`$-left, $`R`$-right):

$$A_N^L=\frac{1}{P_Bcos\varphi }\frac{N_L^uN_L^d}{N_L^u+N_L^d},$$
(3)
$$A_N^R=\frac{1}{P_Bcos\varphi }\frac{N_R^uN_R^d}{N_R^u+N_R^d},$$
(4)
where $`P_B`$ is the average beam polarization (39%), and $`N_{L(R)}^{u(d)}`$ is the number of normalized events, for protons, polarized up ($`u`$) or down ($`d`$) for the left ($`L`$) or right ($`R`$) arm. The so called canonical helicity frame is used here for asymmetry definition, as it is in . Asymmetries for two arms have been averaged, which decreased some systematic errors and doubled the statistics. Only statistical errors are taken into account in the results. The relative systematic error in $`A_N`$ due to possible errors in beam polarization, beam monitoring and $`cos\varphi `$ measurement is estimated to be less than 20%.

Along with the asymmetries for single particles the asymmetries for particle ratios have been calculated as well since they are free from beam monitoring problems.

Asymmetries for $`\pi ^\pm `$, $`K^\pm `$, protons and antiprotons are presented as functions of $`p_T`$ in Table 12 and 3, respectively.

Asymmetry for $`\pi ^+`$ is shown in Fig.3, as a function of $`x_T`$. The lines in the figure represent fit by expression:

$$A_N=A_0(x_TX_0),$$
(5)
where slope ($`A_0`$) and sign change point ($`X_0`$) are free parameters.

The bin width for $`p_T`$ is 0.25 GeV/c and mean $`p_T`$ varies from 0.66 to 3.37 GeV/c. The plotted values of $`x_T`$ correspond to the mean values within the bins. Since the laboratory angle of the FODS-2 spectrometer is $`9^o`$, mean values of corresponding $`x_F`$ at 40 GeV beam energy are not exactly zero and are rising from 0.02 to 0.10 with the rise of $`p_T`$. Asymmetry for $`\pi ^+`$ is rising with $`x_T`$ and changes sign near $`x_T=0.37`$.

The data for the BNL experiment  are also shown in Fig.3 for comparison. It is seen from Fig.3 that $`A_N`$ changes sign at the same $`x_T`$ for three different beam energies (13.3, 18.5 and 40 GeV). The cut $`x_F0.18`$ is used for the BNL data to compare with this experiment.

Asymmetry for $`\pi ^{}`$ is shown in Fig.4, where BNL data are also shown for comparison. Asymmetry for $`\pi ^{}`$ is negative and it is rising in absolute value with $`x_T`$ increase. $`A_N`$ at the BNL energies is close to zero with some indication for positive slope $`A_0`$.

Asymmetry for $`K^+`$ (Fig.5) is rising with $`x_T`$ similar to $`\pi ^+`$ data and changes sign at $`x_T`$ near 0.37.

Asymmetries for $`K^{}`$ and antiprotons (Fig.5) have an indication of the rise with $`x_T`$ similar to $`K^+`$ data but statistical errors are too high to make this conclusion final. The last three points for p̄ from Table 3 are combined in Fig.5 to decrese error.

The comparison of $`A_N`$ for protons with the BNL data are shown in Fig.6, as a function of $`x_T`$. For both experiments $`A_N`$ does not depend on $`x_T`$ and at 40 GeV the mean value of $`A_N`$ is $`0.050\pm 0.009`$. The mean value of $`A_N`$ within explored range of $`x_T`$ is shown in Fig.7 for three energies as a function of $`\sqrt{s}`$. It is seen from Fig.7, that the mean value of $`A_N`$ is well described by a linear function of $`\sqrt{s}`$:

$$A_N=(0.0113\pm 0.0027)(4.34\pm 0.45\sqrt{s}/GeV).$$
(6)
This result, if it is confirmed by other experiments, can be used to measure beam polarization in a wide range of energies.

As was mentioned above asymmetry for particle ratios, which under the first approximation equals the difference of asymmetries for two types of particles, is not sensitive to beam monitoring and thus has less systematics. Predictions of models in many cases are also better defined for particle ratios than for each particle separately. As an example of asymmetries for particle ratios Fig.8 shows $`A_N`$ for $`K^+/\pi ^+`$ and $`K^{}/\pi ^{}`$ ratios. For $`K^+/\pi ^+`$-ratio $`A_N`$ is consistent with zero, as could be expected, since the fragmentation of scattered valence $`u`$-quark is a main production source of both mesons . Asymmetry for $`K^{}/\pi ^{}`$ is increasing with $`x_T`$ and it is positive for $`x_T`$ 0.28.

Fit parameters for the data, shown in Figs. 3-8, are presented in Table 4. As is seen from Table 4 there is a significant slope ($`A_0`$) in dependence on $`A_N`$ vs $`x_T`$ (5) for all particles except protons. The sign of the slope is positive for particles containing valence $`u`$-quark of colliding hadrons (for $`\pi ^+,K^+`$), or sea quark (for $`K^{}`$,p̄). For particles containing valence $`d`$-quark from colliding hadrons (for $`\pi ^{}`$) it is negative.

As was shown above, the asymmetry dependence on $`x_T`$ can be described by linear function (5). To study the dependence of parameters $`A_0`$ and $`X_0`$ of (5) on energy we fitted $`x_T`$ dependence of asymmetry for several experiments in a wide range of energies . Since in some experiments only $`\pi ^0`$-mesons have been detected, their asymmetry is compared below with mean $`A_N`$ for $`\pi ^+`$ and $`\pi ^{}`$-mesons, averaged according to a parton model with weights proportional to their production cross sections :

$$A_N(\pi ^o)=\frac{A_N(\pi ^+)R(x_T)+A_N(\pi ^{})}{R(x_T)+1},$$
(7)
where $`R(x_T)`$ is a ratio of the production cross sections for $`\pi ^+`$ and $`\pi ^{}`$ -mesons. The resulting dependences of $`X_0`$ and $`A_0`$ on $`\sqrt{s}`$ are shown in Figs.9-10, respectively. The dependence of $`X_0`$ on $`\sqrt{s}`$ is consistent with constant value of $`X_0=0.374\pm 0.009`$, while $`A_0`$ decreases with the rise of energy:

$$A_0=(0.029\pm 0.005)(21\pm 2\sqrt{s}/GeV),$$
(8)
which means the decrease of spin effects with the rise of energy. Scaling for $`X_0`$, observed in different experiments, proves that we have no significant systematics for asymmetry in this experiment. One point in Fig.10, corresponding to the PROZA-M experiment result is much higher than other points, which could be related to different beam ($`\pi ^{}`$) for that experiment . Of course, the exact dependence of $`A_0`$ on energy is unavailable and the linear expression is just the simplest one.

It is interesting to compare the observed scaling for parameter $`X_0`$ with predictions of theoretical models. Predictions for asymmetry in the central region for a wide range of energies have been recently made in the framework of QCD . The main role in this model belongs to the orbital angular momentum of the quark-antiquark cloud in the internal structure of constituent quarks. For polarization of the constituent quarks SU(6)-model values $`P_U`$ = 2/3 and $`P_D`$ = -1/3 are used. Predictions of asymmetry for $`\pi ^+`$ and $`\pi ^{}`$-mesons have been made for 70, 200 and 800 GeV. These predictions along with their extrapolation to 40 GeV are shown in Fig.11 as the functions of $`x_T`$. It is seen from Fig.11, that the value of $`x_T`$, where $`A_N=0`$ is decreasing with the energy rise (which corresponds to constant $`p_T=0.75`$ GeV/c), contrary to the experimental data where it takes place at constant $`x_T`$ around 0.37.

The comparison of data and model predictions for 40 GeV is shown in Fig.12. The 40 GeV energy model predictions are in a qualitative agreement with the results of this experiment. The model predicts the rise of $`A_N`$ with $`x_T`$, the opposite sign of $`A_N`$ for $`\pi ^+`$ and $`\pi ^{}`$, and half as much the absolute value of asymmetry for $`\pi ^{}`$ against for $`\pi ^+`$. The model predicts the increase for $`A_0`$ and the decrease for $`X_0`$ with the rise of energy. That is in some contradiction with the results shown in Figs.9-10 for $`\pi ^0`$-mesons. Asymmetry for $`\pi ^0`$-meson at 200 GeV is very small , which is possible to explain if $`P_U=P_D`$ and hence $`A_N^{\pi +}`$ = -$`A_N^\pi `$, which is not seen at 40 GeV energy and below it.

## 7 Conclusion

The first asymmetry measurements have been performed for $`\pi ^\pm `$, $`K^\pm `$, p and p̄ with 40 GeV/c IHEP polarized proton beam. Asymmetries for $`K^\pm `$ and p̄ have never been measured before. Asymmetries for $`\pi ^\pm `$, $`K^\pm `$, p̄ show approximately linear dependence on $`x_T`$ and change the sign near 0.37. The slope sign of this dependence is negative for $`\pi ^{}`$ and positive for other hadrons, except protons. For protons the negative asymmetry, independent of $`x_T`$ has been found, which grows with energy in absolute value as the comparison with the other experiment indicates.

The SU(6) model predicts asymmetry that is in a qualitative agreement with our data for $`\pi `$-mesons.

A significant increase in statistics is possible (by a factor around 40) in future runs which will allow us to make definite conclusion about asymmetries in a wide range of $`x_T`$ and $`x_F`$.

We gratefully acknowledge the assistance and support of the IHEP staff and directorate. The research was supported in part by the Russian Fund for Fundamental Research.


text:
# H𝛼 Distance Constraints for High Velocity Clouds in the Galactic Halo

## 1. Introduction

Observations of the Galactic halo make a compelling case that the formation of halos continues to the present day (Wyse 1999). The halo appears to have built up through a process of accretion and merging of low-mass structures which is still going on at a low level. Hierarchical cold dark matter (CDM) simulations, however, predict that the Galactic halo should have many more satellites than are actually observed (Klypin et al. 1999; Moore et al.1999).

As much as 40% of the sky is peppered with high-velocity H$`\mathrm{I}`$ clouds (HVCs) which do not conform to orderly Galactic rotation (Wakker & van Woerden 1991; Putman 2000). These are interesting accretion candidates $``$ particularly if they are associated with dark matter ‘mini halos’ $``$ except that their distances, $`d`$, are unknown for all but a few sources. As a result, fundamental physical quantities $``$ size ($`d`$) and mass ($`d^2`$) $``$ are unconstrained which has encouraged wide speculation as to the nature of HVCs (Wakker & van Woerden 1997).

Indeed, the current renaissance in HVC studies can be traced in part to an interesting suggestion by Blitz et al. (1999). They showed that the velocity centroids and groupings of positive/negative velocity clouds on the sky may be understood within a reference frame centered on the Local Group barycenter (cf. Zwaan & Briggs 2000). They interpret HVCs as gas clouds accreting onto the Local Group over a megaparsec sphere. Braun & Burton (2000) identify specific examples of compact clouds that have ‘rotation curves’ consistent with CDM mass profiles. For sources at 700 kpc, the kinematic signatures imply a high dark-to-visible mass ratio of 10$``$50.

This paper addresses the distance problem. In earlier papers (e.g.,  Bland-Hawthorn et al. 1998; Bland-Hawthorn & Maloney 1999 $``$ BM99), we showed that faint H$`\alpha `$ measures from distant H$`\mathrm{I}`$ clouds could, in principle, be used to estimate crude distances to the clouds. The H$`\alpha `$ emission measure from any cloud which can be detected at 21 cm is a direct measure of the Lyc (Lyman continuum) radiation field, independent of distance. In order to interpret the H$`\alpha `$ emission, we require a realistic model of the Galactic halo ionizing field. Here, we present some recent developments in deriving H$`\alpha `$ distances which significantly extends our earlier work.

## 2. The Escape of UV radiation from the Galaxy

Since we have discussed the UV escape fraction ($`\widehat{f}_{\mathrm{esc}}`$) from the Galaxy at length in earlier articles (q.v.,  Bland-Hawthorn & Putman 2001), only a short discussion is given here. The likely value of $`\widehat{f}_{\mathrm{esc}}`$ is still uncertain but there is increasing evidence that UV does radiate far from HII regions. To cite one of several recent demonstrations, the warm (diffuse) ionized medium in external galaxies clearly shows the hallmark of the HII regions even after these have been masked out of the H$`\alpha `$ image (e.g.,  Zurita et al. 2001; Cianci 2001).

Further evidence that UV must escape the Galaxy comes from the measured electron density profile from halo pulsars. Manchester & Taylor (2000; see also Nordgren, Cordes & Terzian 1992) have modelled this with a scale height of 800 pc which exceeds or is comparable to the scale height of the diffuse H$`\mathrm{I}`$ (warm neutral medium; Lockman 1984). Without fine-tuning, it is unlikely that the Reynolds layer represents a radiation-bounded medium within a co-extensive H$`\mathrm{I}`$ envelope. We know that the radiation field must be soft from the weakness of HeI$`\lambda `$5876 and non-detection of \[HeII\]$`\lambda `$4686 (Reynolds & Tufte 1995). Furthermore, the observed weakness of \[OI\]$`\lambda `$6300/H$`\alpha `$ indicates two things: (i) the ionization fraction must be high (Reynolds 1989), (ii) all of the UV photons produced in the disk cannot be absorbed in radiation-bounded H$`\mathrm{II}`$ regions (Domgörgen & Mathis 1994).

It is important to appreciate that the Fabry-Perot ‘staring’ method is so sensitive that we only need $`f_{\mathrm{esc}}`$ $``$ 1% for the H$`\alpha `$ distance method to be useful. In fact, our preliminary analysis of H$`\alpha `$ levels from H$`\mathrm{I}`$ clouds with stellar distance brackets suggests that indeed $`f_{\mathrm{esc}}`$ $``$ 1$``$2% (see § 6). At distances of 300 kpc along the polar axis, the Galaxy field is comparable to the cosmic UV field. The expected H$`\alpha `$ levels are of order $`12`$ mR and it may be possible to reach these levels2221 milliRayleigh $`=`$ $`10^3/4\pi `$ phot cm-2 s-1 sr-1 $`=`$ $`2.4\times 10^{10}`$ erg cm-2 s-1 sr-1 at H$`\alpha `$. with new differential techniques (e.g.,  Glazebrook & Bland-Hawthorn 2001).

## 3. The importance of spiral arms

Our early attempts to derive H$`\alpha `$ distances made use of a smooth exponential distribution of ionizing sources (BM99). For HVC distances within 10 kpc of the plane, we need a more realistic distribution of ionizing sources than the exponential disk model. This would be straightforward if we knew the exact location of all O stars, and the precise dust distribution throughout the Galaxy. But this will not be possible until the GAIA astrometric mission flies in 2010. However, most studies of spirality in the Galaxy agree that the tangent points of the spiral arms are well defined over a wide range of methods and wavelengths, in particular, the distribution of pulsar dispersion measures with galactic longitude (Taylor & Cordes 1993).

For this reason, the non-axisymmetric component of the diskhalo ionization model links its fortunes to the standard model for determining pulsar distances. Rough distances to pulsars are deduced from the dispersion (and scattering) measure due to warm electrons along the line of sight. Early attempts used a smooth distribution of electrons (e.g., Manchester & Taylor 1981) although Lyne, Manchester & Taylor (1985) showed that typical distance estimates have random errors as large as a factor of two. After the inclusion of smooth spiral arms, Taylor & Cordes (1993) predict that most distances should be good to $``$20%. This level of accuracy is somewhat surprising when one examines face-on spirals in the Ultraviolet Imaging Telescope (UIT) database. But the distance model is largely borne out by lower limits derived from pulsar sight lines which show H$`\mathrm{I}`$ in absorption.

## 4. Stellar bulge contribution

A component largely overlooked to date is the UV field arising from the Galactic spheroid (O’Connell 1999). A ‘UV upturn’ component is now well established for ellipticals and S0s (Greggio & Renzini 1999; Macchetto et al. 1996; Bica et al. 1996). Binette et al. (1994) suggest this may arise from post-AGB stars evolving from an old stellar population. A rough estimate is $`7\times 10^{50}(M_{\mathrm{bulge}}/10^{10}M_{})`$ ionizing phot s-1. $`M_{\mathrm{bulge}}`$ is the mass of the Galactic halo bulge which is uncertain to a factor of two, i.e.,  $`M_{\mathrm{bulge}}12\times 10^{10}M_{}`$ (Dwek et al. 1995; Zhao 1996). Since much of the bulge stars lie outside most of the absorbing ISM, the UV flux which escapes into the halo may be non-negligible and more smoothly distributed compared to flux which escapes the disk. When a bulge component is discussed below, we include an isothermal distribution of sources with a total luminosity of $`7\times 10^{50}`$ phot s-1. We assume that all of the UV escapes, although the dusty disk blocks half of the radiation. The hot corona component is discussed in BM99. Both sources of UV are examined briefly in the next section.

## 5. The disk-halo ionization model

The diskhalo model includes five basic components: the spiral arm (or exponential) disk defined by the OB star population; stellar bulge; hot coronal halo; LMC (and SMC); cosmic background. In addition, the resulting radiation field can be moderated by the presence of a plane parallel opacity law, and projected H$`\alpha `$ emission measures can be derived with a variety of schemes.

In Fig. 1, we show how the spiral arms are represented in the ionization model. In the near field (r $`<`$ 10 kpc), the distribution is essentially identical to the Taylor-Cordes model. However, the spiral arms are incomplete on the far side of the Galaxy which required us to extend the arm coverage by 60%. This was done by splicing the Taylor-Cordes model onto the basic model parameters of Ortiz & Lepine (1993). The spiral arms cover a total area of 100 kpc2 within a circle with 12 kpc radius, compared with the exponential model which extends out to a radius of 25 kpc (although most of the emissivity lies within a 12 kpc radius).

The diskhalo manual (Bland-Hawthorn & Maloney 2001) describes how the electron density field is inverted to produce the surface density of ionizing photons. Fig. 2 shows how the three different disk distributions compare as we approach the disk along the polar axis of the Galaxy. Fig. 3 presents a cross section through the halo field and shows that the form of the dusty spiral is very different to that of the exponential disk within 10 kpc, which can also be seen in Fig. 2.

An important prediction of the spiral arm and the exponential models is that ionizing flux seen by a cloud ‘inside’ the Solar Circle should be much higher than the flux seen by a cloud on the ‘outside’. We have tried to illustrate this effect in Fig. 4. For this trend to be evident, a large nearby sample of clouds is required over a wide solid angle. One possible population which might reveal this effect are the intermediate velocity clouds (IVCs). In the next section, we show the prediction for one such cloud, Complex K.

Fig. 4 shows how a stellar UV bulge (or, for that matter, a hot UV corona) can wash out some of the contrast produced by spiral arms (see also Fig. 8). The presence of a UV-bright stellar bulge can remove the near-field solution which typically arises with the spiral arm models. This phenomenon is seen in the exponential disk predictions in Fig. 5. The hot galactic corona is not as effective as the stellar bulge at washing out the spiral structure since, even though it has a scale length a factor of 5 larger than the bulge, the expected flux levels are lower by about the same factor.

## 6. Predicted H$`\alpha `$ distances to H$`\mathrm{I}`$ clouds

Here we present predictions for HVCs with measured H$`\alpha `$ fluxes. We emphasize that the method is only intended as a statistical constraint. Even with large numbers of H$`\alpha `$ detections for each cloud complex, it may only be useful to within a factor of a few and even this level of confidence needs to be tested.

We used the diskhalo code with the dusty spiral disk ($`\tau _{LL}`$$`=`$2.8). We use a conservative model with a hot corona but no stellar bulge component. The choice of $`\tau _{LL}`$ is discussed in detail elsewhere but largely arises from the WHAM detections of Complexes A, C, K and M (see below). Weiner et al. (2001, these proceedings) derive a similar value from an analysis of their own observations. Since galactic coordinates have not yet been published for some of the detections, we have averaged our calculations over the observed H$`\mathrm{I}`$.

A major benefit of the spiral arm model is the order of magnitude larger contrast in H$`\alpha `$ that one obtains compared to the exponential disk model. A good example of this comes from the Complex L sight line presented in Fig. 5. This cloud is found to be very bright in H$`\alpha `$ by both the Las Campanas (Weiner et al. 2001) and TAURUS teams (Putman et al. 2001, in prep.). The exponential model fails by a large factor to explain the signal whereas this is a natural consequence of the spiral arm model if the complex is 1 kpc from the disk, in which case the cloud system lies directly above a spiral arm. Independent support for this comes from the observation by both teams that the \[NII\]/H$`\alpha `$ ratios are very elevated in this system, an effect seen at low latitudes in external edge-on galaxies. The same temperature effect (Reynolds, Haffner & Tufte 1999) is seen in the Smith Cloud which is also thought to lie within a few kpc of the disk (see below).

But the contrast effect comes at a price. For most of the modelled sight lines in Fig. 5, there is a near-field and a far-field prediction. This is because in almost any direction away from the Sun, the radius vector crosses a spiral arm. Multiple solutions are generally avoided at latitudes higher than 20 because as we move away from the disk, the details of the distribution become less important. An example is the Smith Cloud (GCP) discussed in Bland-Hawthorn et al. (1998). This cloud is predicted to lie at 1 kpc or roughly 15 kpc. Two of the sight lines (Complex L, Co-Rotate) are multiple valued for the observed H$`\alpha `$.

The WHAM detections of Complexes A, M and C are also modelled in Fig. 5. These clouds are particularly important as they have distance bounds from the stellar absorption line (SAL) method. If a star of known distance lies beyond the cloud, and another star of known distance falls in front, we obtain a distance bracket for the cloud if it is seen in absorption against the distant star, but not in the nearer star. Complex M (specifically, Cloud MII) has an upper limit of 4 kpc on its distance (Ryans et al. 1997); no lower bound exists. The far-field limit is only consistent at the 1.3$`\sigma `$ level. Complex A lies between 4 and 10 kpc from the Sun (van Woerden et al. 1999) which appears consistent with the far field limit. There is a suggestion that this range could be tightened to 8 and 10 kpc should the complex not be detected on the spectrum of PG 0832+675, in which case the far field prediction is a factor of two too small. Complex C has a secure lower limit on its distance of 1 kpc, based upon five stellar probes; a much weaker limit of $`>`$6 kpc is provided by the non-detection of the cloud in BS 16034-0114 (Wakker 2001). Both near- and far-field solutions are consistent with a firm distance limit ($`>`$1 kpc); only the far-field solution is consistent with the weaker limit ($`>`$6 kpc).

The WHAM team have recently published H$`\alpha `$ observations for the IVC Complex K (Haffner et al. 2001). There is little absorption line data in this direction although the cloud does appear to lie between 0.3 and 7.7 kpc. The sight line to Complex K lies directly over the tangent point of a spiral arm which explains the broad distribution in Fig. 5. The WHAM team present three averaged H$`\alpha `$ measurements for this cloud. In the absence of a characteristic H$`\alpha `$ measurement for a fully mapped cloud, we take the peak H$`\alpha `$ detection towards the center of a cloud as the most important for distance determination.333This emphasizes a weakness in the H$`\alpha `$ method $``$ see Bland-Hawthorn & Maloney 1999$`b`$ for a detailed discussion $``$ in that a detailed mapping is normally required to reliably interpret what part of the cloud is being lit up in H$`\alpha `$. However, this is not always a problem. Weiner et al. (2001, these proceedings) and Tufte et al. (1998) find that the internal H$`\alpha `$ dispersion in a significant number of HVCs can be small. The brighter detections are entirely consistent with the SAL bounds, although the weaker detections would push the far-field limit out to 10 kpc or more.

In summary, for HVCs with well defined distance bounds, we find (see also Weiner et al. 2001) that the observed H$`\alpha `$ is roughly consistent with $`\widehat{f}_{\mathrm{esc}}`$ $``$ 6% ($`\tau _{LL}`$$`=`$ 2.8) although the present uncertainties are about a factor of 2. Note that this escape value $`\widehat{f}_{\mathrm{esc}}`$ is defined orthogonally to the disk plane; the solid-angle averaged value is $`f_{\mathrm{esc}}`$ $``$ 1$``$2%.

Finally, the Las Campanas team (Weiner et al. 2001) have achieved some faint detections to clouds which are thought to be farther afield. We show the predictions for three of these (Co-Rotate, GCN and Population N). We have run the code for all HVCs with H$`\alpha `$ detections to date. Population N appears to lie at the greatest distance (30$``$45 kpc) if the H$`\alpha `$ emission arises from the disk radiation field.

## 7. The Magellanic System

The LMC has several highly active star forming regions, particularly regions of very recent star formation (Shapley III) and of ongoing star formation (30 Doradus). The basic ionizing requirement of the LMC from combined UV, optical and radio studies appears to be $`5\times 10^{51}\mathrm{phot}\mathrm{cm}^2\mathrm{s}^1`$. Within a factor of two, this is consistent with OB star counts (Walborn 1984; Parker 1993), radio continuum observations (McGee, Brooks & Batchelor 1972; Israel & Koornneef 1979), and vacuum ultraviolet observations (Smith et al. 1987) of the LMC. However, the total number of ionizing photons produced by the LMC H$`\mathrm{II}`$ regions, spread over a 5 kpc region, may be as high as $`1.53\times 10^{52}\mathrm{phot}\mathrm{s}^1`$ (Smith et al. 1987). OB star counts around 30 Dor (e.g., Parker 1993) could well underestimate the total ionizing flux by a substantial factor. Kennicutt et al. (1995) suggest that fully one third of the ionizing radiation in the LMC arises from within 0.5 of 30 Dor. The ground-based results may suffer from crowding which means that the total number of stars is underestimated.

We now test whether the Magellanic H$`\mathrm{I}`$ Bridge can constrain the escape fraction of UV photons from the Magellanic Clouds. Fig. 6 shows one such model prediction for the Bridge. Fujimoto & Sofue (1976) give specific positions for the LMC and SMC in Galactic coordinates. Here, we take the total ionizing flux from the LMC to be $`1\times 10^{52}\mathrm{phot}\mathrm{s}^1`$ and the SMC to be an order of magnitude smaller. For both galaxies, we assume that 15% of the UV escapes. The expected level in the Bridge is of order $`50100`$ mR, an order of magnitude higher than levels produced by the Galaxy at that large polar angle.

There are two published claims of H$`\alpha `$ in the Magellanic Bridge. Johnson, Meaburn & Osman (1982) claim to see diffuse H$`\alpha `$ across the entire HI bridge extending over 30 or more at an observed surface brightness of 8 R$`\pm `$4 R. Marcelin, Boulesteix & Georgelin (1985) detected the Shapley wing of the SMC in H$`\alpha `$ at the level of 4 R. These H$`\alpha `$ detections, which are partly confirmed by the UKST H$`\alpha `$ survey (Zealey 2001, personal communication), are more than an order of magnitude stronger than expected in our simple model proposed. It seems that some of the bright H$`\alpha `$ regions arise from embedded UV sources, e.g.,  UV-bright stars in the Shapley wing. However, this does not appear to be the explanation for recent detections of bright H$`\alpha `$ along the Magellanic Stream (Weiner & Williams 1996; Weiner et al. 2001, these proceedings).

Since the Magellanic Stream passes directly over the South Galactic Pole (illustrated in Fig. 7), BM99 had hoped this would constrain $`\widehat{f}_{\mathrm{esc}}`$ from the Galaxy but their original model underestimates the required flux by at least a factor of 5. Weiner et al. (2001) claim detections to some points of the Stream as high as 1 R. The Galaxy fails to produce this level of ionization by at least an order of magnitude. The expected H$`\alpha `$ levels produced by the Galaxy are shown in Fig. 8 where we also show the additional contribution from the LMC. Since the LMC and much of the Stream lie close to $`Y=0`$, the model predictions are given in this plane. Bland-Hawthorn & Putman (2001) discuss various dynamical scenarios which may account for the bright Stream detections.

In summary, the Magellanic Stream poses a clear problem for the H$`\alpha `$ distance method since, in a few locations, it appears to be an order of magnitude brighter than predicted by the model. The Magellanic Bridge is also H$`\alpha `$ bright in places, but this may be due to internal UV sources. However, the HVCs appear to be broadly consistent with the spiral arm ionization model, and we find that most HVCs detected to date are scattered throughout the halo up to distances of 50 kpc from the Sun. Our suspicion is that most of the HVCs, like the Magellanic Bridge and the Magellanic Stream, result from galaxy interactions, specifically dwarf galaxies disrupted by the Galaxy. An interesting prospect that we have been investigating with B.K. Gibson is that some of the HVCs have been dislodged from the outer disk due to the disruptive passage of a nearby dwarf. There are several observations which support this: (i) the metallicities are consistent with the outer disk; (ii) the outer HI disk looks highly disturbed (Burton & Te Lintel Hekkert 1986); (iii) in some instances the HVCs show continuity in velocity with the outer disk (e.g.,  Gibson et al. 2001); and (iv) the H$`\alpha `$ distances place many of the clouds on 10 kpc scales. We have begun to carry out hydrodynamical simulations in order to test this idea further.

## References

Bica, E. et al. 1996, A&A, 313, 405

Binette, L. et al. 1994, A&A, 292, 13

Bland-Hawthorn, J. et al. 1998, MNRAS, 299, 611

Bland-Hawthorn, J. & Maloney, P.R. 1999$`b`$, ASP Conf. 166, In Stromlo Workshop on HVCs, eds. B.K. Gibson & M.E. Putman, 212

Bland-Hawthorn, J. & Maloney, P.R. 1999, ApJ, 510, L33 (BM99)

Bland-Hawthorn, J. & Putman, M.E. 2001, In Gas & Galaxy Evolution, eds. Hibbard, Rupen & van Gorkom, 240, 369–382 (astro-ph/0110043).

Blitz, L., Spergel, D. N., Teuben, P. J., Hartmann, D. & Burton, W.B. 1999, ApJ, 514, 818

Burton, W.B. & Te Lintel Hekkert, P. 1986, A&AS, 65, 427

Braun, R. & Burton, W.B. 2000, A&A, 354, 853

Cianci, S. 2001, PhD, Sydney University

Domgörgen, H. & Mathis, J.S., 1994, ApJ, 428, 647

Dwek, E. et al. 1995, ApJ, 445, 716

Fujimoto, M. & Sofue, Y. 1976, A&A, 47, 263

Glazebrook & Bland-Hawthorn, J. 2001, PASP, 113, 197

Greggio, L. & Renzini, A. 1999, Mem. Soc. Astron. Ital., 70, 691

Haffner, L.M., Reynolds, R.J. & Tufte, S.L. 1999, ApJ, 523, 223

Haffner, L.M., Reynolds, R.J. & Tufte, S.L. 2001, ApJ, 556, L33

Israel, F.P. & Koornneef, J. 1979, ApJ, 230, 390

Johnson, P.G., Meaburn, J. & Osman, A.M.I. 1982, MNRAS, 198, 985

Kennicutt, R.C. et al. 1995, AJ, 109, 594

Klypin, A., Kravtsov, A.V., Valenzuela, O. & Prada, F. 1999, ApJ, 522, 82

Lockman, F.J. 1984, ApJ, 283, 90

Lyne, A.G., Manchester, R.N. & Taylor, J.H. 1985, MNRAS, 213, 613

Macchetto, F. et al. 1996, A&AS, 120, 463

Manchester, R.N. & Taylor, J.H. 2001, Pulsars (W.H. Freeman: San Francisco).

Marcelin, M., Boulesteix, J. & Georgelin, Y. 1985, Nat, 316, 705

Mathewson, D.S. & Ford, V.L. 1984, in Structure and Evolution of the Magellanic Clouds, IAU Symp., 108, 125

McGee, R.X., Brooks, S.W. & Batchelor, R.A. 1972, Aust. J. Phys., 25, 581

Moore, B., Ghigna, S., Governato, F., Lake, G., Quinn, T., Stadel, J. & Tozzi, P. 1999, ApJ, 524, L19

Nordgren, T.E., Cordes, J.M. & Terzian, Y. 1992, AJ, 104, 1465

O’Connell, R.W. 1999, ARAA, 37, 603

Ortiz, O. & Lepine, J.R.D. 1993, A&A, 279, 90

Parker, J.W. 1993, AJ, 106, 560

Putman, M.E. 2000, PhD, Australia National University

Reynolds, R.J. 1984, ApJ, 282, 191

Reynolds, R.J. 1989, ApJ, 345, 811

Reynolds, R.J. & Tufte, S.L. 1995, ApJ, 439, L17

Reynolds, R.J., Haffner, L.M. & Tufte, S.L. 1999, ApJ 525, L21

Ryans, R.S.I. et al. 1997, MNRAS, 289, 83

Smith, A.P., Cornett, D.H., & Hill, R.S. 1987, ApJ, 320, 609

Taylor, J.H. & Cordes, J.M. 1993, ApJ, 411, 674

Tufte, S.L., Reynolds, R.J. & Haffner, L.M. 1998, ApJ, 504, 773

van Woerden, H. et al. 1999, Nat, 400, 138

Wakker, B.P. & van Woerden, H. 1991, A&A 250, 509

Wakker, B.P. & van Woerden, H. 1997, ARAA, 35, 217

Wakker, B.P. 2001, ApJS, in press (astro-ph/0102147)

Walborn, N.R. 1984, in Structure and Evolution of the Magellanic Clouds, IAU 108, 243

Weiner, B.J. & Williams, T.B. 1996, AJ, 111, 1156

Weiner, B.J., Vogel, S.N. & Williams, T.B. 2001, these proceedings(astro-ph/0109055)

Wyse, R.F.G. 1999, ASP Conf. Ser. 165: The Third Stromlo Symposium: The Galactic Halo, 1

Zhao, H. 1996, MNRAS, 283, 149

Zurita, A. et al. 2001, A&A, in press (astro-ph/0108135)

Zwaan, M. & Briggs, F. 2000, A&A, 530, L61


text:
# Stationary Regime of Random Resistor Networks Under Biased Percolation

## I Introduction and model

Electrical breakdown of disordered media has been widely studied in the last twenty years -. This is due to its relevant implications on material technology and on fundamental aspects related to the understanding of the response properties of disordered systems to high external stresses . It is well known that the application of a finite stress (electrical or mechanical) to a disordered material generally implies a nonlinear response, which leads to a catastrophic behavior in the high stress limit . Percolation theory provides a powerful approach for studying breakdown phenomena of disordered media . In this framework, several models have been proposed to describe the electrical breakdown of granular metals and of conductor-insulator composites in terms of critical phenomena near the percolation threshold -. Furthermore, the critical exponents characterizing these behaviors have been both theoretically calculated and measured in severals materials ,-, . Nevertheless, few attempts have been made so far to describe the behavior of disordered media over the full range of applied stress . The understanding of breakdown phenomena in the full dynamical regime is thus unsatisfactory to the present time. On the other hand, important information is expected from such a study, like precursor effects, rôle of disorder, predictability of breakdown, etc. .

The aim of this paper is to present a percolative model of sufficient generality to address the above issues. To this purpose, we consider a random resistor network (RRN)  in which two competing percolations are present, defect generation and defect recovery, which determine the values of the elementary network resistances. Both processes are driven by an external current and by the heat exchange between the network and the thermal bath. Monte Carlo simulations are performed to explore the network evolution in a wide range of bias values. A stationary state or an irreversible breakdown of the RRN can be reached depending on the value of the applied current. By focusing on the steady state, we study the resistance and the resistance noise properties. We found that the average network resistance and the relative resistance noise scale with the ratio of the applied current to the current value corresponding to the onset of nonlinearity. These results are discussed in connection with measurements in composite materials and in conducting polymers .

We consider a two-dimensional, square-lattice RRN of total resistance $`R`$, made of $`N_{tot}`$ resistors with resistance $`r_n`$. We take a square geometry $`N\times N`$, where $`N`$ determines the linear size of the network. A constant current $`I`$ is applied through electrical contacts realized by perfectly conducting bars at the left and right sides of the network. As a consequence, a current $`i_n`$ is flowing in the nth 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 temperature $`T_n`$, according to:

$$r_n(T_n)=r_0[1+\alpha (T_nT_0)]$$
(1)
where $`r_0`$ is the resistance value of the elementary resistor at the temperature $`T_0`$ and $`\alpha `$ is the temperature coefficient of the resistance. The local temperatures are calculated as in Ref. :

$$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)
In this expression, $`N_{neig}`$ is the number of first neighbours of 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 notice that Eq. (2) implies an instantaneous thermalization of each resistor at the value $`T_n`$, and then, by adopting Eq. (2), we are neglecting for simplicity time dependent effects which are discussed in Ref. .

In the initial state of the network, $`I=0`$, $`T_nT_0`$ and all the resistors are identical $`r_nr_0`$. The evolution of the RRN arises from the presence of two competing percolations. The first consists of generating fully insulating defects (broken resistors). This process occurs with probability:

$$W_D=\mathrm{exp}[E_D/K_BT_n]$$
(3)
where $`E_D`$ is a characteristic activation energy and $`K_B`$ the Boltzmann constant . The second percolation consists of recovering the insulating defects. This process occurs with a probability $`W_R`$ expressed as in Eq. (3) but with a different activation energy $`E_R`$. As a result of the competition between these two percolations, depending on the parameters which specify the physical properties of the system and depending on the external conditions (bias current and bath temperature), the RRN can reach a steady state or the percolation threshold. In the first case, $`R`$ fluctuates around its average value $`<R>`$, while in the second case, an irreversible breakdown occurs, 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 . By focusing on the effect of the external current, we define $`I_b`$ as the greatest current value for which the RNN is stationary. We notice that for biased percolation the following expression, $`E_R<E_D+k_BT\mathrm{ln}[1+\mathrm{exp}(E_D/K_BT_0)]`$, represents a necessary condition for the existence of a steady state . Monte Carlo simulations are performed according to the following procedure: (i) starting from the perfect lattice with given local currents and temperatures, $`i_n`$ and $`T_n`$, (ii) resistances $`r_n`$ are changed according to Eq. (1) and defects are generated with probability $`W_D`$; (iii) the currents $`i_n`$ are calculated by solving Kirchhoff’s loop equations; the local temperatures are updated; (iv) the temperature dependence of the resistances $`r_n`$ is again accounted for and defects are recovered with probability $`W_R`$; (v) $`R`$, $`i_n`$ and $`T_n`$ are finally calculated and this procedure is iterated from (ii) until electrical breakdown or steady state is achieved. In the last case the iteration runs long enough to allow a fluctuation analysis to be performed. Each iteration step can be associated with an elementary time step on an appropriate time scale (to be calibrated with experiments). As reasonable values of the parameters, simulations have been performed by taking: $`N=75`$ (except when differently specified), $`T_0=300`$ (K), $`\alpha =10^3`$ (K<sup>-1</sup>), $`A=5\times 10^5`$ (K/W), $`E_D=0.17`$ (eV). Several values of $`E_R`$ and $`r_0`$ have been considered: $`0.026E_R0.16`$ (eV) and $`1r_010`$ ($`\mathrm{\Omega }`$). The values of the external current range between $`0.001`$ and $`3`$ (A).

## II Results

We show in Fig. 1 a picture of the RRN near the electrical breakdown. In this case we have taken: $`N=45`$, $`r_0=1`$ ($`\mathrm{\Omega }`$), $`E_R=0.10`$ (eV) and $`I=0.5`$ (A), i.e. $`I>I_b=0.45`$ (A). The different levels of gray correspond to different values of $`r_n`$. We can clearly see that, with respect to the initial state (perfect network), the network has evolved to a disordered state associated with the formation and growth of a channel of broken resistors elongated in the direction perpendicular to the applied current. This kind of damage pattern well reproduces the experimental pattern observed in metallic lines failed as a consequence of electromigration . Typical evolutions of $`R`$ are shown in Fig. 2 at increasing bias values. In this case $`N=75`$ while all the other parameters are the same of Fig. 1. Thinner curves refer to steady state regime while the thicker curve refers to a RRN undergoing electrical breakdown ($`I>I_b=0.75`$). In the steady state two regimes can be identified: an Ohmic regime (lower two curves) and a nonlinear regime characterized by a significant increase of resistance (remaining curves). By focusing on the steady state regime, we report in Fig. 3 the average resistance $`<R>`$ as a function of the applied current. The different curves correspond to different values of $`r_0`$, i.e. to RRNs of different initial resistance, while the recovery activation energy is $`E_R=0.026`$ (eV). Each value have been calculated by considering the time average on a single realization and then making the ensemble average over 20 realizations. At low biases the average resistance takes a constant value $`<R>_0`$ which represents the intrinsic linear response property of the network (Ohmic regime). When $`I`$ is above a threshold value $`I_0`$ (threshold for the nonlinearity onset), $`<R>`$ increases with bias until the applied current reaches the $`I_b`$ value, above which the RRN undergoes electrical breakdown. Thus, in the following, we indicate with $`<R>_b`$ the average value of $`R`$ at $`I=I_b`$, i.e. the last stable value of the resistance. Figure 3 also shows that by increasing $`r_0`$ and thus the initial network resistance, both $`I_0`$ and $`I_b`$ decrease. Precisely, we have found: $`I_bR_0^\delta `$ and $`I_0R_0^\delta `$ with $`\delta =0.51\pm 0.01`$. Therefore the ratio $`I_b/I_0=28\pm 1`$ is independent of the initial network resistance. Moreover, we have also found $`<R>_b/<R>_0=1.85\pm 0.08`$. The effect of the recovery activation energy on the steady state is shown in Fig. 4, which reports the ratio $`<R>/<R>_0`$ as a function of the applied current for different values of $`E_R`$ (in this case all the curves are obtained for $`r_0=1`$ ($`\mathrm{\Omega }`$)). The overall behavior is similar to that shown in Fig. 3: an Ohmic regime at low bias is followed by a nonlinear regime for $`I>I_0`$. Moreover, by increasing $`E_R`$ both $`I_0`$ and $`I_b`$ decrease and the stability region is thus strongly reduced. Nevertheless, an important difference between the effect of varying the initial network resistance and that of varying $`E_R`$ is that, in the last case, the ratio $`I_b/I_0`$, exhibits a significant dependence on $`E_R`$, as shown in Fig. 5. To investigate the existence of scaling relations and their universality, Fig. 6 reports the log-log plot of the relative variation of the average resistance, $`(<R><R>_0)/<R>_0`$, as a function of the ratio $`I/I_0`$ for different values of $`r_0`$ and $`E_R`$. The plot shows that all these curves collapse onto a single one and that the relative variation of $`<R>`$ as a function of $`I/I_0`$ exhibits a power law behavior. We conclude that, the average resistance follows the scaling relation:

$$\frac{<R>}{<R>_0}=g(I/I_0),g(I/I_0)1+(I/I_0)^\theta $$
(4)
with the scaling exponent $`\theta =2.1\pm 0.1`$ being independent of both the initial resistance of the RRN and the recovery activation energy. Other simulations, performed on rectangular networks, confirm the same value for $`\theta `$. The quadratic dependence of $`<R>`$ on $`I`$, found here, can be understood in the spirit of mean-field theory when we consider that $`\mathrm{\Delta }R\alpha R_0\mathrm{\Delta }T`$ and $`\mathrm{\Delta }TAR_0I^2`$. Moreover, recalling previous results concerning the ratio $`I/I_0`$, Eq.(4) explains the independence of the ratio $`<R>_b/<R>_0`$ on the initial RRN resistance (Fig. 3) and, by contrast, its significant dependence on $`E_R`$, as shown in Fig. 4. All these results well agree with recent measurements in the Joule regime of carbon high-density polyethylene composites reported in Ref. .

The resistance fluctuations are analyzed for different values of $`E_R`$ and $`r_0`$. Figure 7 reports the relative variance of resistance fluctuations, $`\mathrm{\Sigma }<\mathrm{\Delta }R^2>/<R>^2`$, as a function of the external current. Curves 1, 2 and 3 (with full circles) show $`\mathrm{\Sigma }`$ for $`r_0=1`$ ($`\mathrm{\Omega }`$) and $`E_R=0.060`$, $`0.043`$, $`0.026`$ (eV) respectively, while the curves belonging to the set 3 are obtained for $`E_R=0.026`$ (eV) and different values of $`r_0`$. Figure 7 points out the existence of two different noise regimes. The first regime occurs for $`I<I_0`$, i.e. when Joule heating effects are negligible. This noise arises from two random percolations and represents an intrinsic noise of the RRN, depending only on the values of $`E_D`$ and $`E_R`$ . The second regime occurs when $`I>I_0`$ and the value of $`\mathrm{\Sigma }`$ is found to become strongly dependent on the external current. By plotting $`\mathrm{\Sigma }/\mathrm{\Sigma }_0`$ as a function of $`I/I_0`$ we have found that all the data of Fig. 7 collapse onto a single curve, as shown in Fig. 8. Moreover, a power law behavior is observed in the pre-breakdown region. Therefore, we can conclude that the relative variance of resistance fluctuations follows the scaling relation:

$$\frac{\mathrm{\Sigma }}{\mathrm{\Sigma }_0}=f(I/I_0),f(I/I_0)1+(I/I_0)^\eta $$
(5)
where the scaling exponent is $`\eta =4.1\pm 0.1`$. This value of $`\eta `$ agrees with the values obtained from electrical noise measurements in conducting polymers .

In conclusion, we have studied by Monte Carlo simulations the stationary regime of RRNs resulting from the simultaneous evolutions of two competing percolations. The two percolations consist of generating (recovering) fully insulating defects which are driven by an external current and by the heat exchange with a thermal bath. We have analyzed the behavior of the average resistance and of the relative variance of resistance fluctuations over a wide range of the applied current and as a function of different model parameters. We have found that both these quantities follow a scaling relation in terms of the ratio between the applied current and the current value corresponding to the nonlinearity onset. Both scaling exponents are found to be independent of the model parameters. These results compare well with resistance measurements in composite materials performed in the Joule regime up to breakdown  and with noise measurements in conducting polymers .

###### Acknowledgements.

This research is performed within the STATE project of INFM. Partial support is also provided by ASI project, contract I/R/056/01.

FIGURE CAPTIONS

Fig. 1 - Network evolution near the electrical breakdown. Different gray levels from black to white correspond to increasing values of $`r_n`$ from 1 to 3 ($`\mathrm{\Omega }`$).

Fig. 2 - Resistance as a function of time at increasing bias values. Thinner curves correspond to steady state regime and they are obtained, going from bottom to top, for $`I=0.01`$, $`0.05`$, $`0.10`$, $`0.35`$, $`0.70`$, $`0.75`$ (A). The thicker curve is obtained for $`I=0.78`$ (A) and it corresponds to a RRN undergoing electrical breakdown.

Fig. 3 - Normalized average resistance versus external bias. We take $`E_R=0.026`$ (eV) and $`E_D=0.167`$ (eV), while the value of $`r_0`$ ranges between 1 $`\mathrm{\Omega }`$ and 10 $`\mathrm{\Omega }`$.

Fig. 4 - Normalized average resistance versus external bias. We take $`r_0=1`$ $`(\mathrm{\Omega })`$, $`E_D=0.167`$ (eV) while the values of $`E_R`$ range between $`0.026`$ (eV) and $`0.155`$ (eV).

Fig. 5 - Plot of the ratio $`I_b/I_0`$ as function of the recovery activation energy. The curve is a fit with a power law: $`I_b/I_0E_R^{0.36}`$.

Fig. 6 - Log-log plot of the relative variation of resistance versus $`I/I_0`$. Data shown in this figure are the same of those reported in Figs. 3 and 4.

Fig. 7 - Relative variance of resistance fluctuations as a function of the external bias. Curves 1, 2, 3 refer to $`r_0=1.0`$ ($`\mathrm{\Omega }`$) and $`E_R=0.060`$, $`0.043`$, $`0.026`$ eV, respectively. The five curves belonging to set 3 are obtained with $`r_0=1.0`$ ($`\mathrm{\Omega }`$) (full circles), $`2.5`$ ($`\mathrm{\Omega }`$) (open squares), $`5.0`$ ($`\mathrm{\Omega }`$) (full triangles), $`7.5`$ ($`\mathrm{\Omega }`$) (open triangles), $`10.0`$ ($`\mathrm{\Omega }`$) (full diamonds).

Fig. 8 - Log-log plot of the relative variance of resistance fluctuations normalized to the same quantity calculated in the linear regime versus $`I/I_0`$. A power law fit is shown in the pre-breakdown regime.