text:
# The photometric-amplitude and mass-ratio distributions of contact binary stars

## 1 INTRODUCTION

Contact binary stars (also called W UMa-type variable stars) have a unique property among eclipsing binaries in that their geometrical effects are much more important in defining the shapes of their light curves than the atmospheric properties of their components. Because of the similarity of effective temperatures of the components, the eclipse depths are practically independent of the effective temperatures, but depend strongly on geometrical parameters such as the orbital inclination, $`i`$, the degree of contact, $`f`$, and the mass ratio, $`q`$. This situation is very much different than for detached binaries where – in most cases – differences in component temperatures strongly affect the eclipse depths. Also, for a fixed mass ratio, just one parameter (the potential) defines the relative size of both components in a contact binary, in place of two independent radii as in detached binaries, so that the light curves are described by fewer free parameters.

The simplicity of the contact-binary light curves and lack of features in them (absence of external – and in most cases – internal eclipse contacts) mean that the light curves carry little information. Extensive calculations of large grids of light curves (Rucinski, 1993a, Paper I) have shown that – in the general case – only two of the three main geometrical parameters are well constrained by the light curve shapes. Only when radial-velocity information on the mass ratio $`q`$ is available can the light curve analysis also yield the $`(i,f)`$ pair and the whole set of the parameters can be determined. From time to time one sees attempts at determination of the mass ratio through analysis of the $`\chi ^2(q)`$ curve of the light-curve fit. This is a very dangerous approach because such results on $`q`$ are very poorly constrained and, in fact, frequently plainly wrong, as the new radial velocity data for many contact binaries (Lu & Rucinski, 1999; Rucinski & Lu, 1999; Rucinski et al., 2000; Lu et al., 2001) have shown: In numerous cases a spectroscopic (i.e. correct) value of the mass ratio is far away from the “best” photometric value of $`q`$, frequently far beyond the stated errors estimated from the shape of the $`\chi ^2(q)`$ curve. An exception to this indeterminacy is the special case of the totally eclipsing contact systems. Such light curves show characteristic inner contacts with duration of totality setting a strong constraint on the $`(q,i)`$ pair (Mochnacki & Doughty, 1972a, b).

Instead of considering results of individual light curve solutions, augmented by radial velocity studies, one can attempt to utilize statistical distributions of the observed parameters for large samples of binaries. The simplest to obtain and most obvious among such distribution is the photometric amplitude distribution. However, it has never been used, probably because it is generally recognized that such an observational distribution must be very closely related to (if not observationally defined by) discovery selection effects. Now the situation is changing: New data from microlensing surveys and from massive, deep sky surveys, aimed at other goals, but extensive and statistically rigorous, start providing a wealth of sound “by-product” information for contact binaries (Rucinski, 1997a, b; Alcock et al., 1997; Rucinski, 1998a, b, 1999a; Maceroni & Rucinski, 1999; Rucinski & Maceroni, 2001). Such surveys, with well controlled selection biases, are far superior to any statistical inferences based on individual light and radial velocity solutions as the latter are very heavily biased by discovery and observational selection effects and by observer preferences.

This paper is basically a much expanded version of a preliminary discussion in Section 5 of Rucinski (1997b). It attempts to characterize and analyze information contained in light variation amplitudes of contact binaries. By definition, we call the amplitude the depth – in magnitudes – of the deeper eclipse. The amplitude distribution is considered here as a tool to study two important issues: (1) How many contact systems are missed and remain unknown due to their amplitudes being smaller than the detection threshold, and (2) Can such a distribution be used to derive properties of contact systems, especially the distribution of the mass ratios? Both questions can be addressed for the idealized case of isolated contact binaries. However, because of the blending problems which affect the currently-best statistics for the Baade’s Window data and – even more importantly – because contact binaries frequently have close, unresolved companions, the results of this paper must be considered as preliminary and illustrative.

The paper continues the approach of Paper I in that we analyze the main, global properties of the light curves by covering the entire range of relevant parameters. While the approach may appear to be coarse, the intention is to uncover the main relations and dependencies. We feel that this is justified in view of lack of similar global approaches in this field.

We give an expanded background in Section 2 while the results of the calculations are presented in Section 3. Applications of the theoretical results to observational statistics based on the data for the OGLE Baade’s Window sample, to derive the mass-ratio distribution, are discussed in Section 4. This distribution cannot yet be analyzed at its low mass-ratio end because of the poor statistics for small-amplitude systems. However, we predict that the expected cut-off to the mass-ratio distribution at $`q_{min}`$ will strongly modify the amplitude distribution at its low end; the matter of the smallest possible mass ratios, $`q_{min}`$, is discussed in Section 5. The results of the paper and their implications for our understanding of the evolution of contact binary systems are discussed in Section 6. The conclusions of the paper are summarized in Section 7.

## 2 BACKGROUND ON PREDICTION OF AMPLITUDE DISTRIBUTIONS

### 2.1 The contact binary model

Current light-curve synthesis models permit computation of light curves with almost arbitrarily high accuracy. The accuracy depends basically only on how many integration points one is ready to use. Obviously, the model must represent real contact binaries, a requirement which – with the improved numerical accuracy of the models – is becoming increasingly difficult to fulfill in view of many complications such as activity-induced star spots or possible deviations from plane-parallel atmosphere models in the “neck” region between components of a contact binary. In this paper, we neglect all complications of spots and stellar activity and concentrate on the basic dependencies of light-variations on the geometrical parameters. In order to do so, we fix the atmospheric parameters, such as the limb and gravity darkening coefficients. It may be argued that they are relatively less important than the geometrical parameters because – due to constancy of the effective temperature over the whole common surface – these coefficients are practically the same everywhere. Thus, as in Paper I, we consider only one representative wavelength (matched to the $`V`$-band filter) and one effective temperature for a typical solar-type star (5770 K). Extensive tests have shown that such a combination gives reasonable representation of the observed light curves of contact binaries for yellow-red spectral region (in this paper we utilize observational amplitude distributions from the OGLE sample in the photometric $`I`$ band) and for a wide range of effective temperatures at which the contact binaries are most common, i.e. with spectral types from early F through G to early K.

The geometrical elements which define the shapes of the contact-binary light curves are designated as $`(i,f,q)`$ (Rucinski, 1985, 1993b). The orbital inclination, $`0i90^{}`$, is conventionally defined in such a way that $`i=0^{}`$ corresponds to the orbit in the plane of the sky and $`i=90^{}`$ corresponds to the orbit seen edge on. The simplest assumption concerning the orbits is that they are randomly oriented in space. This means that probability of detection of a binary with a given $`i`$ is proportional to $`\mathrm{sin}i`$. This is because for $`i0^{}`$, the available range of orientations of the orbital ascending node goes to zero. The degree of contact is usually defined through the potential of the common equipotential surface. In this paper we consider only three discrete values of the degree-of-contact parameter, $`f=(C_1C)/(C_1C_2)`$, where $`C`$ is the Jacobi constant for a common equipotential: $`f=0`$, $`f=0.25`$ and $`f=0.5`$. The inner contact with the components just touching corresponds to $`f=0`$. The frequency distribution of the parameter $`f`$, $`F(f)`$, is currently unknown. However, analysis of a large number of contact systems in the OGLE sample (Rucinski, 1997b) has confirmed the earlier indications (Lucy, 1973; Rucinski, 1973) that $`f`$ is typically small, within $`0<f<0.5`$, with a tendency for a (poorly defined) maximum around $`f=0.25`$. Large values of $`f1`$, corresponding to the outer contact, are not observed.

The third geometrical parameter is the mass ratio. From the point of view of the contact binary formation and evolution this is astrophysically the most important parameter among those that define the light curve shape. We define the mass ratio as $`q=M_2/M_11`$. From many radial velocity studies – which tend to preferentially favor equal-mass combinations – we know that contact systems definitely avoid the equal-mass situation of $`q=1`$ and this fact has a good theoretical explanation (see Section 6). Systems with $`q`$ approaching unity do exist; an extreme example is the recently discovered system with very good radial velocity data, V753 Mon, which has $`q=0.970\pm 0.003`$, Rucinski et al. (2000). However – in spite of the relative ease of detection of large mass-ratio systems – we see very few systems with $`q>0.8`$, so they must be very rare in space. Currently, the only property of the mass-ratio distribution, $`Q(q)`$, that we can be sure of is its tendency to increase for $`q0`$, but the rate of the increase is not known. This paper attempts to determine $`Q(q)`$ in Sections 35; the results will be related to the theoretical evolutionary predictions in Section 6.

As the reader may have already noticed, we use capital letters to denote the distributions of parameters denoted by lower case letters. In particular, the amplitude distribution can be written as $`A(a)da=(dN/da)da`$ so that $`A(a)`$ takes the significance of the number of binary systems observed in the elemental increment $`da`$.

### 2.2 Derivation of the amplitude distribution from light curve synthesis models

The amplitudes predicted from a light-curve model can be written as three-dimensional functions: $`a=a(i,q,f)`$. In practice, we use arrays calculated for fine grids of the parameters $`i`$, $`q`$, $`f`$. To simplify handling of such large arrays, and because we know that the degree-of-contact is not very important in defining the light curve shapes, we have computed such arrays for only a few fixed values of $`f`$. Each amplitude distribution for a fixed value of $`q`$, $`A_q(a)`$, can be also considered individually; such normalized distributions will be used in an attempt to derive the mass-ratio distribution $`Q(q)`$ (Section 4). For the cases in which discrete values of $`f`$ and $`q`$ are used, we will use the notation $`a_q^f`$. In particular, we will consider the functional dependence of the amplitude on the orbital inclination, $`a_q^f=a_q^f(i)`$ as well as the distribution of the amplitudes, $`A_q^f=A_q^f(a)`$.

We derive the expected amplitude distribution, $`A_q^f(a)`$, by eliminating the orbital-inclination dependence under the assumption of random orientation of the orbital planes in space. Two approaches are possible, both utilizing light curves computed with the light curves synthesis codes. One is a semi-analytical approach: We start with the dependence $`a=a(i)`$ computed by varying the inclination in the models (for a fixed combination of $`q`$ and $`f`$). For a random distribution of orientations, the number of discovered systems should scale with the inclination angle $`i`$ as: $`dN=\mathrm{sin}idi`$, which can be written as $`(dN/da)(da/di)=\mathrm{sin}i`$, hence

$$A(a)=(dN/da)=\mathrm{sin}i(a)(da/di)^1$$
(1)
The amplitudes monotonically increase with $`i`$ so that there is no difficulty in inverting the function $`a(i)`$ into $`i(a)`$. However, the reciprocal of the derivative $`da/di`$ must be computed numerically. This is a complication because the light-curve synthesis results on $`a(i)`$ carry numerical noise which becomes amplified in the numerical differentiation.

We can get some insight into the behavior of $`da/di`$ by analyzing one specific case, as shown in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars. The left panel of the figure shows 45 light curves computed with inclinations varied in 2-degree steps for a representative case of $`q=0.35`$ and $`f=0.25`$. The corresponding variations of $`a(i)`$ and $`da/di`$ are shown in the right panel. Notice the progressive increase of the amplitude with inclination to about $`i70^{}`$ where the derivative has its peak and then flattening of the $`a(i)`$ curve in the region of the total eclipses. From the shape of the $`(da/di)`$ curve, we can expect two regions at both ends of the distribution where $`(da/di)^1`$ will be particularly difficult to evaluate, one for $`i0^{}`$ and one for $`i90^{}`$. The region at $`i90^{}`$ is usually small, but its size depends on the mass ratio and grows for small values of $`q`$. This region is populated by systems showing total eclipses of similar depths for orbital inclinations close to $`90^{}`$. At the other end, for $`i0^{}`$, the light curves have very small amplitudes which change very little with $`i`$. The resulting increase $`(da/di)^1\mathrm{}`$ is expected to be moderated in Equation 1 by the term $`\mathrm{sin}i0`$, so that the result is difficult to predict. We show later on in Section 3 that – somewhat unexpectedly – for $`a0`$, the amplitude distributions for all values of $`q`$ appear to converge to mass-ratio dependent constant values.

A much simpler second approach to evaluate $`A(a)`$, which avoids the numerical difficulties of computation of $`(da/di)^1`$, is through a Monte Carlo simulation. In such an experiment, the inclinations are drawn randomly with a distribution $`I(i)\mathrm{sin}i`$. These are then interpolated in the amplitude arrays $`a_q^f=a_q^f(i)`$, to form the amplitude distributions by simply binning the resulting distributions, $`A(a)=(dN/da)`$. This was the approach actually used in this paper. It has the advantage of direct modeling of the difficult regions at $`i0^{}`$ and $`i90^{}`$, and entirely avoids evaluation of $`(da/di)^1`$.

We show in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars a representative case of the amplitude distribution $`A(a)`$ computed using the Monte-Carlo approach for the same light curves as in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars, that is for the mass ratio $`q=0.35`$ and for the degree-of contact parameter $`f=0.25`$ (this figure shows also the corresponding distributions for $`f=0`$ and $`f=0.5`$). We will discuss an amplitude distribution like this one in the next Section 3, after summarizing the details and assumptions of the Monte Carlo computations.

## 3 CALCULATIONS OF THE AMPLITUDE DISTRIBUTION

### 3.1 Details of the Monte Carlo calculations

For consistency with previous calculations of the large grid of light curves in Paper I, the light-curve grid used to predict the amplitude distributions utilized the same assumptions on the temperatures, limb and gravity darkening laws and relative fluxes. We suggest the reader consult that paper for details. The only differences are as follows: (1) We improved the sampling of the degree-of-contact parameter, $`f`$, by adding the most likely value of $`f=0.25`$, in addition to the two bracketing values considered before, $`f=0`$ (the inner contact) and $`f=0.5`$; we have also dropped the unobserved case of $`f=1`$ (the outer contact); (2) The amplitude calculations have been extended over the whole width of the inclination angles, $`0^{}i90^{}`$ (the previous calculations applied only to $`i30^{}`$); (3) While the parameter $`f`$ has been sampled rather coarsely, the parameters $`q`$ and $`i`$ have been sampled with fine grids of $`\mathrm{\Delta }q=0.01`$ and $`\mathrm{\Delta }i=1^{}`$; (4) The amplitude distributions have been computed using the magnitude scale (note that the light intensity was used in Paper I), with a fine grid of $`\mathrm{\Delta }a=0.01`$. Later on, for practical applications, larger bins of $`\mathrm{\Delta }a=0.05`$ have been used.

In the Monte Carlo experiments, the inclinations were drawn using uniformly distributed random numbers in the $`x[0,1]`$ range, mapped into the $`\mathrm{sin}i`$-distributed inclination angles through $`i=\mathrm{arccos}x`$. Each experiment for a fixed pair of $`(f,q)`$ involved $`4\times 10^6`$ samples of $`a_q^f=a_q^f(i)`$ which were then binned into $`A_q^f(a)`$ distributions over $`0a1.1`$, with 110 bins of $`\mathrm{\Delta }a=0.01`$. Thus, typically, each bin $`\mathrm{\Delta }a`$ contained $`4\times 10^4`$ samples giving an error at the level of about 0.5 percent. Results of more accurate computations of the normalized $`A_q(a)`$ with larger bins of $`\mathrm{\Delta }a=0.05`$ and for 10 values of the mass ratio in intervals of $`\mathrm{\Delta }q=0.1`$ are given in Table 1 for the case of $`f=0.25`$. The wider bins in $`a`$ and $`q`$ result in an increase of accuracy by a factor of about 4 times so that the $`A_q(a)`$ distributions in Table 1 should be accurate to about 0.1 percent. The distributions for $`f=0`$ and $`f=0.5`$, have been also computed, but we do not give them here for economy of space111The tables similar to Table 1 for $`f=0`$ and $`f=0.5`$, as well as the detailed tables sampled at $`\mathrm{\Delta }q=0.01`$ and $`\mathrm{\Delta }a=0.01`$ are available over the Internet from http://www.astro.utoronto.ca/$``$rucinski..

Figure The photometric-amplitude and mass-ratio distributions of contact binary stars shows the detailed behavior of the amplitude distributions $`A_{0.35}^{0.0}(a)`$, $`A_{0.35}^{0.25}(a)`$, $`A_{0.35}^{0.5}(a)`$, i.e. for a case of $`q=0.35`$, and for $`f=0`$, 0.25, 0.5. The narrow interval of $`\mathrm{\Delta }q=0.01`$ around $`q=0.35`$ has been used in the figure for illustrative purposes because this way we can clearly demonstrate the narrow peak of $`A(a)`$ at large amplitudes which corresponds to the range of inclinations giving total eclipses. When larger bins in $`q`$ are considered, as in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars (see the panel labeled “0.35” for a direct comparison with Figure The photometric-amplitude and mass-ratio distributions of contact binary stars), the peak becomes broader and less prominent. Later on, in actual applications in Sections 4 and 5, we use $`\mathrm{\Delta }q=0.1`$ or $`\mathrm{\Delta }q=0.2`$. The location of the total-eclipse peak weakly depends on the degree-of-contact, $`f`$: The amplitudes are larger for contact systems with stronger contact and smallest for binaries with components just touching.

### 3.2 Dependence of the amplitude distribution on the mass ratio

Figure The photometric-amplitude and mass-ratio distributions of contact binary stars shows the ten predicted amplitude distributions $`A_q^{0.25}(a)`$, sampled at $`\mathrm{\Delta }a=0.05`$, for $`f=0.25`$, with $`q`$ spanning the whole range $`0<q<1`$ in intervals of $`\mathrm{\Delta }q=0.1`$. The individual panels of the figure correspond to columns in Table 1; each $`A_q(a)`$ has been separately normalized to unity. Note the same details as for the specific case of $`A_{0.35}(a)`$ in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars, but now displayed for the whole range of the mass ratio, in wider bins in $`q`$ and $`a`$. The figure illustrates the strong dependence of $`A_q(a)`$ on the mass ratio. While all fixed-$`q`$ distributions show a rise towards small amplitudes, the range of the amplitudes is very different: For large $`q`$, all amplitudes are possible, while for small $`q`$, only small ones are permitted. If one particular value of the mass ratio were to dominate in the mass-ratio distribution $`Q(q)`$, the characteristic features of its corresponding $`A_q(a)`$ would be present in the combined $`A(a)`$. The observed amplitude distributions appear to be featureless, so that the corresponding mass-ratio distributions must be smooth, too. As we discuss in Section 4, the only characteristic feature that we can be sure of is lack of large amplitude systems indicating a strong weighting towards small mass ratios.

In addition to all individual distributions of $`A_q(a)`$ in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars, we also show distributions resulting from a simple addition of all distributions over the whole range of $`q`$. This, in effect, corresponds to the case of $`Q(q)=const`$. Such a distribution also monotonously increases towards low amplitudes, as do the individual $`A_q(a)`$ distributions. The local peaks at the characteristic maximum amplitudes in individual $`A_q(a)`$ disappear in the combined $`A(a)`$ through “dilution” when summed with other distributions. Bearing in mind the independent normalizations of each $`A_q(a)`$ and for the combined distribution, we can see that all $`A(a)`$ distributions look very similar in the central portions of the amplitude range and show the strongest dependence of the shape on $`q`$ at low and large values of the amplitudes.

### 3.3 The detection threshold

Of particular importance to practical applications of our calculations is the behavior of $`A_q(a)`$ for small amplitudes, because this determines how many systems would be lost in photometric surveys with specific detection thresholds. As we can see in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars for the case of $`q=0.35`$, the amplitude distribution tends to converge to constant value, $`A(a0)0.020.03`$, per the $`\mathrm{\Delta }a=0.01`$ bin. The same convergence is observed for the wider bins used in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars, for all values of $`q`$. Thus, it appears that, for $`i0^{}`$, the $`\mathrm{sin}i`$ term in Equation 1 does not force $`A(a)0`$ for $`a0`$. In practical terms, it means that the first bin in any amplitude histogram, $`[0,\mathrm{\Delta }a]`$, will always contain some objects, irrespectively of the size of the bin. This is an important new result indicating that a substantial number of contact binaries may have small amplitudes below detection thresholds. For the specific case of $`q=0.35`$ illustrated in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars, if the threshold were at the (currently unachievable) $`a_{min}=0.01`$ mag., then the loss would be about 2 – 3 percent.

Figure The photometric-amplitude and mass-ratio distributions of contact binary stars shows – as a function of $`q`$ – the fraction of systems falling below an assumed value of the threshold amplitude, $`\mathrm{\Sigma }A_q(a<a_{min})`$. For the detection threshold of $`a_{min}=0.01`$ mag., only about 2 to 5 percent of systems would be typically lost; this can be verified by inspecting the first bin $`0<a<0.01`$ in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars. However, Figure The photometric-amplitude and mass-ratio distributions of contact binary stars indicates that for a more representative value of 0.1 mag., some 20 to 40 percent would remain undetected; for larger $`a_{min}`$ the fractional loss would be larger. Note that for small values of $`a_{min}`$, there exists a simple proportionality between the fraction of the undetected systems and the height of the threshold, $`\mathrm{\Sigma }A_q(a<a_{min})a_{min}`$, which is a direct consequence of the convergence of $`A_q(a)`$ to a constant value for $`a0`$.

The detection losses strongly depend on the mass ratio. Thus, the main question is, what are the mass ratios and how they distribute the losses due to the detection thresholds? As we will see later, the mass-ratio distribution appears to be steeply rising for $`q0`$ which makes the left hand edge of Figure The photometric-amplitude and mass-ratio distributions of contact binary stars particularly important for estimation of the fraction of undetected systems.

### 3.4 Maximum amplitudes

As one can see in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars, for each mass ratio, there exists a characteristic maximum amplitude, $`a_{max}`$. This maximum amplitude grows with $`q`$. The mere presence of large amplitudes in an observed distribution indicates that among contact binaries of the sample there exist ones with large mass ratios.

The maximum amplitudes are tabulated as a function of $`q`$ for the three values of the fill-out parameter $`f`$ in Table 2. They are also shown in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars. We will discuss in the next section that the well established and most trustworthy part of the observed $`A(a)`$ extends only above $`a0.3`$. As a result, nothing can be said about the mass-ratio distribution below $`q0.12`$. This inter-relation between $`a_{max}`$ and the accessible range of mass ratios emphasizes the importance of the good detection statistics at low amplitudes. We will see later in Section 5 that the definition of $`A(a)`$ in the region around $`0.1<a<0.3`$ is particularly important for shedding light on the expected low mass-ratio cut-off in $`Q(q)`$ at $`q_{min}0.070.1`$.

## 4 DETERMINATION OF THE MASS-RATIO DISTRIBUTION

### 4.1 Observational data: 7.5 magnitude-limit sample

The strong dependence of the amplitude distributions on the mass ratio can be utilized to derive the functional shape of $`Q(q)`$ for the observed systems. The most obvious source for $`A(a)`$ are the catalogue data for contact binaries, as listed in the General Catalogue of Variable Stars (Kholopov et al., 1985–1988, GCVS). However, this material must be used judiciously as it is known to be strongly affected by discovery and observational selection biases. We used the most recent, October 2000 version of the GCVS, which is only available electronically222The GCVS has been obtained from ftp://ftp.sai.msu.su/pub/groups/cluster/gcvs/iii/.. It consists of the main catalogue of three volumes, augmented by the “name lists” number 67 to 74. We added name list number 75 from the same source and cross-referenced the combined catalogue with the list of variable stars in the Hipparcos catalogue (ESA, 1997). At that point, the variability types were checked on the basis of the Hipparcos light curves and several incorrect or disputable assignments of types were found.

Our previous experience with the amplitude statistics of contact binaries, based on the GCVS data (Rucinski & Kałużny, 1994; Rucinski, 1997b), indicated a strong dominance of large amplitude variables for the simple reason that large amplitude variables are always easier to discover than small amplitude ones. Such an excess of large amplitude systems seemed implausible, even without confrontation with the results of this paper which predicts dominance of small amplitudes. To avoid the bias of only large-amplitude systems being discovered among fainter systems, an assumption has been made that the sky has been fully inspected for variability by numerous observers and by Hipparcos to a relatively bright level. This level has been selected to coincide with the astrometric completeness limit of the Hipparcos mission at $`V7.5`$. The limiting minimum amplitude for such a sample is unknown. Perhaps it is $`a_{min}0.1`$, although fainter systems were detected down to 0.03 mag. We feel that this is the most complete currently-available sample for the sky-field contact binaries. Deeper samples that we attempted to construct in the same way indicated the presence of discovery selection losses past $`V7.5`$.

An important limitation in the current context is the presence of unresolved physical companions of contact binaries leading to photometric blending and a decrease in observed amplitudes. This is a difficult and murky area as no statistics exist to evaluate importance of the blending. Some observers remarked that close companions surprisingly frequently accompany W UMa binaries (Rucinski & Kałużny, 1988; Chambliss, 1992; Hendry & Mochnacki, 1998). We also keep discovering them in our spectroscopic survey of short-period binary systems which is currently conducted at the David Dunlap Observatory (Lu & Rucinski, 1999; Rucinski & Lu, 1999; Rucinski et al., 2000; Lu et al., 2001). The spectroscopic detection is the most sensitive and least biased of all available techniques, although speckle interferometry surveys have been already more systematic in surveying all bright stars of the sky. The spectroscopic analysis is particularly easy when use is made of broadening-functions (Rucinski, 1999b) which permit separation of components even for rather large difference in brightness through very different spectral signatures of broad and sharp components (Lu et al., 2001).

A “third light” contribution relative to the maximum brightness of the contact system, $`x=L_3/(L_1+L_2)`$, is expected to change the true amplitude $`a`$ to the observed one $`a^{}=2.5\mathrm{log}[(10^{0.4a}+x)/(1+x)]`$. At present it is very difficult to quantify the influence of the unresolved companions on the observational amplitude distribution because we have no idea about frequency of occurrence of the companions and thus about the distribution of the quantity $`x`$. If contact binaries are formed preferentially in a hierarchical process which produces wide orbits first and leaves pre-stellar clouds of low angular momentum to form contact systems, the frequency of occurrence may be higher than for other stars. Further, if the close binary formation process has more to do with random pairing then $`x`$ should be – on the average – small, but if the process tends to prefer equal-mass components, then $`x1`$. We see contact binaries in systems with bright companions, such as 44 Boo, with companions comparable in brightness as in HT Vir, but we also see faint, low mass companions to systems like VW Cep, so that apparently all values of $`x`$ are possible. The nearby systems offer the best chance of detection of the companions so that the Hipparcos sample is probably the best one to start from.

We have created what we call the “7.5 magnitude limit” sample of contact binary systems on the basis of the merged GCVS and Hipparcos data for binary systems with periods shorter than one day. By utilizing partially unpublished spectroscopic results from the David Dunlap Observatory, we could clean the sample of short-period pulsating stars and take into account the presence of close companions. The full discussion of the sample, which consists of 41 close binary systems, in that 10 systems having visual, speckle or spectroscopic companions, will be a subject of a separate investigation. We summarize here only the preliminary conclusions based on this sample concerning the amplitude distribution.

The sample includes all binaries designated in catalogues by codes “EW”, “EB” or “Ell” to the limiting maximum magnitude of 7.5 mag. and with periods shorter than one day. While we wanted to isolate true contact binaries, their separation from the related semi-detached and poor-thermal-contact systems was not easy. For that reason we initailly considered the EW and EB groups together, recognizing that the semi-detached EB systems are on the average brighter than contact systems so that we see them deeper in space; therefore, the EB systems are over-represented in a magnitude-limited sample such as the one considered here. As a first step, we carefully checked the light-curve types and in a few cases exchanged the EW and EB types. Since most systems do not have radial velocity data and the depth-of-eclipse criterion does not always give a unique answer, the separation of systems into the two groups remains preliminary. The “Ell” (ellipsoidal) variables were also included because they are mostly a mixture of the EW or EB systems just seen at low orbital inclination angles. We made an effort to assign them to either EW or EB groups on the basis of light-curve shapes and relative depths of eclipses – equal or unequal – respectively. The final numbers are 27 EW systems and 14 EB systems. Only 13 EW and 10 EB systems have retained the classification as in the Hipparcos catalogue. An additional complication in a separation of the two groups of binaries is the fact that some apparently genuinely contact systems show unequally deep eclipses. There are very few such poor-thermal-contact systems, about 2 percent in the volume-limited OGLE sample (Rucinski, 1997b), but they do have deeper primary eclipses than well behaving contact binaries and would normally be classified as the EB systems.

The EB systems appear to have – on the average – longer periods than EW systems so that different period distributions for both groups can help in assigning the type and separating the two groups. In particular, contact binaries are extremely rare in volume-limited samples for periods longer than about 0.6 – 0.65 day (see Figure 1 in Rucinski (1998b)). The left panel of Figure The photometric-amplitude and mass-ratio distributions of contact binary stars shows the amplitudes plotted versus the orbital period for the whole sample of 41 systems with $`P<1`$ day, with different symbols for EW and EB groups. With the additional constraint of $`P<0.65`$ days, the 7.5 magnitude-limit sample shrinks by about one half, but presumably consists mostly of contact binaries; it consists of 20 EW systems and 4 EB systems (among the latter, two with very small and two with moderate amplitudes, hence not really typical for the EB’s). Figure The photometric-amplitude and mass-ratio distributions of contact binary stars shows how (still uncertain) light contributions from close companions affect the observed amplitude distribution. While the companions are well known in frequently observed systems such as 44 Boo, HT Vir or VW Cep, not all stars of the 7.5 magnitude-limit sample have been scrutinized for presence of companions so that the amplitude distribution for the subsample with $`P<0.65`$ days, shown in the right panel of Figure The photometric-amplitude and mass-ratio distributions of contact binary stars, must be considered as preliminary.

The 24 systems with $`P<0.65`$ days form too small a sample to securely define the amplitude distribution for the mass-ratio determination. Since it is quite unlikely that the 7.5 magnitude-limit sample is missing bright contact systems with amplitudes larger than $`a0.1`$, an increase in numbers for better statistics could be now achieved by deeper systematic searches for contact sky-field systems. (Note, that according to the previous Section, the total loss of “zero-amplitude” or $`a<0.1`$ systems is still substantial at about 20 to 40 percent.) Currently, an extension beyond the 7.5 magnitude limit on the basis of catalogue data would be too risky to attempt: From among the 41 systems of the sample, as many as 17 (i.e. 41 percent) have been discovered by the Hipparcos mission which is complete only to $`V7.27.8`$. There exists no other deeper all-sky survey which would compete with the Hipparcos survey in terms of the systematic temporal coverage and photometric accuracy.

Despite the small size of the 7.5 magnitude-limit sample, we can note in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars the absence of large-amplitude systems (except for HT Vir, but only after its amplitude is corrected for the presence of its companion) and the well-defined rise of the amplitude distribution toward the small amplitudes, as expected by our results in Section 3.

### 4.2 Observational data: The OGLE sample

The results of the OGLE project (Rucinski, 1997a, b), more fully interpreted in Rucinski (1998b), provide sound data on the observed amplitude distribution, $`A_{obs}(a)`$. The statistics are based on two volume-limited samples, to $`d=3`$ kpc and to $`d=5`$ kpc, designated as BW3 and BW5. As discussed in Rucinski (1998b), the samples are complete to the absolute magnitudes $`M_V=5.5`$ and $`M_V=4.5`$, respectively. The sample sizes are 98 systems for BW3 and 238 systems for BW5, with BW5 including BW3. There may exist a dependence of the amplitude on the absolute magnitude: The hotter, brighter systems may show a stronger admixture of EB binaries. Since BW5 consists preferentially of brighter systems seen deeper in space, while BW3 is based on fainter, more local systems, it was felt prudent to consider the two samples separately. The BW3 sample would be in general the preferred one as it is expected to better represent the typical population of contact systems.

The main limitation of the statistics based on the OGLE data is blending of the images in the crowded Baade’s Window area, leading to systematically smaller variability amplitudes. The random-pairing blending occurs on top of the influence of close companions, as for the nearby stars. The difficulty is that in the case of the OGLE survey the stars are not easily accessible to medium-resolution spectroscopy which would not only provide confirmation of binarity (i.e. elimination of $`\delta `$ Sct and RR Lyr stars), but would also permit detection of spectroscopic companions. Since we cannot address the matter of blending, we ignore it entirely. We suggest that our analysis of the OGLE data should be taken as an illustration how blending-corrected data could normally be treated using our approach.

The observed amplitude distributions $`A(a)`$ for BW3 and BW5 are tabulated in Table 3 and are shown in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars by shaded histograms. The distributions are given with the amplitude bins of $`\mathrm{\Delta }a=0.05`$, centered on the values given in the first column of Table 3. The amplitudes are in the photometric $`I`$-band. Because the $`I`$-band amplitudes are typically only 3–5 percent shallower than in the $`V`$-band, we disregarded the small difference which is immaterial in view of the current, poor statistics. However, the matter of the band matching may have to be addressed in future, by more advanced investigations.

The completeness threshold for discovery of the OGLE sample was estimated in Rucinski (1997a) at about 0.3 mag. As stated in Section 3, the basis for this estimate was the absence of a further increase in numbers of detected systems for $`a<0.3`$ mag. The assumption that the OGLE sample is complete for $`a>0.3`$ may be conservative, but provides full assurance of an unbiased statistics of the amplitudes. Also, even when measurement errors ($`\sigma `$) are at the level of 0.01 – 0.03 mag., as for the OGLE project, detection of variability requires a signal several times larger, say $`5\sigma `$. To characterize the variability and estimate the variability type requires still more margin. All in all, the full completeness limit of 0.3 mag. is not at all unrealistic in such a case. The OGLE project has in fact discovered several contact binaries with $`0.1<a<0.3`$, but we suspect that not all contact systems have been discovered in this interval.

The continuous lines in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars give the predicted amplitude distribution $`A(a)`$ calculated assuming $`f=0.25`$ and a flat distribution $`Q(q)=const`$ (this is the same as the one marked by the thin line in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars), superimposed on the OGLE distributions. Absence of large-amplitude systems in $`A_{obs}(a)`$ is striking. This may be because the mass-ratio values close to unity are extremely rare ($`Q(q)`$ rising for $`q0`$) or because of strong blending of images or – most likely – both. Since we cannot estimate the blending effects, we assume that the shape of $`Q(q)`$ is reflected in $`A(a)`$. The results of this assumption are described below.

### 4.3 Determination of the mass-ratio distribution

An observed amplitude distribution can be modeled by adjusting the $`Q(q)`$ distribution. We can predict the $`A(a)`$ distribution by utilizing the computed distributions $`A_q(a)`$ (Section 3):

$$A_{pred}(a)=\mathrm{\Sigma }Q(q_i)A_{q_i}(a)$$
(2)
Strictly speaking, we should denote the fact that we use a specific value of $`f`$, so that a proper superscript would be in order. For clarity, in what follows, we will assume $`f=0.25`$ unless noted otherwise. The functions $`A_{q_i}(a)`$ are each normalized to unity (for each interval of $`q`$); this permits expressing $`A_{pred}(a)`$ and $`Q(q)`$ in the actual numbers of systems so that uncertainties can be simply estimated from the Poisson statistics. We attempted to determine $`Q(q)`$ by representing it by five independent bins $`\mathrm{\Delta }q=0.2`$ wide and by two-parameter power-laws, as described below.

#### 4.3.1 Random-search fits

As the first step in estimating $`Q(q)`$ on the basis of the OGLE amplitude distribution, a simple random search for a best fit of $`A_{pred}`$ to $`A_{obs}`$ was conducted by using five bins of $`\mathrm{\Delta }q=0.2`$. Each bin of $`Q(q_i)`$ was considered independent, without any assumption of smoothness or continuity constraints on $`Q(q)`$. The solution was obtained by an extensive random trial search, iterated until the smallest value of the “quality-of-fit” measure, $`\chi ^2`$, defined as $`\chi ^2=(A_{obs}A_{pred})^2/\sigma _A^2`$ was found. Poissonian estimates $`\sigma _A=\sqrt{A_{obs}}`$ for each bin were used for the standard errors. Because of the small number of filled bins in $`A_{obs}`$ for $`a>0.3`$ (7 and 8 for BW3 and BW5), the 5-parameter description of $`Q(q)`$ obviously could be considered only as indicative, yet perhaps useful as the first stage.

The detailed results on $`Q(q)`$, expressed as the number of systems per a bin of $`\mathrm{\Delta }q=0.2`$, are given numerically and graphically (continuous line histogram) in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars. Note that the first bin $`0q0.2`$ senses the amplitude distribution only between our low limit of $`a=0.3`$ and $`a0.43`$ which corresponds to $`q=0.2`$. Nevertheless, this is the region where the observational $`A(a)`$ for the OGLE samples were best determined.

While the main solution was done with $`A_q^{0.25}`$, i.e. for the case of $`f=0.25`$, as tabulated in Table 1, we made also fits for $`f=0`$ and $`f=0.5`$ ($`A_q^{0.0}`$ and $`A_q^{0.5}`$). These solutions are shown in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars by dotted and broken line histograms. They were important in establishing the sensitivity of the results to the presently poorly constrained value of $`f`$. Because of the large number of the free parameters (5 bins in $`Q(q)`$) relative to the number of independent data (7 or 8 bins in $`A(a)`$ for both BW samples), the individual uncertainties for each bin of $`Q(q)`$ were very large, in fact much larger than the Poissonian uncertainties. Still, it did not prevent us from iterating the random search to a unique and stable solution for each value of $`f`$. All solutions turned out to be very similar for all three values of the fill-out parameter, in spite of the very poorly constrained search. This leads us to conclude that, irrespectively of the assumed value of $`f`$, the mass-ratio distribution appears to steeply rise for very small values of $`q`$; at the present level of accuracy, the matter of the degree-of-contact is of secondary importance. For that reason, the subsequent analysis will consider only the case of $`f=0.25`$.

#### 4.3.2 Power law approximations

The next assumption is that the mass-ratio distribution can be represented by a power law. Here we have a choice of expression, either $`Q_a(q)=a_0(1q)^{a_1}`$, with $`a_1>0`$ or $`Q_b(q)=b_0q^{b_1}`$, with $`b_1<0`$. The first form is a bit more convenient because $`a_0`$ multiplies a factor which is confined between 0 and 1, so that one has a clean separation of the shape dependence, controlled by $`a_1`$, from the normalization, controlled by $`a_0`$. The latter form is preferable for comparison with theory which usually involves straight power-law expressions in $`q`$. The function $`Q_b`$ tends to infinity for $`q0`$; however, there exists a low limit to $`q`$ that prevents the divergence, as we discuss in Section 5. Both forms involve only two parameters. We found that we cannot, at the present time, generalize them by addition of an additive parameter. Tests of the significance of such a third parameter indicate an insignificant decrease of $`\chi ^2`$ so that two-parameter fits must currently suffice. An absence of the additive term in $`Q(q)`$ agrees with the expectation that $`Q(q)0`$ for $`q1`$, an effect which is caused by a thermal instability at $`q=1`$ (Lucy, 1976; Flannery, 1976); we discuss this further in Section 6. The calculations of the predicted $`A(a)`$ were made with 10 bins of $`\mathrm{\Delta }q=0.1`$ with the distributions $`A_q(a)`$ as given in Table 1.

The solutions for $`Q_a`$ and $`Q_b`$ are given in Table 4. In terms of the quality of fit $`\chi ^2`$, the solutions for $`Q_a`$ are slightly better than those for $`Q_b`$. Figure The photometric-amplitude and mass-ratio distributions of contact binary stars shows the $`\chi ^2`$ contours corresponding to the one-sigma errors of both parameters for both solutions based on admissible range of $`\mathrm{\Delta }\chi ^2=2.3`$ above $`\chi _{min}^2`$ for the 68 percent significance level. Because of the very small number of data in the observational $`A(a)`$ and of the correlation between the multiplicative and power parameters, the errors of the parameters are large. At this point, we are unable to decide which power law is the correct one.

Figure The photometric-amplitude and mass-ratio distributions of contact binary stars illustrates the $`Q(q)`$ solutions for both power laws (left panels) and the implied amplitude distributions (right panels). The continuous and broken lines show the best fitting $`Q_a(q)`$ and $`Q_b(q)`$ distributions and the resulting $`A(a)`$ distributions. Qualitatively, the two forms of $`Q(q)`$ appear to be similar in the range $`0.12q1`$ which maps into $`A(a>0.3)`$. However, the two power laws do differ at the quantitative level when specific predictions for populations of individual bins are compared. For example, we can follow Van’t Veer (1978), who compared the ratios of the populations in the bins $`0.1<q<0.2`$ and $`0.8<q<0.9`$ (see below in Section 6). We find that the better fitting law $`Q_a`$ predicts the ratio of populations of these bins of $`20\times 10^3`$ while $`Q_b`$ predicts the ratio of 37. Unfortunately, with only 50 systems spread over 7 amplitude bins for BW3 and with 120 systems spread over 8 amplitude bins for BW5, the distributions $`A(a)`$ cannot constrain the results any better.

For further considerations we will simplify the results to $`Q_a(q)(1q)^{a_1}`$, with $`a_1=6\pm 2`$, and $`Q_b(q)q^{b_1}`$, with $`b_1=2\pm 0.5`$. This simplification is justified in view of (1) the large parameter errors for both power laws so that our results are only very preliminary, (2) the difference in the results for the BW3 and BW5 samples, and (3) our preference for the BW3 sample which is a very small, yet is more rigorously defined.

### 4.4 Expected discovery selection at low amplitudes

Having the predictions of the amplitude distributions which best fit the amplitude range $`a>0.3`$, we can check how many systems would be predicted over the whole range of amplitudes, including those below this amplitude limit. These estimates depend very strongly on the shape of $`Q(q)`$ for small mass ratios. We give the predictions for the power-law approximations of $`Q(q)`$ in Table 4 in the line $`\mathrm{\Sigma }A_{pred}`$. The numbers are of all systems expected over the whole range of amplitudes. For the OGLE sample, we can compare them with the actual numbers, including the systems with small amplitudes below $`a=0.3`$. The number of observed systems is 98 and 238, for BW3 and BE5, respectively. For the $`Q_a(q)`$ distribution, the ratio $`N_{pred}/N_{obs}=3.1`$ is identical for both BW samples, but for the $`Q_b(q)`$ distribution, the ratio is 8.8 and 12.9. Obviously, to a large degree, these estimates measure the amount of divergence of $`A_{pred}(a)`$ for $`a<0.3`$. They cannot be used to address the important issue of the conversion of the apparent frequency to actual (spatial) frequency of contact binaries. There is one effect which prevents the conversion factor from becoming uncomfortably high. It is the low limit on the mass ratio which is described in the next section.

## 5 THE MINIMUM MASS RATIO

Webbink (1976) pointed out that stability of a contact binary is compromised by a redistribution of angular momentum for very small values of the mass ratio. For a mass ratio smaller than a threshold value $`q_{min}`$, the system will find it easier to store its angular momentum in one star rather than in an extreme mass-ratio binary so it will quickly (in a dynamical time scale comparable to one orbital revolution) merge into a single, rapidly rotating star. Rasio (1995) re-analyzed this tidal instability process and lifted the expected number from the very small value suggested by Webbink (1976) to $`q_{min}0.09`$. The exact location of the limit may depend on the evolutionary state of the stars. This can in fact explain the existence of such a well-known system as AW UMa with $`q=0.075`$.

The distributions $`A_q(a)`$ sampled at $`\mathrm{\Delta }q=0.1`$, as in Table 1 and Figure The photometric-amplitude and mass-ratio distributions of contact binary stars, are obviously useless in predicting $`A_{pred}`$ in the presence of the cut-off at $`0.07q_{min}0.09`$. However, we can use these distributions to obtain a preliminary (upper limit) estimate on $`A_{pred}`$ by simply setting $`Q(q)=0`$ for the first bin $`0<q<0.1`$. The predicted number of all systems is then substantially reduced, for the $`Q_a`$ law from 301 and 732 systems to 143 and 350 and, for $`Q_b`$ law from 850 and 3082 to 147 and 359, for BW3 and BW5, respectively. Thus, the ratio $`N_{pred}/N_{obs}`$, which was as large as 3 to 12 with the full $`Q(q)`$ distributions discussed in the previous section, is reduced to more acceptable levels: $`N_{pred}/N_{obs}=1.46`$ and 1.47 for $`Q_a(q)`$ and $`N_{pred}/N_{obs}=1.50`$ and 1.51 for $`Q_b(q)`$. Similarity of these numbers at about the level of about 1.5 attests to the fact that the power-laws $`Q_a`$ and $`Q_b`$ are both equally successful when only the amplitudes $`a>0.3`$ and the mass ratios $`q>0.1`$ are considered.

The above estimate is approximate because the cut-off is almost certainly below $`q_{min}=0.1`$. To reproduce the shape of $`Q(q)`$ below this point for use in Equation 2, we must resort to the original fine grid of $`A_q(a)`$ calculated with small bins of $`\mathrm{\Delta }q=0.01`$. The results for the specific case of the BW3 sample and for both power laws are shown in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars. Note the particularly strong influence of the cut-off at $`q_{min}`$ on the $`Q_b`$ results when the very steep increase in $`A(a)`$, due to the divergence of $`Q_b`$ for $`q0`$, is avoided.

One can calculate integrals of the curves in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars and take their ratio. Values of the ratio $`N_{pred}/N_{obs}`$ for $`q_{min}=0.10`$, 0.08 and 0.06 (and for an imaginary extension down to $`q=0.01`$) are given in Table 5. $`N_{pred}/N_{obs}`$ is the correction factor which can be used in converting the apparent (inclination uncorrected) frequency of occurrence of W UMa systems – as derived on the basis of the BW3 sample – into the true spatial frequency of occurrence. The uncertainty with the value of $`q_{min}`$ prevents us from establishing this factor to any better than about 1.5 to 2.0. Thus, the currently best estimate of the inclination-uncorrected frequency of contact binaries in the old disk population of about 1/130 relative to FGK dwarfs (Rucinski, 1998b) would translate into the spatial frequency of about 1/80 – 1/65.

The last two columns of Table 5 give the predicted ratio of the total number of systems to those detectable by surveys fully complete for amplitudes $`a>0.1`$, for both power laws. The missed fraction for such surveys would be very similar, with the correction factor at the level of 1.4 to 1.5 times. The closeness of the estimates of $`N_{pred}/N_{obs}`$ for the BW3 sample (about $`1.52.0`$) and for a fully complete sample down to $`a=0.1`$ (about $`1.41.5`$) is due to (1) a flattening of $`A(a)`$ in the region $`0.1<a<0.3`$ caused by the mass-ratio cut-off, and (2) the fact that the OGLE sample, although probably not complete below $`a=0.3`$, contains a fair number of systems with $`0.1<a<0.3`$.

The shape of a well defined amplitude distribution $`A(a)`$ down to $`a=0.1`$ is expected to sensitively reflect the location of $`q_{min}`$. We can see in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars that a local maximum is expected to form in $`A(a)`$ around $`a0.20.25`$. This maximum is better defined for the $`Q_b(q)`$ power law because the rise of the predicted $`A(a)`$ below $`a=0.3`$ is steeper for this law so that the effects of the cut-off in $`q`$ are stronger. We actually see a maximum in the observed $`A(a)`$ for the OGLE sample exactly in this interval, but we suspect that this feature is simply due to the detection selection effect setting in for $`a<0.3`$. To be sure of the existence of the local maximum, we should see indications of a small minimum beyond the peak and of the further increase in $`A(a)`$ for $`a<0.150.2`$. As we have said above, an extension of completeness down to $`a=0.1`$ will still leave some 40 to 50 percent of all systems below the detection level.

## 6 THE MASS-RATIO DISTRIBUTION IN THE CONTEXT OF CONTACT BINARY EVOLUTION

The mass-ratio distribution for W UMa-type systems and the evolution of this distribution over time are related to entirely different processes than those of star formation producing an almost flat $`Q(q)`$ for detached binaries (Mazeh et al., 1992). Contact binaries have the freedom of exchanging mass between components through a complex interplay of energy exchange, mass exchange and angular-momentum loss (AML) processes. The first major re-structuring takes place at the moment when the contact system forms from two detached components; from that point, further, more gradual changes in the mass distribution are expected as the system evolves over time.

Following the pioneering works of Lucy (1976) and Flannery (1976) who showed that contact systems are inherently thermally unstable and will evolve away from $`q=1`$ to small mass ratios, several theoretical models explored in detail the thermal-relaxation oscillations and ways of preventing them, either through nuclear evolution or through AML, or perhaps combination of both processes (Robertson & Eggleton, 1977; Rahunen, 1981, 1982, 1983). At this moment, the unified scenario of the contact binary formation and evolution presented by Vilhu (1982), appears to be still valid. Among its unexplained features, the most mysterious remains a coupling and/or feedback process between the degree of contact and the amount of AML which prevents rapid coalescence on one hand and disruption of contact (a semi-detached phase) on the other hand. Our understanding of these processes crucially depends on the estimates of the relevant evolutionary time scales which can be estimated from the statistics of $`Q(q)`$.

The discussion of $`Q(q)`$ in Vilhu (1981) revolved around the (then) only available observational derivation of Van’t Veer (1978), based on very meager, combined photometric and radial-velocity data. These results most probably contained strong discovery and observer preference selection effects, in spite of heroic attempts to estimate their size. This distribution was much less steep than the power laws derived in the current paper: Van’t Veer (1978) estimated that his distribution implies ten times more systems in the $`0.1<q<0.2`$ bin than in the $`0.8<q<0.9`$ bin. Our power laws lead to a much larger disparity in the population of these bins: The ratio predicted by $`Q_a(1q)^6`$ is $`20\times 10^3`$ times (there would be almost no large mass-ratio systems), while $`Q_bq^2`$ predicts the ratio of 37 times. While the large difference in the predictions will eventually help in selecting the correct shape for $`Q(q)`$, we are not at present in a position to prefer one power law over the other as both give very similar fits to the observed $`A(a)`$ for $`a>0.3`$.

As discussed by Vilhu (1981) (see Figure 4 in this paper), the mass-ratio evolution in contact is driven mainly by the less-massive component and its thermal (Kelvin–Helmholtz) time scale. When the evolution reaches a steady state condition, the number of systems in a particular evolutionary stage should scale as $`N\tau _{sec}M_{sec}^{1\beta }`$, where $`\beta `$ is the exponent in the mass–luminosity relation, $`LM^\beta `$. For the lower main sequence, $`\beta 4.5`$, so that $`\tau _{sec}M_{sec}^{3.5}`$. Since $`M_{sec}=M_{tot}q/(1+q)`$, for small values of $`q`$, we can expect $`Q(q)q^{3.5}`$. Thus, as the secondary components become less massive, their evolutionary time scale becomes progressively longer resulting in a pile up of contact systems at low mass ratios. This pile-up is limited by the tidal instability at $`q_{min}`$, as discussed in Section 5.

The theoretical expectations described above very well agree with our results which we present in a schematic form in Figure The photometric-amplitude and mass-ratio distributions of contact binary stars. We note that location of the tidal instability at $`q_{min}0.070.1`$ is a relatively new development (Rasio, 1995) and could not be included in the general discussion of Vilhu (1981). Its presence is actually crucial in preventing the problem of embarrassingly too many small mass-ratio systems, if the $`Q(q)q^{3.5}`$ distribution were to continue below $`q0.070.1`$.

Finally, a cautionary note: If the evolution were really to slow down as the mass of the secondary component decreases, i.e. as $`\tau _{sec}M_{sec}^{3.5}`$ then, for very small $`M_{sec}`$, it would take longer than the Hubble time. Assuming the thermal time scale for the Sun, $`\tau _{}3\times 10^7`$ years, then for $`M_{sec}=0.1M_{}`$, $`\tau _{sec}10^{11}`$ years. This would lead to a very inefficient, practically insignificant formation of single, rapidly-rotating stars from contact binaries. But – more likely – the nuclear or AML evolution of primary components, with the associated shorter time-scales, will become more important first. Thus, we have no idea about the rate of production of single stars at the cut-off at $`q_{min}`$, but it may be not as low as the thermal evolution of secondary components would imply.

## 7 CONCUSIONS

This paper discusses the expected amplitude distribution $`A(a)`$ for contact binary stars. The strong dependence of $`A(a)`$ on the mass-ratio distribution, $`Q(q)`$, has been shown to be useful for shedding light on the latter distribution which has a considerable astrophysical significance. We attempted to simplify the details of the approach and to concentrate on the main properties of both distributions. In particular, while the results do depend on the degree-of-contact parameter, $`f`$, the dependence is weak and – for simplicity – most of the discussion has been presented for the most likely value of $`f=0.25`$.

The main limitation of the paper are problems created by presence of “third light” in photometry of contact binaries. Both, presence of unresolved visual companions and of blending in crowded areas such as Baade’s Window, are expected to produce distorted amplitude distributions with a diminished representation of systems with large amplitudes. The degree of such a distortion is difficult to quantify, primarily because the unknown frequency of occurrence of companions to contact systems; this frequency may be different than for wider binaries. The data for the 7.5 magnitude-limit sky-field sample indicate that the main conclusions of the paper are valid even after accounting for presence of close companions.

The main results of the paper are summarized below with references to Figure The photometric-amplitude and mass-ratio distributions of contact binary stars.

1. The two distributions, of the mass ratio, $`Q(q)`$, and of the photometric variability amplitude, $`A(a)`$, are very closely inter-related. Since the $`Q(q)`$ distribution is expected to contain a record of the contact binary evolution, studies of $`A(a)`$ can help in resolving the still poorly understood details of time scales of formation and evolution of contact systems.
2. The amplitude distribution $`A(a)`$ is expected to rise for small amplitudes almost irrespectively of the shape of $`Q(q)`$. This rise is expected for a flat distribution, $`Q(q)=const`$, or even for any of imaginary “monochromatic” distributions with all systems having just one mass ratio, $`Q(q)=\delta (qq_0)`$. The rise becomes even stronger if $`Q(q)`$ steeply increases for small $`q`$, as it appears to be the case.
3. The increase of $`A(a)`$ for $`a0`$ continues to zero amplitude and leads to a convergence to a constant (mass-ratio dependent) value: $`A_q(a0)C(q)`$.
4. Two samples of contact binaries have been considered: The sample of bright systems to 7.5 magnitude and the sample of systems discovered in the Baade’s Window by the OGLE project. While the former appears to be complete to $`a0.1`$ and has been corrected for the presence of currently known companions, it is too small for derivation of $`Q(q)`$ from $`A(a)`$. The latter sample is marginally sufficient in size (98 or 238 systems depending on the spatial depth) and gives a moderately well-defined amplitude distribution, but is only complete for $`a>0.3`$ and is certainly subject to the influence of blending problems which tend to depress the large amplitude end of $`A(a)`$. The conclusions below are preliminary on account of the neglected photometric blending for the OGLE sample.
5. Thanks to the non-linearity of the relation between $`A`$ and $`Q`$, an amplitude distribution complete for $`a>0.3`$, maps into $`Q(q)`$ within $`0.12q1`$, so that the accessible range of $`q`$ in $`Q(q)`$ is respectable. Figure The photometric-amplitude and mass-ratio distributions of contact binary stars shows our best power-law estimates of $`Q(q)`$.
6. The $`Q(q)`$ distribution derived from the OGLE distribution $`A(a)`$ for the interval $`0.12q1`$ climbs very steeply for $`q0`$; it can be approximated by $`(1q)^6`$ or $`q^2`$. The values of the exponents in both expressions are very preliminary, not only because of large statistical errors but – more importantly – because of the distortions to $`A(a)`$ introduced by the blending effects.
7. The steep increase of $`Q(q)`$ is expected to be abruptly terminated at $`q_{min}0.070.1`$ by the process of tidal instability. This alleviates the danger of huge numbers of very small amplitude systems which would be hiding below the detection thresholds, if the approximate power-law relationships were to continue to $`q0`$.
8. The most common contact systems in the interval between the cut-off at $`q_{min}0.070.10`$ and the steep power-law drop at $`q0.30.4`$ are expected to generate a local maximum in the amplitude distribution in the vicinity of $`a0.200.25`$. The exact location of this maximum and the rate of increase of $`A(a)`$ for $`a0`$ will help to establish the value of $`q_{min}`$ which is currently poorly established.
9. It is expected that future well-determined amplitude distributions, good down at least to $`a0.1`$ and with fully characterized blending will define the exact shape of $`Q(q)`$ in the vicinity of the cutoff at $`q_{min}`$. The current sky-field sample to $`V=7.5`$ contains too few systems for a good definition of $`A(a)`$ at small amplitudes; a deeper sample is needed. A complete sky-field sample has a potential of a better control over the problem of unresolved companions than the OGLE sample because of the accessibility to spectroscopic studies.
10. The previous analysis of the OGLE sample led to an estimate of the inclination uncorrected frequency of contact binaries of about 1/130 relative to FGK dwarfs. The OGLE sample contains contact systems with the smallest amplitudes of about 0.1 mag. and appears to be fully complete for $`a>0.3`$ mag. At present, we estimate that a correction factor to convert the OGLE apparent frequency into the spatial frequency is about 1.5 to 2.0, but the exact value sensitively depends on the value of $`q_{min}`$. Thus, the inclination-corrected spatial frequency is one contact binary per 1/80 to 1/65 Disk Population FGK dwarfs.

I would like to thank Stefan Mochnacki for useful discussions about various subjects related to statistics of contact binary stars, to Greg Stachowski for careful reading of the manuscript and for linguistic corrections and to Janusz Kałużny who – acting as a referee – provided several useful suggestions for improvement of the final manuscript. The author acknowledges with gratitude the fact that this paper utilizes the data obtained and made available for public access by the OGLE project. Support of the Natural Sciences and Engineering Council of Canada is acknowledged with gratitude.

Captions to figures:


format: text

text:
# Anomalous finite-size effect in superconducting Josephson junction arrays

## Abstract

We report large-scale simulations of the resistively-shunted Josephson junction array in strip geometry. As the strip width increases, the voltage first decreases following the dynamic scaling ansatz proposed by Minnhagen et al. \[Phys. Rev. Lett. 74, 3672 (1995)\], and then rises towards the asymptotic value predicted by Ambegaokar et al. \[Phys. Rev. Lett. 40, 783 (1978)\]. The nonmonotonic size-dependence is attributed to shortened life time of free vortices in narrow strips, and points to the danger of single-scale analysis applied to a charge-neutral superfluid state.

Dissipation of the supercurrent through superfluid and superconducting films is caused by the flow of free vortices transverse to the current. This phenomenon was analyzed in detail by Ambegaokar, Halperin, Nelson, and Siggia (AHNS). In the context of a uniform superconducting film in zero magnetic field, their theory predicts a power-law current-voltage ($`I`$-$`V`$) relationship with a temperature-dependent exponent. Despite the relatively simple construct of the theory, its experimental and numerical verification has remained controversial . It has been noted that boundary and finite-size effects can dominate the measured voltage drop across the sample at sufficiently low temperatures, making it difficult to achieve an unambiguous comparison between theory and experiments.

An alternative theory of vortex flow dissipation which also yields nonlinear $`I`$-$`V`$ curves in the superconducting state derives from the dynamic scaling hypothesis. For superconducting films and Josephson junction arrays (JJA), the two approaches differ in their predictions of the exponent $`a(T)`$ characterizing the power-law behavior $`VI^a`$ below the Kosterlitz-Thouless-Berezinskii (KTB) transition temperature $`T_{\mathrm{KT}}`$. According to the AHNS theory,

$$a_{\mathrm{AHNS}}(T)=\chi +1,$$
(1)
where $`\chi =\pi J_R/k_BT`$ and $`J_R`$ is the renormalized spin-wave stiffness. On the other hand, the dynamic scaling analysis of Minnhagen et al. suggests

$$a_{\mathrm{MWJO}}(T)=2\chi 1.$$
(2)
The two expressions coincide at $`T_{\mathrm{KT}}`$ where $`\chi =2`$, but the difference grows rapidly as one moves to low temperatures. While there seem to be ample numerical support for the Minnhagen et al. scaling, agreement with the AHNS theory has also been reported. To rationalize the two scenarios, Bormann made an interesting suggestion that strong current creates a dense set of vortex-antivortex pairs, thereby invalidating the AHNS treatment. Nevertheless, the AHNS theory should still apply at sufficiently weak currents. This picture has not been borne out by a recent numerical study which shows persistently larger value of $`a`$ than Eq. (1) predicts.

In this Letter, we present $`I`$-$`V`$ results from large-scale simulations of the resistively shunted Josephson-junction (RSJ) array. A rectangular strip geometry is used to decouple boundary effects introduced by the current leads from finite-size effects arising from vortex/antivortex motion transverse to the current $`I`$. The main finding of our work is the existence of three distinct regimes as the strip width $`L_y`$ is varied. For $`L_y<l_\mathrm{b}I^1`$, one is in the linear-response regime where the near-equilibrium dynamic scaling hypothesis applies. In this regime, the voltage $`V`$ decreases as $`L_y`$ increases, and reaches a minimum value predicted by Minnhagen et al. When $`L_y`$ exceeds a second length scale $`l_\mathrm{r}I^{\chi /2}`$, $`V`$ saturates to the asymptotic value predicted by the AHNS theory. The intermediate regime $`l_\mathrm{b}<L_y<l_\mathrm{r}`$ is characterized by a rising $`V`$ as $`L_y`$ increases. This behavior is explained in terms of “self-recombination” of vortex-antivortex pairs under periodic boundary conditions.

We begin with the model Hamiltonian of a two-dimensional (2D) JJA in zero magnetic field,

$$H=J\underset{ij}{}\mathrm{cos}(\varphi _i\varphi _j).$$
(3)
Here $`\varphi _i`$ denotes the phase of the superconducting order parameter on grain $`i`$, and $`J`$ sets the strength of the Josephson coupling between neighbouring superconducting grains. The summation in Eq. (3) is over all nearest neighbor grain pairs. In the RSJ dynamics, the network is insulated from the substrate, but normal currents are allowed between neighboring grains, in addition to the supercurrent. The dynamics of the phase variables follows the Josephson relations plus the Kirchoff law for current conservation at each node,

$$\frac{\sigma \mathrm{}}{2e}\underset{j\mathrm{n}.\mathrm{n}.\mathrm{of}i}{}\left(\frac{d\varphi _i}{dt}\frac{d\varphi _j}{dt}\right)=\frac{2e}{\mathrm{}}\frac{H}{\varphi _i}+I_{\mathrm{ext}}.$$
(4)
Here $`\sigma `$ is the conductivity of normal links, and $`I_{\mathrm{ext}}`$ is the external current source which is nonvanishing only at the boundaries of the network. To model the effect of temperature $`T`$, a gaussian noise current $`\eta _{ij}`$ is added to each link between two adjacent nodes $`i`$ and $`j`$, with $`\eta _{ij}(t)\eta _{i^{}j^{}}(t^{})=2\sigma k_BT\delta _{ij,i^{}j^{}}\delta (tt^{})`$.

To integrate the coupled equations (4), one needs to first solve for $`V_i=(\mathrm{}/2e)d\varphi _i/dt`$. This can be done efficiently using a pseudo-spectral algorithm. In our simulations, the network is chosen to be a strip of $`L_x`$ columns, with $`L_y`$ nodes in each column. Open boundary conditions are used in the $`x`$-direction (along the strip) to allow for current inputs and outputs. Periodic boundary conditions (PBC) are used in the $`y`$-direction perpendicular to the strip. After performing the Fourier transform in the $`y`$-direction, the left-hand-side of Eq. (4) reduces to a tridiagonal form which can be easily solved by the Gauss elimination method. The time stepping is done using a second-order Runge-Kutta method with $`\mathrm{\Delta }t=0.1`$ (in units of $`\sigma \mathrm{}/2eI_c`$, where $`I_c=2eJ/\mathrm{}`$ is the critical current through each link). To minimize the boundary effects, voltage drop across the sample is measured in the interior of the network. For the data presented below, we typically choose $`L_x=2048+2\times 128`$ and discard 128 columns at each end for this purpose.

When the external current $`I`$ into each boundary node is set to zero, we recover a dynamical version of the XY model whose equilibrium properties have been well characterized. In particular, the KTB transition temperature has been determined to be $`T_{\mathrm{KT}}0.89J/k_B`$. We have confirmed this value in the simulation by examining the correlation function of the superconducting order parameter and the helicity modulus (or $`J_R`$), which exhibit universal behavior at the transition. Direct measurement of the $`I`$-$`V`$ curve at this temperature has also confirmed the predicted exponent $`a(T_{\mathrm{KT}})=3`$.

We now describe the main results of our simulation. Figures 1(a) and 1(b) present the $`I`$-$`V`$ data at $`T=0.8J/k_B`$ and $`T=0.7J/k_B`$, respectively, for four different array sizes: $`L_y=8,32,128`$, and 512. Here the current $`I`$ is measured in units of $`I_c`$, and $`V`$, the voltage drop per column, in units of $`\mathrm{}/2e`$. As the current $`I`$ decreases, the finite size effect becomes more and more evident. Examining the data at the smallest $`I`$ shown, we see a clear nonmonotonic size dependence: for small values of $`L_y`$, $`V`$ decreases with increasing $`L_y`$. However, this trend is reversed as $`L_y`$ is increased further. The solid and dashed lines in the figure correspond to the power-law behavior predicted by AHNS \[Eq. (1)\] and Minnhagen et al. \[Eq. (2)\], respectively, using $`\chi =2.8`$ at $`T=0.8J/k_B`$ and $`\chi =3.5`$ at $`T=0.7J/k_B`$. These values of $`\chi `$ are obtained from measurements of the helicity modulus at the respective temperatures, and agree with previous studies. For $`T=0.8J/k_B`$, the $`I`$-$`V`$ data at $`L_y=512`$ is well-described by the AHNS theory. On the other hand, the envelop of data sets at smaller values of $`L_y`$ follows the Minnhagen et al. formula (2). For $`T=0.7J/k_B`$, we are only able to obtain reliable data for $`I0.1I_c`$ due to the rapid decrease of $`V`$ with decreasing $`I`$. Nevertheless, better agreement with the AHNS prediction is apparent on the low current side for the $`L_y=512`$ array. The larger apparent slope of the data at higher currents may be attributed to the larger value of $`J_{\mathrm{eff}}`$ due to the smaller length scale probed, which yields a higher $`\chi _{\mathrm{eff}}=\pi J_{\mathrm{eff}}/k_BT`$.

Figure 2 shows three sets of voltage data against transverse array size $`L_y`$. Data at $`T=0.8J/k_B`$ are represented by open ($`I=0.1I_c`$) and solid circles ($`I=0.05I_c`$), while the set at $`T=0.7J/k_B`$ and $`I=0.1I_c`$ are presented by open squares. Statistical error of each data point is smaller or comparable to the symbol size. As $`I`$ decreases, the size $`l_\mathrm{m}`$ where the minimum value $`V=V_m`$ is reached increases. $`l_\mathrm{m}`$ also varies with $`T`$, though this dependence may be weak. Since the change of $`V`$ with $`L_y`$ is slow in the neighborhood of $`l_\mathrm{m}`$, it may give the false impression that the finite-size effect is already saturated at $`l_\mathrm{m}`$. As our data shows, the true asymptopic value of $`V`$ is reached only at much larger values of $`L_y`$. At $`T=0.7J/k_B`$ and $`I=0.1I_c`$, for example, the asymptopic value is larger than $`V_m`$ by five-fold.

In the remaining part of the paper, we propose an explanation of the unusual finite-size effect, based on the standard picture of creation and annhilation of free vortices under an applied current $`I`$. For small $`I`$, the supercurrent through the network is maintained by an average phase jump $`\mathrm{\Delta }\varphi =I/I_c`$ per column. Passage of a free vortex across the strip (i.e., in the $`y`$-direction) decreases the overall phase advance along the strip by an amount $`2\pi `$. In a steady-state situation, this loss is compensated by a difference in the phase velocities at the two opposing ends of the strip, and hence a voltage drop across the system. Simple counting then yields the following expression for the average voltage drop per column in a strip of length $`L_x1`$ (see, e.g. Ref.),

$$V=\frac{2\pi }{L_xL_y}\left(N_+\overline{v}_+N_{}\overline{v}_{}\right).$$
(5)
Here $`N_\pm `$ are the number of $`\pm `$ vortices in the system, and $`\overline{v}_\pm `$ are their respective mean velocities in the $`y`$-direction. Although Eq. (5) does not discriminate between bound and free vortices, it is clear that correlated drift of a vortex-antivortex pair does not contribute to $`V`$.

In the AHNS theory, the energy of a vortex/antivortex pair of size $`R`$ and oriented perpendicular to the external current is given by,

$$E(R)=2E_c+2\pi J\left(\mathrm{ln}R\frac{I}{I_c}R\right),$$
(6)
where $`E_c`$ is the core energy of a single vortex. Equation (6) has a maximum $`E=E_B=2E_c+2\pi J\mathrm{ln}(I_c/I)2\pi J`$ at $`R=l_\mathrm{b}=I_c/I`$. This is the activation energy for pair breaking. The rate for such processes follows the Arrhenius law,

$$\mathrm{\Gamma }e^{E_B/k_BT}(I/I_c)^{2\chi }.$$
(7)
In the second step of the argument, one writes down a mean-field rate equation governing the density $`\rho =\rho _\pm `$ of free vortices,

$$\dot{\rho }=\mathrm{\Gamma }\gamma \rho ^2,$$
(8)
where $`\gamma `$ is the diffusion constant. The steady-state vortex density is then given by $`\rho _{\mathrm{}}=(\mathrm{\Gamma }/\gamma )^{1/2}`$. Going back to Eq. (5), one finds

$$V_{\mathrm{}}\rho _{\mathrm{}}\overline{v}_\pm I^{\chi +1},$$
(9)
where we have assumed the drift velocity $`\overline{v}_\pm `$ of a free vortex to be proportional to $`I`$.

In the RSJ dynamics, the spin-waves have a wavelength-independent relaxation time, of the order of the basic time unit $`\sigma \mathrm{}/2eI_c`$. Hence the equilibrium interaction between vortices in (6) is un-retarded in the dynamic case. Equation (8), on the other hand, is based on the assumption that there is no spatial correlation between vortices and antivortices. It turns out that this assumption is violated when the width of the strip is less than a certain characteristic size $`l_\mathrm{r}=\rho _{\mathrm{}}^{1/2}(I_c/I)^{\chi /2}`$. In a bulk system, there is typically one free vortex/antivortex pair in an area of linear size $`l_\mathrm{r}`$ at any given time. This is no longer the case when $`L_y<l_\mathrm{r}`$. Under periodic boundary conditions in the $`y`$-direction, a vortex may annhilate with an antivortex created from the same bound pair. Consequently, the life time of the two are shortened, giving rise to a spurious lower voltage as seen in Fig. 2.

To obtain a semi-quantitative estimate of the voltage reduction, let us consider a simplified model of vortex-antivortex recombination. In a coarse-grained description, we take $`l_\mathrm{b}`$ to be the basic unit of length. The time for a single vortex to traverse the width of the strip is $`\tau =L_y/\overline{v}_\pm =L_yl_\mathrm{b}`$. During this period, the typical displacement of the vortex in the $`x`$-direction is $`l=\tau ^{1/2}`$. The probability for the vortex to be within a distance $`l_\mathrm{b}`$ from its initial position after one round is thus $`p_1=l_\mathrm{b}/l=(l_\mathrm{b}/L_y)^{1/2}`$. This is the probability for a vortex-antivortex pair to recombine in one round across the strip. Assuming $`p_11`$, the probability for this to happen in two rounds is $`p_2(l_\mathrm{b}/2L_y)^{1/2}`$. In general, $`p_np_1n^{1/2}`$ provided $`p_1n^{1/2}<1`$. The survival probability after $`n`$ rounds is $`1_{i=1}^np_i1p_1n^{1/2}`$. This yields a self-recombination time,

$$\tau _sp_1^2\tau =L_y^2.$$
(10)
The density of free vortices in this case is $`\rho _s\mathrm{\Gamma }\tau _sL_y^2(I/I_c)^{2\chi }`$. From Eq. (5) we obtain,

$$V_\mathrm{s}\rho _s\overline{v}_\pm L_y^2(I/I_c)^{2\chi +1},$$
(11)
which increases with increasing $`L_y`$.

The typical displacement of the vortex in the $`x`$-direction over $`\tau _s`$ is $`\tau _s^{1/2}=L_y`$. To preempt self-recombination, we need $`\mathrm{\Gamma }L_y^2\tau _s1`$, i.e., creation of the second vortex/antivortex pair in the same region in space within the time $`\tau _s`$. This yields $`L_y=(I_c/I)^{\chi /2}=l_\mathrm{r}`$! Therefore Eq. (11) applies in the intermediate regime $`l_\mathrm{b}L_yl_\mathrm{r}`$.

The situation is somewhat different when $`L_y<l_\mathrm{b}`$. In this case, the pair-breaking interaction in (6) due to the external current is always weaker than the vortex-antivortex attraction and hence can be treated as a perturbation. The $`I`$-$`V`$ curves can in principle be calculated in a linear response theory. Qualitatively, though, we may estimate the finite-size effect from the equilibrium relaxation time of a vortex-antivortex pair on scale $`L_y`$,

$$\tau _{\mathrm{eq}}L_y^z,$$
(12)
where $`z=2\chi 2`$ is the dynamic exponent. We now divide the strip into square segments of $`L_y`$ columns each. The rate of creating a vortex-antivortex pair of size $`L_y`$ in such a segment is $`\tau _{\mathrm{eq}}^1`$. In thermal equilibrium, the vortex/antivortex may drift in either direction around the strip and recombine to generate a $`\pm 2\pi `$ phase jump. This symmetry is broken by the external current which enhances the $`2\pi `$ jump (along the current direction) by a factor equal to the ratio between the energy difference of the two possibilities and $`k_BT`$. This leads to the estimate,

$$V_{\mathrm{Ohm}}\frac{1}{\tau _{\mathrm{eq}}L_y}\left(\frac{IJL_y}{I_ck_BT}\right)(I/I_c)L_y^z.$$
(13)
In the Ohmic regime, $`V`$ decreases rapidly with increasing $`L_y`$, as seen in Fig. 2.

At $`L_y=l_\mathrm{b}`$, the two expressions (11) and (13) coincide to give the minimum value $`V_m=(I/I_c)^{1+z}`$. This is indeed the power-law form obtained by Minnhagen et al. by assuming $`l_\mathrm{b}`$ as the only relevant length scale.

A direct consequence of the self-recombination picture proposed above is the intermittent phase difference jumps across a square array of size less than $`l_\mathrm{r}`$. In the steady state, the phase difference $`\mathrm{\Delta }\varphi `$ across the system, averaged in the direction perpendicular to the current, is expected to grow linearly with time. However, if there is only at most one pair of free vortex/antivortex present in the system at any given time, the growth of $`\mathrm{\Delta }\varphi `$ is intermittent. This is indeed the case for $`L_y<l_\mathrm{r}`$, where the life time of a free vortex/antivortex pair is less than the time it takes to create the pair. The step-wise growth of $`\mathrm{\Delta }\varphi `$ has been observed in a previous numerical study using a square array which we have confirmed. A similar phenomenon occurs in the $`L_y<l_\mathrm{b}`$ regime, although in this case there are occasionally reverse steps. For a long strip, however, $`\mathrm{\Delta }\varphi `$ is a much smoother function of time due to averaging along the strip.

In summary, extensive simulations of the Josephson junction array in strip geometry reveal the existence of three size regimes with distinct current-voltage relationships. Arrays of size less than $`l_\mathrm{b}I^1`$ are in the Ohmic regime and obey dynamic scaling as proposed by Minnhagen et al. Between $`l_\mathrm{b}`$ and a second scale $`l_\mathrm{r}I^{\chi /2}`$, a reversed size-effect is observed. The power-law $`I`$-$`V`$ dependence of AHNS is shown to hold in the asymptotic regime. Our study reconciles a long-standing dispute between two theoretical approaches, and provides a framework for systematic analysis of finite-size effects in 2D superfluid and superconducting systems with equal number of vortices and antivortices.

This work is supported by the Research Grant Council (RGC) of the Hong Kong SAR under grant HKBU2061, by the Hong Kong Baptist University under grant FRG/98-99/II-74, and by the National Science Foundation under Grant No. PHY99-07949. One of us (QHC) is supported in part by the National Nature Science Foundation of China under Grant No. 10075039.