text:
# Confirmation of non-thermal hard X–ray excess in the Coma cluster from two epoch observations

## 1  Introduction

In the hierarchical scenario of structure formation, clusters of galaxies form by the gravitational merger of sub-clusters and groups. Numerical simulations of large scale structure formation indicate that this class of objects undergo several merger processes as they form (West, Villumsen, & Dekel 1991; Katz & White 1993). Cluster mergers are highly energetic events and the associated large scale shocks and turbulence could provide the ingredients necessary to the formation of extended radio regions (radio halos or relics) detected so far in a limited number of clusters, namely a magnetic field amplification and particle re-acceleration (Tribble 1993; Roettiger, Burns, & Stone 1999; Roettiger, Stone, & Burns 1999; Cavaliere, Menci, & Tozzi 1999; Brunetti et al. 2001; Fujita, Takizawa, & Sarazin 2003). Clusters of galaxies with detected diffuse radio emission indeed show significant evidence of merger activity. The existence of Mpc-scale radio halos or relics combined with the relatively short radiative lifetimes of the electrons ($`10^8`$ yrs) suggests an in-situ electron re-acceleration induced by a very recent or current merger event (Markevitch & Vikhlinin 2001).

The presence of large radio regions could be related to the origin of the non-thermal hard X–ray (HXR) emission detected in the Coma cluster (Fusco-Femiano et al. 1999; Rephaeli, Gruber, & Blanco 1999; Rephaeli & Gruber 2002) and Abell 2256 (catalog ) (Fusco-Femiano et al. 2000; Rephaeli & Gruber 2003) by BeppoSAX and RXTE and, at lower confidence level, by BeppoSAX in A754 (catalog ) (Fusco-Femiano et al. 2003). The most likely interpretation of the non-thermal HXR radiation is inverse Compton (IC) emission by the same radio synchrotron electrons scattering the cosmic microwave background (CMB) photons. However, one cannot completely exclude the possibility that the non-thermal emission detected by BeppoSAX and RXTE may be due to the presence of obscured sources in the field of view of the detectors. The BeppoSAX/MECS images test this possibility only in the cluster central region of size $`30^{}`$ in radius while the BeppoSAX/PDS (Frontera et al. 1997), which is able to detected HXR radiation in the energy range 15–200 keV, has a larger field of view (FWHM $`1.3^{}`$). The probability to find obscured sources, like Circinus (catalog ) (Matt et al. 1999) very active at high energies, in the field of view of the PDS is estimated to be of the order of 10% (Kaastra 1999; Fusco-Femiano et al. 2002). Future deep observations by IBIS on-board INTEGRAL with a spatial resolution of $`12^{}`$ can definitely resolve this uncertainty. A different mechanism given by non-thermal bremsstrahlung from supra-thermal electrons formed through the current acceleration of the thermal gas (En$`ß`$lin, Lieu, & Biermann 1999; Dogiel 2000; Dolag & En$`ß`$lin 2000; Sarazin & Kempner 2000; Blasi 2000; Liang, Dogiel, & Birkinshaw 2002) requires an unrealistically high energy input in order to maintain the HXR emission for more than $`10^8`$ yrs (Petrosian 2001;2002). In conclusion, the origin of the detected non-thermal HXR excesses in some clusters of galaxies seems to be restricted to a diffuse IC emission or to the presence of point sources in the field of view of the PDS.

In this Letter, we present the results of a long BeppoSAX observation of $``$ 300 ksec that confirms the presence of a non-thermal HXR tail in the spectrum of the Coma cluster in agreement with the detection obtained by a previous shorter observation of $``$ 91 ksec, after data re-analysis of the first observation slightly modified the numerical results reported in Fusco-Femiano et al. 1999, but not the general conclusions. Finally, we show the combined spectrum obtained by summing the spectra of the two observations.

Throughout this Letter we assume a Hubble constant of $`H_o`$ = 50 km s-1 Mpc-1 $`h_{50}`$ and $`q_0=1/2`$, so that an angular distance of $`1^{}`$ corresponds to 40.6 kpc ($`z_{Coma}=0.0232`$). Quoted confidence intervals are at $`90\%`$ level, if not otherwise specified.

## 2 PDS Data Reduction and Results

The Coma cluster was observed for the first time in December 1997 for $``$91 ksec and re-observed in December 2000 for $``$300 ksec. The pointing coordinates of BeppoSAX are at J(2000): $`\alpha :12^h58^m52^s`$; $`\delta :+27^{}58^{}54^{\prime \prime }`$. The total effective exposure times of the PDS in the two observations were 44.5 ksec and 122.2 ksec, respectively (hereafter OBS1 and OBS2).

The PDS spectra of both the observations were extracted using the XAS v2.1 package (Chiappetti & Dal Fiume 1997). The choice of using this software is dictated by the non standard pipeline needed to extract the net count spectra (see below). Because our Institute (i.e., IASF/Bologna) was in charge of the design, construction and maintenance of the PDS, and we developed and tested the XAS package, specifically created to handle the PDS peculiarities (while the SAXDAS package, used for the standard analysis, is more suitable for handling imaging instruments, like MECS and LECS), we felt more confident in using XAS for the PDS analysis.

Since the source is rather faint in the PDS band ($``$5 mCrab in 15–100 keV) a careful check of the background subtraction must be performed. The background sampling was performed by making use of the default rocking law of the two PDS collimators that samples ON/+OFF, ON/–OFF fields for each collimator with a dwell time of 96 sec (Frontera et al. 1997a). When one collimator is pointing ON source, the other collimator is pointing toward one of the two OFF positions. Initially, we used the standard procedure to obtain PDS spectra (Dal Fiume et al. 1997); this procedure consists of extracting one accumulated spectrum for each unit for each collimator position. We then checked the two independently accumulated background spectra in the two different +/–OFF sky directions, offset by $`210^{}`$ with respect to the on-axis pointing direction (+OFF pointing: $`\alpha :12^h58^m57.8^s`$; $`\delta :+24^{}28^{}55^{\prime \prime }.1`$ –OFF pointing: $`\alpha :12^h58^m47.0^s`$; $`\delta :+31^{}28^{}54^{\prime \prime }.7`$).

The comparison between the two accumulated backgrounds (difference between the +OFF and –OFF count rate spectra) shows that for OBS1 the difference is compatible with zero ($`0.044\pm 0.047`$ cts/s for a background level of $`21.66\pm 0.02`$ cts/s in 15–100 keV), while for the longer, more sensitive OBS2, there is an excess of $`0.064\pm 0.021`$ cts/s (background $`16.76\pm 0.01`$ cts/s)111The $``$20% variation in the PDS background is due to the BeppoSAX orbital decay: the lower orbit for OBS2 increased the shielding to ambient particles, therefore lowering the diffuse background.. A careful check of possible variable sources in the PDS offset fields lead our attention to the BL Lac source 1ES 1255+244 (catalog ), present in the +OFF field. Luckily, this same source was observed by BeppoSAX on May 1998 in the framework of a spectral survey of BL Lacs by Beckmann et al. (2002), who, somewhat surprisingly, state ”for 1ES 1255+244 there are no PDS data”. Indeed, we retrieved the raw data from the ASI Scientific Data Center, extracted the PDS spectrum (the background has to be evaluated only on one offset field because the other is pointed exactly on Coma — the two sources are contaminating each other!) and found that the source is quite faint, consistent with zero flux being detected. Because of the very short exposure time ($``$ 3 ksec) it is only possible to give a 2$`\sigma `$ upper limit of 0.26 cts/s in 15–100 keV, corresponding to 1.6 mCrab, however compatible with the background excess measured in OBS2.

It is worth noticing that it is possible to derive a more stringent upper limit on the X–ray flux from 1ES 1255+44 at the time of our Coma OBS2 by assuming that all the contamination in the +OFF field comes from this source. It is straightforward to show that the $`0.064\pm 0.021`$ cts/s excess translates into $`428\pm 424`$ net counts when one takes into account the much shorter BL Lac observation (3340 sec) compared to our OBS2 and an upward factor two correction due to the $`40^{}`$ off-axis position of the source, corresponding to a 50% intensity reduction because of the triangular collimator response (Frontera et al. 1997b). At its face value this is consistent with the run of PDS counts as compared to the MECS counts in the BL Lac sample studied by Beckmann et al. (2002), barring time variability effects.

Returning to the Coma observations, to remain on the safe side we decided to exclude the +OFF field in the background evaluation, and consider only the –OFF field as the ”un-contaminated” background for both the Coma observations. Moreover, just in the center of the +OFF field is also present the extremely weak ROSAT source RX J125847.1+242741 (catalog ). However, in section 3 we also report the level of confidence of the non-thermal excess considering the average of the measured backgrounds in the two positions.

The observed count rate of OBS1 is 0.78$`\pm `$0.03 cts/s in the 15–100 keV energy range, at a confidence level of $`26\sigma `$. In the first analysis the non-thermal excess with respect to the thermal bremsstrahlung emission reported by the PDS was at the confidence level of $`4.5\sigma `$. The derived non-thermal flux was $`2.2\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ in the 20–80 keV energy band (assuming a photon index $`\mathrm{\Gamma }_X`$ = 1.5). By relating the radio and the non-thermal fluxes it was then possible to estimate a volume-averaged intracluster magnetic field $`B_X0.15\mu G`$, using only observables (Fusco-Femiano et al. 1999). It should be pointed out that a re-analysis of these data has evidenced a trivial mistake in the previous data analysis (summing three spectra, one of them was summed twice), so that the correct spectrum obtained from OBS1 shows an excess with respect to the thermal component with the average gas temperature measured by Ginga (8.11$`\pm `$0.07, 90%; David et al. 1993) at a somewhat lower confidence level of $`3.4\sigma `$ (see Table 1). The fit with a single temperature gives $`9.9_{1.1}^{+1.3}`$ keV, above the average gas temperature measured by Ginga (with a field of view comparable to that of the PDS), implying the presence of a second spectral component. The fit with two thermal components (one fixed at 8.1 keV) requires an unrealistic second temperature ($`>`$ 50 keV) that strongly supports a non-thermal mechanism for the additional component present in the spectrum of the Coma cluster. If we consider a power-law for the second component, the PDS data are not able to fix the photon index, but the non-thermal flux is rather stable against index variations. We assume a photon index $`\mathrm{\Gamma }_X`$ = 2.0 to derive the non-thermal flux that results to be $`(2.3\pm 1.0)\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ in the 20–80 keV energy range.

The observed count rate of OBS2 is 0.72$`\pm `$0.02 cts/s in the 15–100 keV energy range, at the confidence level of $`36\sigma `$. At energies above 20 keV the spectrum shows an excess with respect to the thermal emission (kT = 8.1 keV) at a confidence level of $`3.4\sigma `$ (see Table 1). The fit with a single temperature gives $`9.5_{0.6}^{+0.8}`$ keV. Also OBS2 indicates the presence of an additional spectral feature and also in this case the fit with a second thermal component requires unrealistic values for the temperature. The non-thermal flux ($`\mathrm{\Gamma }_X`$=2.0) is $`(1.3_{0.6}^{+0.5})\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ in the band 20–80 keV, consistent with the flux reported in OBS1. The non-thermal fluxes are (marginally) consistent also at a 68% confidence level: $`(2.3\pm 0.7)\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ in the first observation and $`(1.3_{0.4}^{+0.3})\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ in the second one.

The combined spectrum is obtained by summing the spectra of the two observations (see Fig. 1). The total count rate is $`0.740\pm 0.017`$ cts/s in the 15–100 keV energy range, at the confidence level of $`44\sigma `$. At energies $`>`$20 keV the HXR excess is at the confidence level of $`4.8\sigma `$ (see Table 1). Even the inclusion of a 1% systematic to the data, necessary for sources with high S/N ratio but not for faint sources like Coma (Frontera et al. 1997b), does not change the significance of our non-thermal HXR excess. The fit with a single thermal component gives 9.7$`\pm 0.6`$ keV, well above the average gas temperature measured by Ginga, with a statistically unacceptable $`\chi ^2`$ value (=2.1 for 8 d.o.f.). The presence of a second component is more evident from the $`\chi ^2`$ value that has a significant decrement when a second component, a power law, is added to the thermal component with kT=8.1 keV. The improvement passing from the first model ($`\chi _\nu ^2`$ = 4.10 for 9 d.o.f.) to the second one ($`\chi _\nu ^2`$ = 1.2 for 7 d.o.f.) is significant at more than 99.4% confidence level, according to the F-test. Also, the combined spectrum cannot be fitted with a second thermal component unless an unrealistic value for the temperature is assumed, thus supporting the non-thermal origin for this additional spectral feature. The non-thermal flux for $`\mathrm{\Gamma }_X`$=2.0 is $`(1.5\pm 0.5)\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ in the 20–80 keV energy range and it is rather stable assuming reasonable different values of the photon index. In fact, for $`\mathrm{\Gamma }`$ =1.5 the flux is $``$ 6% lower and for $`\mathrm{\Gamma }`$ = 2.5 is $``$ 15% higher.

## 3 Discussion

Non-thermal HXR emission has been reported in two BeppoSAX observations of the Coma cluster performed with a time interval of about three years. Both the observations indicate the presence of a non-thermal excess with respect to the thermal emission at a confidence level of $`3.4\sigma `$. The combined spectrum obtained by adding up the two spectra gives an excess at the confidence level of $`4.8\sigma `$. The spectra of the two observations have been obtained using only the uncontaminated background-accumulated pointing at the \[–OFF\] field. However, by considering the average of the background measurements in the two sky directions in both the observations, the excess is still significant at the level of $`3.9\sigma `$ (observed count rate = 0.324$`\pm `$0.013 cts/s, model predicted rate = 0.273 cts/s). The non-thermal fluxes measured in the two observations are consistent at the 90% confidence level and marginally at the 68% c.l..

The non-thermal flux derived from the combined spectrum, $`(1.5\pm 0.5)\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ in the 20–80 keV energy range, is consistent with our previously published detection of $`(2.2\pm 0.8)\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ and with the value of $`(1.2\pm 0.3)\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ measured by RXTE in the same energy band (Rephaeli, Gruber, & Blanco 1999) and confirmed by a second deeper observation (Rephaeli & Gruber 2002). The non-thermal flux value and the relative c.l. of the excess have been obtained using the Ginga measurement of 8.11$`\pm `$0.07 keV (David et al. 1993) for the average gas temperature that is in good agreement with the XMM-Newton determination (Arnaud et al. 2001) of 8.25$`\pm `$0.10 keV in the central region (R$`<10^{}`$) of the cluster. In fact, all the X-ray observations of Coma (Hughes, Fabricant, & Gorenstein 1998a; Hughes et al. 1988b; Watt et al. 1992; Hughes et al. 1993) have indicated a clear pattern: the larger the field of view, the lower the measured temperature. Besides, RXTE reports a best-fit temperature of 7.90$`\pm `$0.03 keV (Rephaeli & Gruber 2002) in a field of view of $`1^{}`$ comparable to that of the PDS. However, also considering kT=8.25 keV for the average gas temperature in the field of view of the PDS, the non-thermal excess is at the level of $`4.6\sigma `$ (observed count rate = 0.349$`\pm `$0.015 cts/s, model predicted rate = 0.280 cts/s) and the derived non-thermal flux has a negligible variation.

Recent PDS data analysis of Coma performed with the SAXDAS software have lead to controversial results: an analysis of both observations has not reported evidence for a non-thermal excess (Rossetti & Molendi 2003), while the analysis of the first one (OBS1) by Nevalainen et al. (2003) confirms our published HXR detection, albeit at lower confidence level for the systematic uncertainties of their work. A systematic comparison between PDS spectra extracted by means of the two software packages (XAS, used here, and SAXDAS) is under way (Landi et al. , 2003). Preliminary results on the analysis performed on sources of different luminosities show that the spectral parameters do not change when computed with different packages. On the other hand, the errors associated to the spectral parameters are smaller when using XAS. This effect is more accentuated for faint sources. We suspect that this could be due to differences in filtering of good data and/or in the spectral equalization (i.e., conversion from spectral channels to energy channels) that for SAXDAS is performed after summing the four PDS units while for XAS is performed before. These effects will be discussed into details in a forthcoming paper (Landi et al. in preparation).

As discussed in the Introduction, the likely origin of the non-thermal HXR excesses detected by BeppoSAX and RXTE is IC emission by the same relativistic electrons responsible for the diffuse radio emission scattering the CMB photons, as predicted in the 1970‘s (see Perola & Reinhardt 1972; Rephaeli 1979). The probability to find obscured sources in the field of view of the PDS may be lower than $``$ 10%, as estimated by Kaastra et al. (1999) and Fusco-Femiano et al. (2002), considering that non-thermal HXR emission has been detected in at least two clusters, Coma and A2256, both showing extended radio emission. A detailed search by Nevalainen et al. (2003) excludes a significant contamination from obscured AGN present in the FOV of Coma and A2256 and supports an indication for an extended non thermal emission. If we take into account the estimated contribution by AGN in Coma: (9$`\pm `$6)$`\times `$ 10-3 cts/s in the 20–80 keV energy range, the confidence level of the non-thermal excess is in the interval $`(3.94.7)\sigma `$.

In the framework of the IC model the combination of the radio and non-thermal X–ray fluxes allows an estimate for a volume-averaged intracluster magnetic field $`B_X`$ of $`0.2\mu G`$ (see Fusco-Femiano et al. 1999). This value seems to be in contrast with the line of sight magnetic field derived from the Faraday rotation of polarized radiation of sources through the ICM ($`B_{FR}6\mu G`$; Feretti et al. 1995). Newman, Newman & Rephaeli (2002) have recently pointed out that many and large uncertainties are associated with the determination of $`B_{FR}`$ (see also Govoni et al. 2002 and Govoni & Murgia 2003). However, this discrepancy can be attenuated by considering models that include the effects of more realistic electron spectra, spatial profiles of the magnetic fields and anisotropies in the pitch angle distribution of the electrons (Goldshmidt & Rephaeli 1993; Brunetti et al. 2001; Petrosian 2001; Kuo, Hwang, & Ip 2003). The present value of the non-thermal HXR flux is slightly lower than that reported in Fusco-Femiano et al. (1999) favoring a central magnetic field strength of 1–2 $`\mu `$G in the two-phase model of Brunetti et al. (2001), more consistent with the $`B_{FR}`$ values.

The alternative between primary and secondary electrons as responsible for non-thermal phenomena in clusters of galaxies is discussed in many papers. Primary electrons may be injected in the ICM of the Coma cluster by some processes (starbursts, AGNs, shocks, turbulence) during a first phase and re-accelerated during a second phase (Brunetti et al. 2001). Secondary electrons may be due to decay of charged pions generated in cosmic ray collisions within the ICM (Dennison 1980; Blasi & Colafrancesco 1999; Dolag & En$`ß`$lin 2000; Miniati et al. 2001; Miniati 2003). Radio and HXR spectral properties of Coma provide observational constraints able to discriminate between these two different populations of electrons (Brunetti 2002). In particular, the derived volume-averaged intracluster magnetic field of $`0.2\mu G`$ implies relativistic electrons at energies $`\gamma 10^4`$ to explain the observed diffuse synchrotron emission. At these energies IC losses may determine a cutoff in the spectrum of the accelerated electrons as supported by the radio spectral cutoff observed in Coma (Deiss et al. 1997). The cutoff in the electron spectrum may be naturally accounted for in the context of re-acceleration models, while it is not expected if the radio emission is due to a continuous production of secondary electrons. More recently, a radio spectral cutoff has been found also in the case of A754 by relating the VLA observation at 1.4 GHz (Bacchi et al. 2003) to the observations of Kassim et al. (2001) at lower frequencies. This cluster also shows non-thermal HXR radiation detected at a confidence level slightly above 3$`\sigma `$ by BeppoSAX (Fusco-Femiano et al. 2003) and the derived value of the magnetic field is of the same order of that determined in Coma. The PDS detection should be confirmed by a deeper observation with imaging instruments for the presence of the radio galaxy 26W20 located at a distance of $`27^{}`$ from the BeppoSAX pointing. IBIS on-board INTEGRAL with its spatial resolution of $`12^{}`$ has the possibility to eliminate this ambiguity and to detect the excess at a higher confidence level with respect to that obtained by BeppoSAX.

## 4 Acknowledgments

We wish to thank F. Frontera for stimulating discussions and the referee for the useful suggestions. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France, and of data retrieved from the ASI Scientific Data Center operated at the ESA establishment of ESRIN, Frascati, Italy.


format: text

text:
# How to Claim a Discovery

## I Introduction

Setting limits for new particles or decay modes has been an active research area for many years. In high energy physics it received renewed interest with the unified method by Feldman and Cousins Cousins-Feldman (1). Giunti Giunti (2) and Roe and Woodroofe Roe-Woodroofe (3) gave variations of the unified method, trying to resolve an apparent anomaly when there are fewer events in the signal region than expected. They all discuss the problem of setting limits for the case of a known background rate. The case of an unknown background rate was discussed in a conference talk by Feldman Feldman (4) and a method for handling this case was developed by Rolke and López Rolke and Lopez (5). Little work has been done though on the question of claiming a discovery. This problem could be handled by finding a confidence interval and claiming a discovery if the lower limit is positive. Instead the question of a discovery should be done separately, by performing a hypothesis test with the null hypothesis $`H_o`$:”There is no signal present”. Rejecting this hypothesis will then lead to a claim for a new discovery. In carrying out a hypothesis test one needs to decide on the type I error probability $`\alpha `$, the probability of falsely rejecting the null hypothesis. This is of course equivalent to the major mistake to be guarded against, namely that of falsely claiming a discovery.

In practice a hypothesis test is often carried out by finding the p-value. This is the probability that an identical experiment will yield a result as extreme (with respect to the null hypothesis) or even more so given that the null hypothesis is true. Then if $`p<\alpha `$ we reject $`H_0`$; otherwise we fail to do so. For the test discussed here it is not possible to compute the p-value analytically, and therefore we will find the p-value via Monte Carlo.

Maybe the most important decision in carrying out a hypothesis test is the choice of $`\alpha `$, or what we might call the discovery threshold. As we shall see, this decision is made much easier by the method described here because we will need only one threshold, regardless of how the analysis was done. What a proper discovery threshold should be in high energy physics is a question outside the scope of this paper, although we might suggest $`\alpha =0.001`$ (roughly equivalent to $`3\sigma `$). Sinervo Sinervo (6) argues for a much stricter standard of $`5\sigma `$, or $`\alpha =2.910^7`$. We believe that such extreme values were used in the past because it was felt that the calculated p values were biased downward by the analysis process, and a small $`\alpha `$ was needed in order to compensate for any unwittingly introduced biases. If we were to trust that our p-value is in fact correct, a $`1`$ in $`1000`$ error rate should to be acceptable.

A general introduction to hypothesis testing with applications to high energy physics is given in Sinervo Sinervo (6). A classic reference for the theory of hypothesis testing is Lehmann Lehmann (7).

## II The Signal Test

Our test uses $`T=xb`$ or $`T=xy/\tau `$ as the test statistic, depending on whether the background rate $`b`$ is assumed to be known or not. Here $`x`$ is the number of observations in the signal region, $`y`$ is the number of observations in the background region and $`\tau `$ is the probability that a background event falls into the background region divided by the probability that it falls into the signal region. Therefore $`y/\tau `$ is the estimated background in the signal region and $`xy/\tau `$ is an estimate for the signal rate $`\lambda `$. $`T`$ is the maximum likelihood estimator of $`\lambda `$, and it is the quantity used in Feldman and Cousins Cousins-Feldman (1) without being set to $`0`$ when $`xy/\tau `$ is negative. This is not necessary here because a negative value of $`xy/\tau `$ will clearly lead to a failure to reject $`H_0`$.

Other choices for the test statistic are of course possible. For example, a measure for the size of a signal that is often used in high energy physics is $`S/\sqrt{b}`$. Under the null hypothesis this statistic is approximately Gaussian, at least if there is sufficient data. Unfortunately the approximation is not sufficiently good in the extreme tails where a new discovery is made, leading to p-values that are much smaller than is warranted. Even when using Monte Carlo to compute the true p-value, this test statistic can be shown to be inferior to the one proposed in our method because it has consistently lower power, that is its probability of detecting a real signal is smaller.

In order to find the p-value of the test we need to know the null distribution. In the simplest case of a known background rate and everything else fixed this is given by the Poisson distribution, but in all other cases it is not possible to compute the null distribution analytically, and we will therefore find it via Monte Carlo. As an illustration consider the following case shown in figure 1: here we have $`100`$ events on the interval $`[0,1]`$, with the signal region a priori set to be $`[0.44,0.56]`$. There are $`25`$ events in the signal region, and the background distribution is known to be flat. Therefore we find $`x=25`$, $`y=75`$, $`\tau =7.33`$ and $`T=14.77`$. Because we know that the background is flat on $`[0,1]`$, and because under the null hypothesis all $`100`$ events are background we can simulate this experiment by drawing $`100`$ observations from a uniform distribution on $`[0,1]`$ and computing $`T`$ for this Monte Carlo data set. Repeating this $`150000`$ times we find the histogram of Monte Carlo $`T`$ values shown in figure 2, case 1. In this simulation $`8`$ of the $`150000`$ simulation runs had a value of $`T`$ greater than $`14.77`$, or $`p=0.000053`$. Using $`\alpha =0.0001`$ we would therefore reject the null hypothesis and claim a discovery. Note that in addition to rejecting the null hypothesis we can also turn the p-value into a significance by using the Gaussian distribution and claim that this signal is a $`3.87\sigma `$ effect.

How would things change if the signal region had not been fixed a priori but instead was found by searching through all signal regions centered at $`0.5`$ and we would have accepted any signal with a width between $`0.01`$ and $`0.2`$? That is if we had kept the signal location fixed but find the signal width that maximizes $`T`$, the estimate of the number of signal events? Again we can find the null distribution via Monte Carlo, repeating the exact analysis for each simulation run individually. The histogram of $`T`$ values for this case is shown in figure 2, case 2. Here we find a value of $`T`$ larger than $`14.77`$ in $`570`$ of the $`150000`$ runs for a p-value of $`0.0038`$ or $`2.67\sigma `$. At a discovery threshold of $`\alpha =0.001`$ we would therefore not find this signal significant anymore.

Even more, what if we also let the signal location vary, say anywhere in $`[0.2,0.8]`$? That is for any pair of values $`(L,H)`$ we define $`[L,H]`$ as the signal region and $`[0,L),(H,1]`$ as the background region, compute $`T_{L,H}`$ for this pair and then maximize over all possible values of $`L`$ and $`H`$. Note that because $`T_{L,H}`$ is monotonically increasing in $`\tau `$ as long as all the observations stay either in the signal or in the background region, we can find the maximum fairly quickly by letting $`L`$ and $`H`$ be the actual observations. The histogram of $`T_{L,H}`$ values for this case is shown in figure 2, case 3. We find a value of $`T`$ larger than $`14.77`$ in $`9750`$ of the $`150000`$ runs for a p-value of $`0.065`$ or $`1.51\sigma `$, clearly not significant.

It was necessary in the second and third cases above to limit the search region somewhat, to the interval $`[0.2,0.8]`$ and to signals at least $`0.01`$ and at most $`0.2`$ wide, because otherwise the largest value of $`T`$ is almost always found for a very wide signal region, even when a clear narrow signal is present. This restriction will not induce a bias as long as the decision on where to search are made a priori.

In the general situation where the background is not flat on $`[0,1]`$ we can make use of the probability integral transform. Of course this requires knowledge of the background distribution $`F`$, but if it is not known we can estimate it from the data, either using a parametric function fitted to the data or even using a nonparametric density estimator. Again all calculations are done under the null hypothesis so we do not need to worry about the signal or its distribution.

As long as we copy exactly for the Monte Carlo events what was done for the real data we will find the correct p-value. This includes using cuts used to improve the signal to noise ratio, but it then requires the ability to correctly Monte Carlo all the variables used for cutting, including their correlations.

## III Performance of the Method

As an illustration for the performance of the signal test consider the following experiment: we generate $`100`$ events from a linear background on $`[3,5]`$ and (if present) a Gaussian signal at $`3.9`$ with a width of $`0.05`$. Then we find the signal through a variety of situations, from the one extreme where everything is fixed a priori to the other where the largest signal of any width is found. The background density is found by fitting and the background rate is estimated. The power curves are shown in figure 3. No matter what combination of items were fixed a priori or were used to maximize the test statistic, and with it the signal to noise ratio, all cases achieved the desired type I error probability, $`\alpha =0.05`$. Not surprisingly the more items are fixed a priori, the better the power of the test.

## IV Conclusion

We have described a statistical hypothesis test for the presence of a signal. Our test is conceptually simple and very flexible, allowing the researcher a wide variety of choices during the analysis stage. It will yield the correct type I error probability as long as the Monte Carlo used to find the null distribution exactly mirrors the steps taken for the data. Monte Carlo studies have shown that this method has satisfactory power.

## V
Acknowledgments

This work was partially supported by the Division of High Energy Physics of the US Department of Energy.

## VI
Appendix


text:
# Surface effects on nanowire transport: numerical investigation using the Boltzmann equation

## 1 Introduction

The effect of a finite system size on the conductivity of a material is a subject of considerable physical interest, which has recently been lent added relevance and importance by rapid developments in nanowire synthesis and assembly , electrical characterization and transport measurement methods . The development of nanowires currently represents an important part of materials and applied physics research. Hierarchial self-assembly techniques  envisaged in semiconductor nanowires make them promising central elements of future integrated electronics.

Reports of basic functional two- and three-terminal semiconductor nanowire devices including junctions, bipolar transistors and field-effect transistors are now widely found in the literature. The nano-scale transport properties of several important semiconductors including Si, , GaAs , GaN and Ge ), and semi-metals  have been investigated in detail. (Bi has attracted great attention because of its unique combination of interesting properties and its potential for thermoelectric applications .) Bicrystalline nanowires , crossed nanowire structures , functional networks  including ultra-high-density lattices , heterostructures  and superlattice devices  are part of the rapidly growing body of novel nanowire configurations under development.

A central aspect of theoretical enquiry must be the extent to which the conductivity of a nanowire differs from that of the bulk material. At first glance, it seems reasonable to suppose that the conductivity is smaller in a nanowire because of the addition of surface scattering, assuming that the band structure does not change drastically. However, experimental results have yielded conflicting indications on this point, which is presently not well-understood. Particularly intriguing are reports of mobility values higher than their bulk counterparts observed in silicon nanowires ; while in other cases the mobility has been deduced to be orders of magnitude lower.

Such challenging theoretical questions, brought into immediate relevance by the extensive data on electrical transport in nanowires compiled in the last few years, make a thorough quantitative investigation valuable at this point. The widespread pursuit of experiments pertinent to the surface effect on conductivity motivates a generic numerical description of the finite-size effect allowing both freedom and simplicity in the incorporation of nanowire characteristics.

For large- and moderate- sized nanowires operated at room temperature, semi-classical kinetic effects are expected to be as important as quantum mechanical effects like the modification of band structure. Thus a semi-classical study is both necessary and desirable in the common regime where the diameter of the nanowires is much larger than the thermal wavelength. While a quantum mechanical approach is indispensable to investigate conductance in very narrow nanowires (radii of a few nm) , it is neither viable nor suitable for larger nanowires with radii of the order of 100 nm. For the latter the use of the semi-classical Boltzmann equation is the appropriate method to approach a thorough quantitative study of finite-size effects.

These effects have previously been addressed theoretically only with simplistic assumptions regarding the surface. The analytical results most widely quoted are due to Chambers , who used kinetic-theoretical ideas to calculate the modification of the effective mean free path due to the presence of a surface and then used this value to find the conductivity. Essentially the same results were obtained using the Boltzmann equation with suitable boundary conditions in .

We present here the results obtained from a general numerical scheme we have developed to solve the Boltzmann equation in a cylindrical geometry. Such a direct numerical solution offers the requisite freedom in incorporating nanowire characteristics and conditions such as the equilibrium electron density profile, or the presence of defects and impurities. To retain conceptual simplicity, we employ the relaxation-time approximation as a first step towards the systematic modelling and prediction of nanowire conductivity with a view to relating model parameters to experimental data, and possibly casting a light on surface characterization . In this context, comprehensive experimental studies of Bi nanowires are of direct relevance . The introduction of diameter-controlled synthesis of nanowires  provides yet another fruitful context for our study.

We emphasize that surface scattering is just one of many scattering mechanisms that contribute to the resistivity of a nanowire; other mechanisms, especially acoustic phonon scattering, can be produce more dissipation in many circumstances. In addition, for studies of specific nanowires with specific surface defects and impurities, the relaxation time approximation employed here would ideally be replaced by a more exact approach such as Monte Carlo or quantum mechanical simulation as mentioned above. Our intention, however, is to contribute generic intuition about the behavior of surface scattering in large wires, even when the specific surface scattering centers are unknown and even when other resistive effects may be primary. For this reason, we use a generic relaxation time approach and map out its predictions for varied choices of wire characteristics. Despite its approximate character, this approach confers useful insight and is substantially more refined than the often-invoked ”specularity coefficient” model .

Inter alia, the results yield insights into the general problem of transport with a nominal relaxation time that varies with spatial coordinates, which is non-trivial because of the fact that the diffusion of carriers connects different spatial regions, making their properties inter-dependent, and thereby introduces a connectivity to the physical situation. This may be of direct relevance to transport in layered media such as magnetic multi-layers .

The numerical problem resulting from a finite-difference representation of the Boltzmann equation with a simple grid is solved using the conjugate-gradient method, which offers significant computational efficiency. A comparison with analytical results in limiting cases confirms the reliability of the scheme.

The paper is organized as follows: In section 2 we describe our computational framework, and the form of the Boltzmann equation adapted to the problem at hand. Our surface model is described in section 3. We present our results in the simple case of zero–magnetic-field in section 4, showing how they conform to physical expectation, thus validating the numerical scheme used. The important limit of zero-surface width is considered in section 5. Section 6 is devoted to the longitudinal magneto-conductance arising due to the surface effect and is followed by a summary (section 7).

## 2 Boltzmann Equation and Problem Specification

In our approach, we consider a cylindrical conductor which has a finite surface width, so that there is no abrupt change at the boundary to be dealt with through boundary conditions. In particular, this includes treating the unperturbed (i.e. when the external fields are zero) distribution function $`f_0`$ as a function of both space and momenta. In order to use the Boltzmann equation with such a distribution function $`f_0`$ we introduce an effective internal electric field $`𝐄_𝐜(𝐫)`$ in addition to any external field, to account for this spatial variation of $`f_0`$. In a confined cylindrical system where $`f_0`$ decays from its interior magnitude to zero over a finite surface width, the force due to this internal field is exactly analogous to the constraining force that keeps a particle within an enclosure with an abrupt boundary, and is therefore physically expected. With this additional internal field, the Boltzmann equation  for the distribution $`f(𝐩,𝐫)`$ of non-interacting carriers in the relaxation time approximation with a spatially varying relaxation time is

$$f\frac{𝐩}{m}+_𝐏fq(𝐄+𝐄_𝐜+\frac{𝐩}{m}\times 𝐁)=\frac{ff_0}{\tau (𝐫)}$$
(1)
Here E and B are the external electric field and magnetic field, and $`q`$ and $`m`$ are the charge and mass of each carrier. It is convenient to introduce the deviation $`\varphi (𝐩,𝐫)=f(𝐩,𝐫)f_0(𝐩,𝐫)`$ due to the presence of the external fields; the equation for $`\varphi `$ is

$$\varphi =\tau (𝐫)[\varphi \frac{𝐩}{m}+q_P\varphi (𝐄+𝐄_𝐂+\frac{𝐩}{m}\times 𝐁)+q_Pf_0(𝐄+\frac{𝐩}{m}\times 𝐁)].$$
(2)
We consider a long cylinder with no azimuthal or axial inhomogeneity and confine attention to the case where the external fields are uniform and parallel to the axis of the wire, defining a natural axis $`\widehat{z}`$ of reference: $`𝐄=E\widehat{z};𝐁=B\widehat{z}`$ (see Fig. 1.)

Thus the physical problem at hand requires only one spatial degree of freedom $`r=\sqrt{x^2+y^2}`$, although we have to treat all the momentum components. It is convenient to work with the local momentum components $`p_r=𝐩\widehat{r}`$ and $`p_\theta =𝐩\widehat{\theta }`$ rather than the canonical momenta in a cylindrical system, because the canonical momentum conjugate to $`\theta `$ is an angular momentum. Also, to preserve form (2) of the Boltzmann equation, $`f`$, $`f_0`$ and $`\varphi `$ are defined as coordinate densities in the space of $`(p_r,p_\theta ,p_z,r)`$ so that $`fdp_rdp_\theta dp_zdr`$ is the number of carriers in a volume element of this space. The missing volume element factor $`2\pi r`$ is absorbed into $`f`$, changing the internal field $`E_c`$, as we shall see. We solve equation (2) by discretizing on a real space grid of about $`20^4`$ grid points. The resulting matrix equation is solved by the method of conjugate gradients.

The drift velocity, which is the population average of the component of the carrier velocity parallel to the external electric field, is

$$v_z=\frac{p_z}{m}=\frac{1}{m}\frac{d^3p𝑑rf(𝐩,r)p_z}{d^3p𝑑rf(𝐩,r)}=\frac{1}{m}\frac{d^3p𝑑r\varphi (𝐩,r)p_z}{d^3p𝑑rf_0(𝐩,r)}$$
(3)
The effective conductivity $`\sigma `$ is $`qv_zN/lA`$ where $`N`$ is the total number of carriers in the wire, $`l`$ is its length and $`A`$ its cross-sectional area.

An explicit expression for $`E_c`$ can be directly deduced from the spatial variation of $`f_0`$ for an arbitrary form of the latter, since by definition of the unperturbed distribution,

$$f_0\frac{𝐩}{m}+_𝐏f_0q𝐄_𝐜=0$$
(4)
For simplicity, we assume that the spatial dependence of the distribution function is completely separable from the momentum dependence so that $`f_0=F(𝐩)\xi (𝐫)`$ and the internal field $`𝐄_𝐂`$ points only in the direction of $`r`$ : $`𝐄_𝐂=E_C(r)\widehat{r}`$. The computational results described below pertain to the case where the momentum-distribution is Maxwellian:

$$F(𝐩)=\mathrm{exp}(\frac{p^2}{2mkT})$$
(5)
Inserting $`f_0=F(𝐩)\xi (r)`$ and (5) in (4) we see that the simplest consistent form of the internal field is

$$E_c(r)=\frac{kT}{q}\frac{1}{\xi (r)}\frac{d\xi }{dr}.$$
(6)
## 3 Surface Model

We use a simple, continuous model to include the effects of the surface on the conductivity. The fact that there are no carriers beyond the wire radius is accounted for by taking the volume density of carriers to decay spatially as a Gaussian with a width $`w`$ beyond a certain radius $`r_0`$. A Gaussian function is chosen because of its mathematical simplicity and because its qualitative form is physically reasonable – we anticipate qualitatively valid physical conclusions using this form. This corresponds to a physical field that is proportional to the difference $`rr_0`$ for $`r>r_0`$ and is directed towards the centre. An additional term proportional to $`1/r`$ in the effective internal field $`E_c`$ arises on using equation (5) due to the inclusion of the volume factor $`2\pi r`$ in the definition of $`f`$, as mentioned in section 2.

$$f_0=F(𝐩)\xi (r)=\left\{\genfrac{}{}{0pt}{}{r\mathrm{exp}(\frac{p^2}{2mkT})\text{ }r<r_0}{r\mathrm{exp}(\frac{(rr_0)^2}{2w^2})\mathrm{exp}(\frac{p^2}{2mkT})\text{ }rr_0}\right\}$$
(7)
$$E_c(r)=\frac{kT}{qr_0}h(r)=\left\{\genfrac{}{}{0pt}{}{\frac{kT}{qr}\text{ }r<r_0}{\frac{kT(rr_0)}{qw^2}+\frac{kT}{qr}\text{ }rr_0}\right\}$$
(8)
The relative roughness of the surface is described by lower relaxation times $`\tau (r)`$ past a radius $`R`$. The change is also modelled by a Gaussian fall-off, with a width $`w_r`$. For simplicity, the radii $`R`$ and $`r_0`$ are assumed to coincide, and also the widths $`w`$ and $`w_r`$. The radial profile of the relaxation time is shown in Fig. 2.

$$\tau (r)=\tau _0\chi (r)=\left\{\genfrac{}{}{0pt}{}{\tau _0\text{ }r<R}{\tau _0\mathrm{exp}(\frac{(rR)^2}{2w^2})\text{ }rR}\right\}$$
(9)
We take $`w`$ to be of the order of $`0.1\times R`$.

## 4 Conductivity

In the absence of a magnetic field, it is possible to make certain predictions analytically about the behaviour of the conductivity. We work in the limit where $`E`$ is small, so that the drift velocity $`v_z`$ is always small in comparison to the thermal mean velocity $`v_{th}`$. If $`\tau (r)`$ were a constant, $`\tau _0`$, the drift velocity would be given by the familiar result

$$v_z=\frac{qE\tau _0}{m}$$
(10)
When $`\tau `$ acquires a non-trivial radial profile $`\tau (r)`$, its overall scale determines to what extent thin cylindrical shells at different radii affect their neighbours because of the diffusion of carriers from layer to layer. If this scale is so small that the mean free path of particles with velocities of the order of $`v_{th}`$ is much smaller than the length scale $`\zeta `$ over which $`\tau `$ changes significantly (i.e. $`\tau (r)v_{th}\zeta =w`$), the spatial connectivity is negligible, and one may treat the different layers separately. In this case, one can define an effective relaxation time $`\tau _{geom}`$ which depends only on the geometrical distribution of the relaxation time weighted by the relative carrier concentration $`\xi (r)`$. The use of $`\xi (r)`$, which is the spatial factor in the unperturbed distribution function $`f_0`$, here is consistent with the fact that the relaxation-time–scale is small.

$$\tau _{geom}=\frac{𝑑r\xi (r)\tau (r)}{𝑑r\xi (r)}$$
(11)
For the model discussed in section 3, this can be expressed simply in terms of the ratio $`\omega =w/R`$, providing a useful test for the numerical scheme in the limit of small $`\tau _0`$.

$$\tau _{geom}=\tau _0\frac{1+\sqrt{\pi }\omega +\omega ^2}{1+\sqrt{2\pi }\omega +2\omega ^2}$$
(12)
It is therefore useful to consider the essential surface effect as the departure of the actual drift velocity for the nanowire as a whole from the value derived by inserting $`\tau _{geom}`$ in place of $`\tau _0`$ in (3). This, as we have just observed, will be perceptible only when $`\tau _0v_{th}`$ is comparable to the surface width and significant only when comparable to the wire radius.

Now, it is important to note that the effect of spatially varying $`\tau `$ in a conductor as a whole is asymmetric in its action, between smooth (high $`\tau `$) and rough (low $`\tau `$) regions. A rough region is much less affected by the gradient of $`\tau (r)`$ in its neighbourhood than a smooth region. This is precisely because spatial connectivity is enhanced where $`\tau (r)v_{th}`$ is large as explained previously. Thus, in a conductor which has both regions of high $`\tau `$ and low $`\tau `$, the latter separately exhibit geometric (unconnnected) behaviour, showing no effect of the presence of the former.

But the effect of spatial variation on regions with high $`\tau `$ is to decrease mobility there, since carriers there can move enough to sample a significantly rougher region. Therefore, the effect of spatial variation on rough regions being negligible, its effect on a conductor as a whole, too, is a decrease in mobility. This consideration will later appear prominently in explaining the magneto-conductive effect as well (section 6).

Fig. 3 shows the variation of the conductivity as a function of the relaxation-time–scale for different temperatures. The ratio of the surface width to the wire radius $`\omega `$ is fixed at 0.2. As $`\tau _00`$, the conductivity tends to the same fraction of the bulk value at all temperatures, which is seen to be almost exactly equal to the geometric factor 0.8815, the value obtained by putting $`\omega =0.2`$ in equation (12). At higher temperatures, the departure from this ratio is also higher. This is physically expected since a higher temperature makes available higher radial velocities to the carriers in their random motion between collisions, thus increasing the communication between different layers. Thus we see that even if the relaxation time were independent of the temperature, the conductivity would have a temperature-dependence in the low field limit because of the surface effect. Note that to isolate the surface effect, the complex temperature variation of relaxation time in a real medium is deliberately suppressed, though it can be included easily in the computational scheme.

Further, it is seen in the inset that the data can be collapsed on to a single curve by using the transformation $`\tau _0\tau _0\sqrt{T/T_0}`$, where $`T_0`$ is an arbitrary temperature. In other words, when $`B=0`$ the temperature is a reducible parameter in the limit of a low electric field, which can be accounted for exactly by rescaling $`\tau `$ in proportion to the corresponding mean thermal velocity. This can be seen analytically by direct use of the Boltzmann equation, through a transformation to dimensionless variables. Rewriting eqn. (2) with the functional forms in (7), (8) and (9) we have

$$\varphi +qE\tau (r)\frac{\varphi }{p_z}+\frac{kT}{R}h(r)\tau (r)\frac{\varphi }{p_r}+\tau (r)\frac{p_r}{m}\frac{\varphi }{r}=\frac{qE\tau (r)p_z}{mkT}e^{p^2/2mkT}\xi (r)$$
(13)
Now introducing the variables $`𝐮=𝐩/\sqrt{mkT}`$ and $`s=r/R`$, we get

$$\varphi +\frac{qE\tau _0}{\sqrt{mkT}}\chi (s)\frac{\varphi }{u_z}+h(s)\frac{\tau _0}{R}\sqrt{\frac{kT}{m}}\chi (s)\frac{\varphi }{u_r}+\chi (s)u_r\frac{\tau _0}{R}\sqrt{\frac{kT}{m}}\frac{\varphi }{s}=\frac{qE\tau _0\chi (s)\xi (s)u_z}{\sqrt{mkT}}e^{u^2/2}$$
(14)
When E is small, $`\varphi `$ is of the same order of smallness, and therefore the second term in the left hand side of (14) is negligible. This leaves us with a linear differential operator containing $`\tau _0`$ only in the combination $`\tau _0\sqrt{T}`$ acting on $`\varphi `$ in the left hand side. The appearance of $`\tau _0`$ and $`T`$ in the right hand side (inhomogeneous term) of course, merely alters the overall scale of $`\varphi `$. Hence the function $`\varphi (𝐮,r)`$ may be written as $`\frac{qE\tau _0}{\sqrt{mkT}}\times H(𝐮,s;\tau _0\sqrt{T})`$ where H is the solution of (14) with the constant $`\frac{qE\tau _0}{\sqrt{mkT}}`$ absent in the right hand side. The drift velocity is

$`v_z`$ $`=`$ $`\sqrt{{\displaystyle \frac{kT}{m}}}{\displaystyle \frac{d^3u𝑑su_z\varphi (𝐮,s)}{d^3u𝑑s\xi (s)e^{u^2/2}}}`$
$`=`$ $`{\displaystyle \frac{Eq\tau _0}{m}}{\displaystyle \frac{d^3u𝑑su_zH(𝐮,s;\tau _0\sqrt{T})}{d^3u𝑑s\xi (s)e^{u^2/2}}}`$
$`=`$ $`v_z_0{\displaystyle \frac{d^3u𝑑su_zH(𝐮,s;\tau _0\sqrt{T})}{d^3u𝑑s\xi (s)e^{u^2/2}}}`$
Thus $`v_z/v_z_0(=\sigma /\sigma _0)`$ depends on $`\tau _0`$ only through the product $`\tau _0\sqrt{T}`$.

Fig. 4 shows the variation of the conductivity as a fraction of the corresponding bulk value with wire radius, parametrised by $`\tau _0`$. As $`\tau _0`$ decreases, the results approach geometrically determinable values; further, as $`\omega `$ increases, not only does the geometrical factor move further below $`1`$ but also the actual conductivity departs more and more from the geometrical value. Expectably, if the width is fixed, increasing the wire radius diminishes the strength of the surface effect.

The conformity with physical expectation and analytically known limits in the results above indicates the reliability of the numerical scheme and its suitability for other calculations based on the Boltzmann equation.

## 5 Limit of Zero Surface Width

A simple parametrization of the surface in a zero-width surface model is provided by the ”specularity coefficient” $`ϵ`$, which is the non-zero probability of carriers incident on the surface suffering a scattering event there: they undergo diffuse rather than specular reflection . This is often encountered in the literature; for instance, it has also been used in the analogous context of thermal transport by phonons .

We now consider the physically important limit of $`w0`$ which should relate the parameters in our model to the parameter $`ϵ`$. Fig. 5 shows the numerical results obtained by varying the surface width parameter, along with the graph of the analytical geometric factor. Significantly, it is seen that even for large $`\tau _0`$, the departure from the bulk value seems to go down to a small number, possibly zero. This is in contrast to  where the magnitude of the surface effect is unbounded as a function of the relaxation-time scale or mean free path. It may be noted that the computation becomes increasingly more expensive, as the surface width is decreased; we believe that the lowest width used here is low enough to allow us to draw qualitative conclusions from the results.

We thus see an apparent contradiction between the results obtained with two different characterisations of the surface, which demands closer scrutiny. A direct contradiction results only if $`ϵ`$ in the abrupt surface model is considered a free parameter assignable arbitrary values between 0 and 1; one could interpret the results we have seen as placing an upper limit on $`ϵ`$ not far above zero. But this, on the other hand, would imply that such a parametrization is not very meaningful. Indeed, a value of $`ϵ=1`$ is assumed in .

The contrast between the two approaches becomes understandable when one recognizes the correlation between $`ϵ`$ and the fraction of carriers resident within the surface layer. In a finite width model, this fraction is always finite; it approaches zero only in the zero width limit, when $`ϵ0`$. Equivalently, a non-zero value of $`ϵ`$ is only compatible with the existence of a finite fraction of carriers on the surface. Thus if a finite width model is accepted as more realistic, a value of $`ϵ`$ close to 1 is physically possible only when a large part of the total number of carriers resides on the surface even in the absence of a field.

There is another way of pointing out the essential difference between an abrupt surface model and a finite width one. In the former, it is assumed that a scattering at the surface necessarily returns the carrier into the bulk of the surface; whereas this is not the case in the latter, once again bringing to mind the existence of a finite surface population. Finally, in the former, there is no scale of the relaxation-time at which the surface effect saturates, whereas, in the finite width model, saturation occurs where the mean free path exceeds a few times the length scale $`\zeta `$. Saturation is physically expected in the latter case – if in some case, typical carriers can sample regions with significantly different values of $`\tau `$ during a lifetime, any further increase of the lifetime should not make a big difference.

A future investigation of the surface effect resulting from a carrier distribution involving non-trivial features near the surface rather than the simple, monotonic decay studied here is expected to provide further insight. In particular, it may provide an explanation for the apparent experimental evidence  that $`ϵ`$ is close to 1. In any case, the considerations stated in this section challenge the conception of the abrupt surface model, which ought to be re-evaluated in the light of these numerical results.

## 6 Longitudinal Magneto-conductance

The effect of a longitudinal magnetic field ($`B`$) on the conductivity of a large nanowire is of particular interest because the surface provides a classical kinetic mechanism for a magneto-resistive effect. In the large field limit, one would expect the magnetic field to decouple different spatial regions by confining the carriers kinematically, constraining them to helical motion between collisions. Thus, when the cyclotron radius is smaller than $`\zeta `$, the carriers are effectively localised by the magnetic field, so that the conductivity is largely determined by the geometric factor. As we have seen in section 4, the effect of spatial connectivity in a wire with a rough surface is to lower the conductivity below the geometric factor. Since a large $`B`$ field ought to reduce connectivity, the magneto-conductance clearly ought to be positive. Further, this magneto-conductance is expected to saturate with $`B`$ when the cyclotron radius for carriers with typical momentum values (around the thermal momentum $`p_{th}`$) becomes much smaller than $`\zeta `$. On the other hand, the magnetic field has an appreciable effect only when the reciprocal of the cyclotron frequency for carriers with typical momenta is smaller than the relaxation time; the magnetic field is ineffective if the probability for a collision before the carrier velocity turns round once is high.

At intermediate values of $`B`$, it is not obvious whether there can be a case of negative magneto-conductance; the field could conceivably tend to produce a net movement of carriers in some regions down the gradient of $`\tau (r)`$. The results presented in Figs. 6 and 7 suggest that this does not occur; the magneto-conductance is seen to be positive within the accuracy of the calculation in all cases.

For relaxation times corresponding to the bulk mobility of a common semi-conducting material like Si, the effect of the magnetic field is surprisingly small. In particular, the saturation values of $`B`$ exceed realistic laboratory values, whereas the magnitude of the fractional magneto-conductance is still only a few per cent. Again, this is in contrast to the results in . Although it is to be noted that the latter pertain to the degenerate limit (metallic case), the crux of the difference lies in the surface model, as argued in section 5. In short, if the surface effect is small, the magneto-conductive effect must also be limited by the corresponding departure from the geometric factor.

However, there still remains a qualitative similarity between our results and the results in ; for instance, we see that the low $`B`$ results lend themselves to a good parabolic fit, which is consistent with . Further, the arguments we have presented are substantiated by the fact that  systematically predicts a greater surface effect than the experimental results used for comparison, a discrepancy noted there itself. Since the discrepancy is despite the use of best-fit values, the parametrization itself must be regarded as dubious. It is important to note, however, that the results in  and our results are both consistent with recent experimental findings  to the extent that this surface magneto-resistive effect can be experimentally delineated from intrinsic quantum-mechanical effects.

The intriguing fact that the magnetic field increases conduction in all cases examined prompts a second look at the general effect of a spatially varying relaxation time. For one thing, using the same numerical scheme, we find another curious fact – the magneto-conductance does not flip sign when the profile of $`\tau (r)`$ is inverted, i.e. with a thin region of high mobility surrounding a core of low mobility, even in the high $`B`$ limit. More generally, we could not find a profile of $`\tau (r)`$ such as to yield a negative magneto-conductance.

We attribute the absence of a negative magneto-conductance to the asymmetry between smooth (high $`\tau `$) and rough (low $`\tau `$) regions in the effect of a spatial variation of $`\tau `$, discussed in section 4. As explained there, the surface effect, which is the result of connectivity in the model at hand, tends to decrease the conductivity below the geometric value. Thus the introduction of a magnetic field, whose basic action is to undo the transport connection between different regions by confining carriers, in reversing the effect of connectivity, can only increase the conductivity.

One also sees, by the same line of reasoning, that the surface effect must be small if $`\tau (r)`$ is assumed to be continuous, for, the region that contributes most to it is that where both $`\tau `$ and its gradient are large; but this region is of the order of the surface width w, which is expected to be small compared to the wire radius.

## 7 Conclusion

A direct solution of the Boltzmann transport equation offers a powerful approach to transport calculations in large nanowires. We have presented a simple finite width model in a cylindrical geometry and shown that despite its simplicity, its scope is significantly greater than that of an abrupt surface model, from which it exhibits qualitative differences. This challenges the utility of an abrupt surface model, especially when parametrized by the specularity coefficient. Our numerical results show that the classical magneto-conductive effect in large nanowires of materials like Si is limited to a few per cent even with magnetic fields beyond the range of a practical laboratory setup.

The authors thank P. Lammert, T. Mayer, S. Mohney, C. Nisoli, and J. Redwing for useful discussions. We acknowledge support under the NSF NIRT program grant DMR-0103068.


text:
# Scaling determination of the nonlinear 𝐼-𝑉 characteristics for 2D superconducting networks

## I Introduction

The physics of the static Kosterlitz-Thouless (KT) transition ktb  for quasi-two-dimensional (2D) superconductors is well understood in terms of vortex pair unbinding, since long time. minnhagen\_rev  However, in spite of this and in spite of the large interest in quasi-2D superconductors over recent years spurred by the event of high-$`T_c`$ superconductors, blatter\_rev ; minnhagen\_ferm  the dynamics of 2D vortices in the vicinity of the KT transition is still not completely settled. In the present paper we address the specific question of the value of the nonlinear $`I`$-$`V`$ exponent $`a`$ in the low temperature superconducting phase, which has been brought up by Tang et al. in two recent papers. tang1 ; tang2

The first attempt to describe the dynamics close to the KT transition was made by Ambegaokar-Halperin-Nelson-Siggia (AHNS). ahns  This attempt was based on the phenomenological assumption that the vortices could be separated into two distinct categories, i.e., free individual vortices and bound pairs of vortices and antivortices. This reasoning led to the AHNS value $`a=a_{\mathrm{AHNS}}`$ for the nonlinear $`I`$-$`V`$ exponent in $`VI^a`$ in the low temperature phase. ahns  An alternative value was later suggested by Minnhagen et al. minnhagen\_scale  based on the observation that the low temperature phase is quasicritical and as a consequence the exponent $`a`$ should follow from critical scaling. This gave the exponent $`a=a_{\mathrm{scale}}=z+1`$ where $`z`$ is the dynamic critical exponent. However, these two alternative values are different $`a_{\mathrm{scale}}>a_{\mathrm{AHNS}}`$ for $`T<T_{\mathrm{KT}}`$. The question is then which one is correct. This is the question discussed in the present paper.

A series of attempts to settle this issue was based on numerical simulations for various models displaying a KT transition in 2D: The Coulomb gas with Langevin dynamics, holmlund  the lattice Coulomb gas with Monte Carlo (MC) dynamics, weber  the $`XY`$ model with resistively-shunted Josephson junction (RSJ) dynamics. minnhagen\_scale ; kmo  These simulations all suggested that $`a_{\mathrm{scale}}`$ was the correct exponent.

However, for the 2D $`XY`$ model with RSJ dynamics in Ref. simkin, , Simkin and Kosterlitz argued that their data were more consistent with $`a_{\mathrm{AHNS}}`$. This discrepancy between the result of this simulation and the result obtained in the others cited above reflects the general difficulty in the determination of $`a`$: In a log-log plot the power law form of the nonlinear $`I`$-$`V`$ characteristics should give a straight line with the slope $`a`$. However, this form is only valid for small enough $`I`$ and the smaller the $`I`$ the larger system size is needed in order to get a size independent results. Thus great care is required to ensure that the data used in the determination of the slope are really size independent, which in turn makes the simulations demanding.

The experimental determination of $`a`$ is faced with the same type of finite size problems, as recently discussed in e.g. Ref.strachan, . However, in order to experimentally distinguish between the two different predictions of $`a`$ one needs in addition to analyze the temperature dependence of $`a`$ as described in Ref.minnhagen\_ferm, . This makes the experimental route to settle the issue somewhat difficult.

On the other hand, the size dependence of the data can also be turned into an advantage when determining $`a`$ from simulations:kmo  The logic here is that by carefully choosing the boundary condition the data should obey a size scaling from which $`a`$ can be determined provided the scaling assumption leading to $`a_{\mathrm{scale}}`$ is indeed correct. This strategy was used in Ref. kmo,  and provided strong evidence in favor of the critical scaling and $`a_{\mathrm{scale}}`$, as will be discussed in more detail in the present paper.

In an attempt to resolve the issue of the two different results for the value of $`a`$, Bormann in Ref. bormann,  reinvestigated the AHNS reasoning and concluded that it contains both of the results: For small enough current $`I`$ the result $`a=a_{\mathrm{AHNS}}`$ should be correct but as $`I`$ is increased there should be a crossover to a distinct region where $`a=a_{\mathrm{scale}}`$. In the $`a=a_{\mathrm{scale}}`$ regime the scaling should hold, so according to this analysis one should see a breakdown of scaling into an AHNS regime for small enough $`I`$. As found in Ref. kmo,  and which will be further discussed here, the scaling assumption holds for all the data obtained and there is no sign of a crossover to an AHNS regime.

The most recent attempt to settle this issue is by Chen, Tang, and Tong in Refs. tang1, ; tang2, . They again try to estimate the nonlinear $`I`$-$`V`$ exponent $`a`$, ($`VI^a`$) below the KT transition for a 2D Josephson junction array from the slope of $`\mathrm{ln}V`$ versus $`\mathrm{ln}I`$ for small $`I`$. The basic claim made is that an anomalous finite-size effect for small $`I`$ (meaning that in a certain parameter range the voltage for a fixed small $`I`$ increases with size instead of decreases) gives an overestimation of $`a`$. As a consequence it was concluded that $`a`$ is in better agreement with the AHNS prediction ahns instead of the dynamic scaling prediction minnhagen\_scale  concluded in Ref. kmo, ; choi, . However, as shown in Ref. kmo, , such an anomalous finite-size effect is a feature consistent with and emerging from the dynamical scaling and does consequently not affect the reliability of a determination based on finite-size scaling.

Our strategy to settle the issue is based on the observation that the dynamical scaling alternative is very amenable to testing by computer simulations. kmo  This is because the dynamical scaling makes direct predictions of the data obtained for finite size systems and does thus not hinge on any estimate of the asymptotic slope in the limit of small $`I`$ and large size $`L`$. Using this approach, a result consistent with AHNS requires that the dynamical scaling does in fact fail to describe the data. Alternatively, if one wants to verify a crossover to AHNS from a scaling regime, as suggested by Bormann in Ref. bormann, , then one needs to demonstrate that the dynamical scaling breaks down as one passes over into the AHNS regime. As shown in the present paper and in Ref. kmo, , direct tests of the dynamical scaling through simulations of the 2D $`XY`$ model with RSJ dynamics and MC dynamics give excellent agreement with dynamical scaling without any sign of the breakdown required for the AHNS alternative to be valid.

The content of the paper is as follows: In Sec. II we briefly recapitulate the finite-size dynamical scaling and in Sec. III the results of the simulations and the data analysis are presented. In Sec. IV we discuss our results in the context of other attempts to settle the issue and Sec. V, finally, contains some concluding remarks.

## II The dynamical scaling approach

The method used to determine $`a`$ directly from finite-size scaling is described in Ref. kmo, . It is based on the usual scaling form by Fischer et al. in Ref. fischer,  adopted to two dimensions and to finite-size scaling at criticality. Since the low temperature phase is quasicritical the scaling form applies at and below the KT transition. kmo ; dorsey  In order to take maximum advantage of the size scaling we reduce the number of large scales to one by considering a quadratic system with side $`L`$. The size scaling form is then given by kmo

$$v=i^{z+1}\left[\frac{f(Li)}{Li}\right]^z$$
(1)
where $`v=V/L`$ is the voltage per length across the sample, $`i=I/L`$ is the current density, $`z`$ is dynamic critical exponent and $`f(x)`$ is a scaling function which goes to a positive constant for small $`x`$ and is proportional to $`x`$ for large $`x`$. This means that $`vi^{z+1}`$ for a given small enough $`i`$ in the limit of large $`L`$. Consequently, the scaling prediction for the nonlinear $`I`$-$`V`$ exponent $`a`$ ($`VI^a`$) is $`a=z+1`$.

The use of this scaling approach has several advantages when trying to determine the exponent $`a`$ in model simulations: First of all, it does not hinge on the accuracy to calculate the small $`i`$-limit of voltages where at the same time $`Li>>1`$, which is notoriously difficult. Secondly, the existence of a scaling function $`f(x)`$ can be determined from a data collapse using data from all system sizes simultaneously. Thirdly, $`z`$ can be related to and obtained from the equilibrium properties of the system kmo . This gives an independent consistency check on the scaling given by Eq.(1), since this $`z`$-value has to agree with the one obtained directly from the data collapse.

## III Results for the 2D $`XY`$ model

Our main results are for the 2D $`XY`$ model with RSJ dynamics and finite temperatures $`T`$. This model undergoes a KT transition at $`T_c0.89`$ ($`T`$ is measured in units of the Josephson coupling and current in terms of the critical current $`i_c`$ of a single Josephson junction). The method used is described in Ref. kmo, : We use a fluctuating twist boundary condition (FTBC) and a square lattice with size $`L`$. The data presented here are well converged data for sizes up to $`L=256`$. We here analyze data obtained in the low temperature quasicritical phase somewhat below $`T_c`$.

In Fig. 1 we present results for $`T=0.8`$ in the alternative scaling form $`v/i^{z+1}=F(Li)`$ where $`F(x)`$ is related to $`f(x)`$ in Eq. (1) by $`F(x)=\left[f(x)/x\right]^z`$. In this plot we use the value $`z=3.4`$ obtained from equilibrium simulations using the relation minnhagen\_scale  $`z=1/(\stackrel{~}{ϵ}T^{CG})2`$ with the dielectric constant $`\stackrel{~}{ϵ}`$ and the Coulomb gas temperature $`T^{CG}`$, as described in Ref. kmo,  and given in Table I of Ref. kmo, value  As seen from Fig. 1(a) a very good data collapse is obtained. At higher values of $`Li`$ there is a notable systematic deviation from scaling. In order to determine the cause of this deviation the same data are plotted in Fig. 1(b) as a function of $`i`$. As seen the data which in Fig. 1(a) deviate from the scaling curve \[filled symbols in Fig. 1(a) and (b)\] now instead collapse and furthermore the collapse for all sizes starts approximately at the same value of $`i0.3`$. This means that the data for values $`i>0.3`$ are independent of $`L`$. The Josephson junctions become completely resistive for $`i=i_c`$ with a resistance given by normal state $`R_N`$ of the junction and which also means that $`V=R_NI`$ for the complete array. This large current relation is given by the full line in Fig. 1(b). It means that the deviation from the scaling curve in Fig. 1(a) is just the trivial crossover to the normal state resistance which always occurs for large enough $`i`$ because the small scale $`1/i_c`$ becomes relevant and breaks the scaling when $`i`$ becomes of the same order as $`i_c`$. The data in Fig. 1(a) show no deviation from finite-size scaling apart from the trivial crossover to normal state resistance as $`i`$ approaches $`i_c`$.

Figure 2 shows the scaling function $`f(Li)`$ in Eq. (1) obtained from the same data. In addition, we have in this figure also included data from the 2D $`XY`$ model with Monte Carlo (MC) dynamics. As demonstrated in Ref. katya, , MC dynamics gives for the same $`T`$ the same $`I`$-$`V`$ characteristics up to a constant factor and the same method as in Ref. katya,  is used here. The MC data are also for $`T=0.8`$ but also include the larger size $`L=512`$. The data for $`i`$ in the crossover towards $`vi`$ have been excluded in Fig. 2. As seen both sets of data collapse on the scaling curve with $`z=3.4`$. The full line ($`Li`$) corresponds to $`vi^{z+1}`$. Note that the data for a given size $`L`$ which fall on this line have the scaling slope $`a_{\mathrm{scale}}`$ in a log-log plot of $`v`$ versus $`i`$. This fact also gives an idea of the difficulty with a direct determination for a single system size: The larger system the fewer points can in practice be obtained. So although the three data points for $`L=512`$ fall nicely on the scaling function a determination of $`a`$ based on these three points alone would not have the same reliability.

The inset in Fig. 2 shows how the voltage $`v`$ depends on the size $`L`$ for a given fixed current. As seen the voltage for this fixed current increases for larger $`L`$. This is just a reflection of the fact that the scaling curve falls below the full line in Fig. 2 for $`Li`$ between approximately 1 and 10. This means that the increase of the voltage with increasing $`L`$ is a property of the scaling regime where the behavior is controlled by a single relevant large scale. Consequently, it is not connected to the appearance of a second relevant scale, as suggested in Ref.6 in case of the seemingly very similar increase of the voltage with increasing $`L`$ found in this paper.

## IV Comparisons and discussion

Our main result from the 2D $`XY`$ model with RSJ and MC dynamics is that the data collapse on a scaling curve for a certain value of $`z`$. This value of $`z`$ is in agreement with the scaling prediction which connects the value of $`z`$ to equilibrium quantities. minnhagen\_scale ; kmo  The data regime at higher current which do not fall on the scaling curve is controlled by the Josephson junction critical current $`i_c`$ and belong to the crossover to the normal state linear relation $`V=R_NI`$. The data show no crossover to an AHNS behavior.

According to Bormann in Ref. bormann,  such a crossover should be expected for small enough $`i`$. Thus, according to this theory, AHNS should be valid for low enough $`i`$ in the limit of large $`L`$ and there should be a crossover to the scaling result for larger $`i`$. As $`i`$ approaches the critical Josephson junction current there should be a second crossover towards the normal state linear $`I`$-$`V`$ characteristics. As seen from Fig. 2 the data give evidence of a crossover from the scaling regime to the linear Ohmic regime but there is no sign of a crossover towards the AHNS regime for smaller $`i`$. In Fig. 2 such a crossover would mean that the data for larger $`L`$ and smaller $`i`$ would fall above the scaling curve. The crossover current $`i_{\mathrm{cross}}`$ is proportional to $`1/n^{1/(z/21)}`$ where $`n`$ is the vortex density. An attempt towards a more quantitative estimate of the crossover is given in Fig. 1 of Ref. bormann, : Our data for $`i=0.01`$ correspond to a point in this figure $`(x,y)=(z+2,i^1)=(6.4,10^2)`$ which should be far inside the AHNS regime. In fact all the data up to the crossover to the linear regime (at $`i0.3`$) should, according to this estimate, belong to the AHNS regime. However, our data show no deviation from scaling for any size or current, indicating that all the data are instead in the scaling regime and that there is no crossover to an AHNS regime at small currents. However, one can, of course, not entirely rule out that such a crossover might exist for even lower currents and larger sizes than could be reached in the simulations.

Tang et al. in Ref. tang1, ; tang2,  use a rectangular sample with sides $`L_x`$ and $`L_y`$ and the current in the $`x`$ direction. The idea is to use a small ratio $`L_y/L_x`$ in order to minimize the dependence on the scale $`L_x`$ and they in practice use $`0.004L_y/L_x0.25`$ corresponding to $`L_y=8`$ and $`L_y=512`$, respectively. The current injection method used introduces a nonuniform vortex density in the current direction. However this is compensated by skipping a boundary region with the length $`b`$ at both boundaries when measuring the voltage. In contrast the use of FTBC is designed to make the vortex density uniform. Nevertheless, it was shown in Ref. choi,  that the method adopted by Tang et al. tang1 ; tang2  for PBC in the transverse direction gives the same result as the FTBC method provided care is taken to avoid any influence from the two additional length scales $`L_x`$ and $`b`$. Ref. choi,  confirms the scaling prediction for sizes up to $`L_y=64`$. Tang et al. tang2  argues that in the case of a periodic boundary condition (PBC) in the direction transverse to the current one should get a crossover from a scaling regime at lower currents towards an AHNS regime at higher. Note that this is precisely the opposite to the Bormann prediction. bormann  According to Tang et al. such a crossover is supported by their Fig. 6 in Ref. tang2, . This figure for $`T=0.8`$ corresponds precisely to our Fig. 1(a) for the same transverse sizes. The scaling is clearly visible although the quality in their Fig. 6 is not as good as in our Fig. 1a.magnetude  However, whereas our method is specifically designed to test the scaling and only introduces a single scale $`L`$, the method by Tang et al. tang2  have three large scales and in particular the fact that the ratio $`L_y/L_x`$ varies from 0.004 to 0.25 instead of being constant, might well cause deviations from scaling. In addition, there is a crossover for higher currents towards the linear $`I`$-$`V`$ characteristics. However, whereas, as shown above, the crossover in our case can unambiguously be attributed to the trivial crossover to the linear regime starting at around $`i0.3`$ for all sizes \[compare Fig. 1(b)\] the situation in Fig. 6 of Tang et al. tang2  is slightly different. At small values of $`L_y`$ the crossover comes at roughly the same value $`i0.3`$ (up to $`L32`$) and then it decreases to roughly 0.06 at $`L_y=512`$. This decrease of the crossover current with $`L_y`$ indicates that there is another length scale in the problem. Tang et al. suggest that this new length scale is associated with the vortex physics and is given $`L_ri^{(z+2)/4}`$ and that the crossover is to an AHNS regime occurring for $`L_y>L_r`$. The problem with this interpretation is, in the light of Fig. 1(a) and (b), that there is no such size dependent crossover for a square lattice. Since the bulk properties of vortex physics cannot depend on the shape of the sample, the complete absence of the appearance of an additional vortex length in Fig. 1(a) and (b) strongly suggests that the explanation offered by Tang et al. tang2  cannot be the correct one. An obvious candidate is instead the additional ratio $`L_y/L_x`$ introduced by Tang et al. and which is not kept fixed when $`L_y`$ is varied, as required by a proper scaling analysis. The point to note is that it is not the absolute value of $`L_x`$ per se that matters, but the fact that the ratio $`L_y/L_x`$ has to be constant, and as pointed out above it varies from very small to 0.25 in Ref.tang2, . This could explain both why the quality of the size scaling is not as good as for the scaling with fixed $`L_y/L_x=1`$ shown in our Fig. 1amagnetude , as well as the size dependence of the crossover current.

Tang et al. tang2  also obtained the data for an open boundary conditions (OBC) in the transverse direction. The difference with the PBC result is the different finite-size dependence of the voltage. The open boundary causes stronger boundary effects and thus changes the details of the size convergence. The data are plotted as $`v/i^{a_{\mathrm{AHNS}}}`$ versus $`iL_y`$ in Fig. 8 in Ref. tang2, . In Fig. 3(a) our data for FTBC and square lattice are plotted in the same way. As seen in Fig. 3(a) the data over a large region apparently fall on a horizontal line. This line corresponds to $`vi^{a_{\mathrm{AHNS}}}`$ and consequently one might be tempted to interpret this as evidence of an AHNS exponent. The fallacy here is that the data which constitute this line are the size converged data in the crossover region towards $`i_c`$ as shown in Fig. 3(b): The two data points connected with a line in Fig. 3(b) forms a spurious AHNS line by repetition in Fig. 3(a). When going from the steeper slope in the scaling region through the crossover to the linear $`I`$-$`V`$ dependence there will always be some current region where the slope has the AHNS value. However, this continuous crossover to a linear $`I`$-$`V`$ characteristics has nothing to do with the AHNS vortex physics. The important thing to notice in Fig. 3(a) is that the scaling systematically fails for the larger sizes at lower $`i`$ where the data for each size deviate to lower values with decreasing $`i`$ for a fixed size. The same deviation is observed in Fig. 8 in Ref. tang2,  where the data on the horizontal line are a repetition of the two data points in Fig. 7 from Ref. tang2,  for the current sizes $`i=0.3`$ and 0.4 (same currents as in Fig.3a), respectively. Thus for the larger system sizes, which are relevant for the bulk properties of vortex physics, the change of boundary condition to OBC makes little difference. This is of course expected because the boundary condition should not matter at all for large enough samples. For smaller samples the results are different, as apparent when comparing our Fig. 3(a) with Fig. 8 in Ref. tang2, . This is also expected because here the fluctuation associated with the open boundary increases the voltage for a given $`i`$ relative to the PBC boundary. However, the vortex physics in the bulk remains the same and it is these properties which dominate for large samples.

To summarize: Neither the shape of the sample nor the details of the boundary condition can change the vortex physics in the bulk for large enough systems. We find that there is no compelling evidence that the data by Tang et al. tang2  are in contradiction with this statement. Proper scaling analysis with a fixed ratio $`L_y/L_x`$, as shown in the present paper, gives no evidence of an extra scale in addition to $`L`$ and $`i_c`$.

## V Concluding remarks

As discussed in the present paper the $`I`$-$`V`$ characteristics for the 2D $`XY`$ model with RSJ and MC dynamics obey size scaling for a quadratic sample (scaling in $`L_y`$ with fixed at $`L_y/L_x=1`$). The fact that the data obey size scaling means that $`a=z+1`$ where $`z`$ is the dynamic critical exponent. The exponent $`z`$ can be calculated from the equilibrium properties of the system which together with the relation $`a=z+1`$ gives the scaling prediction for $`a`$. As shown this scaling prediction is in agreement with the $`z`$ value determined from the data collapse. The only deviation from the scaling can, as shown here, be linked to the trivial high current crossover to the linear $`I`$-$`V`$ relation.

A prediction for $`a`$ which is different from the scaling prediction can only be correct if the scaling breaks down. As shown in the present paper, there is no such breakdown in the scaling except for the one at high currents associated with the critical current of a single Josephson junction. This does in principle not rule out the possibility that a breakdown could occur at sizes larger and currents smaller than reached in the simulations. However, whereas the scaling prediction has been verified over a large parameter range, the breakdown of scaling, which would have to occur if the AHNS prediction was correct in the ultimate small current limit, has not been supported by any simulation data so far.

In our opinion, the only convincing way to verify that AHNS, or some other alternative prediction is correct in some parameter regime, is to demonstrate that the scaling breaks down in the same regime.

## Acknowledgements

Support from the Swedish Research Council is gratefully acknowledged. B.J.K acknowledges the support from the Korea Science and Engineering Foundation through Grant No. R14-2002-062-01000-0, and the support from Ajou University.