text:
# The Race to Build Supermassive Black Holes

## 1 Introduction

Central supermassive black holes are a common feature to galaxies today, but which came first, the black hole or the galaxy? Conventional thinking would suggest that the the first generation of stars evolved into black holes, which have subsequently settled to the centers of their host galaxies, merged, and accreted gas. But this idea, in which central black holes form inside pre-existing galaxies, has recently earned some scrutiny. First, the discovery of increasingly high redshift quasars requires a surprisingly early formation of the black holes (see, for example, Fan et al. (2001) and Fan et al. (2003)). Second, a large quasar sample shows no evidence of black holes growing in mass with decreasing redshift (Vestergaard, 2002; Dietrich et al., 2002). So we are left to consider the possibility that either the central black holes formed before their host galaxies, or they grew to maturity very quickly within them. Either way, they have grown little since the quasar epoch.

The most distant known quasar lies at $`z=6.41`$, with a central black hole of mass $`M_{}=3\times 10^9\mathrm{M}_{\mathrm{}}`$ (Willott, McLure, & Jarvis, 2003). In the $`\mathrm{\Lambda }`$CDM cosmology observed by WMAP (Bennett et al., 2003), with $`\mathrm{\Omega }_\mathrm{\Lambda }=0.73`$, $`\mathrm{\Omega }_\mathrm{m}=0.27`$, and $`H_0=71\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, this redshift corresponds to a time when the universe was only $`880\mathrm{Myr}`$ old. For the present work, we will take this as the time to beat: $`3`$ billion solar masses in $`880`$ million years.

In the past year, two separate HST studies have cited evidence for intermediate mass black holes (IMBHs) in the centers of globular clusters: a $`4000\mathrm{M}_{\mathrm{}}`$ hole in M15 (van der Marel et al., 2002), and a $`\mathrm{20\hspace{0.17em}000}\mathrm{M}_{\mathrm{}}`$ hole in Andromeda’s G1 cluster (Gebhardt, Rich, & Ho, 2002). This is the lastest and strongest evidence for IMBHs, but there is additional evidence, and good theoretical motivation as well; see van der Marel (2003) for a comprehensive review. IMBHs are widely believed to be a necessary growth stage for SMBHs. In section 2 of this paper, we will review the major proposed routes to the formation of a SMBH, all of which include an IMBH phase, from which time the question is simply one of growth.

We start in Section 2 with a flowchart of avenues leading to the creation of a SMBH. In Section 3, we examine the timescales for each needed process. We conclude in Section 4 by discussing how realistic each avenue is in light of these timescales.

## 2 Flowchart

There are essentially four proposed families of models leading to the formation of IMBHs, and two or three ways to grow them. These approaches are depicted in figure 1 and discussed in turn below.

1. The black holes may be *primordial*, in which case they formed from primordial density variations before big bang nucleosynthesis.

Primordial black holes (PBHs) have been studied extensively, and the most direct mechanism for their creation is the collapse of gaussian random density fluctuations (Carr & Hawking, 1974). These holes come from horizon scale (or smaller) modes, and therefore their masses are determined by their time of formation. In the radiation dominated early universe,

$$M_{}10^5\left(\frac{t}{\mathrm{s}}\right)\mathrm{M}_{\mathrm{}}.$$
(1)
But in order to preserve the successful BBN prediction of light element abundances, there must be no significant rate of PBH formation once nucleosynthesis begins, and therefore the PBHs are capped at intermediate mass. In addition, Carr, Gilbert, & Lidsey (1994) have pointed out that, given a small scalar spectral index $``$ $`n_s=0.93\pm 0.03`$ was recently observed in the CMB (Bennett et al., 2003) $``$ PBHs from density inhomogeneities should only have formed in quantities too small to be of interest.

A more promising, and perhaps inevitable mechanism for forming PBHs also exists, in which the collapse is triggered by “color holes” at the quark-hadron phase transition (Crawford & Schramm, 1982). However, because this occurred at $`10^6\mathrm{s}`$, these PBHs would be smaller than $`1\mathrm{M}_{\mathrm{}}`$ by eq. 1, and would remain as collisionless dark matter today, rather than collecting into larger black holes. (Interestingly, Hawkins (1996) shows evidence for such PBHs in the microlensing of distant quasars, in numbers comparable to that needed to account for dark matter.)

2. Normal *population III stars* formed at $`100\mathrm{M}_{\mathrm{}}`$, evolved to black holes, and merged at the center of their small dark matter halos.

This is perhaps the least exotic way to create IMBHs, and at this point there is very strong support for the process in the form of numerical simulations of structure formation (Abel, Bryan, & Norman, 2000; Yoshida et al., 2003). These simulations include the relevant atomic and molecular processes in the first gas clouds, particularly cooling by rotation of molecular hydrogen, superimposed on CDM halo evolution. They find that $`10`$ stars of $`100\mathrm{M}_{\mathrm{}}`$ form in halos of $`10^6\mathrm{M}_{\mathrm{}}`$, engaging $`1\%`$ of the system’s baryonic matter.

Because the cooling of population III stars is hindered by the lack of metals, these first stars would be more massive than those allowed by fragmented star formation today. Heger & Woosley (2002) have shown that these massive stars will evolve into black holes containing a significant fraction of the star’s initial mass (this fraction varies with the stellar mass, and is of order $`50\%`$), unless the stars are in the mass range $`140\mathrm{M}_{\mathrm{}}<M<260\mathrm{M}_{\mathrm{}}`$, in which case they are completely disrupted when they go supernova. Given the small initial size of such a cosmologically young halo ($`10\mathrm{pc}`$), the holes sink to the center and presumably merge into a single IMBH there.

3. *Supermassive stars* may have been the first baryonic objects to form. If so, they would have evolved rapidly into IMBHs.

Supermassive stars (SMSs), as a class of objects, span from $`10^3`$ to $`10^8\mathrm{M}_{\mathrm{}}`$, although the first generation of them would reside at the lower end of that range (Shapiro & Teukolsky (1983) provides a comprehensive introduction on the topic). A $`10^3\mathrm{M}_{\mathrm{}}`$ SMS has a lifetime of $`10^6`$ years (SMS lifetimes range down to 10 years), at the end of which, it undergoes relativistic collapse to a black hole in a matter of seconds (New & Shapiro, 2001). This collapse sends $`90\%`$ of the star’s original mass into the black hole, with the remaining $`10\%`$ forming an accretion disk (Shapiro & Shibata, 2002; Shibata & Shapiro, 2002).

The radius of a $`10^3\mathrm{M}_{\mathrm{}}`$ SMS, $`R_{\mathrm{SMS}}0.1\mathrm{pc}`$, is around $`10\%`$ of the virial radius for the earliest halos containing gas cool enough to collapse. It is therefore reasonable to suspect that these stars should appear in early halos, especially since population III stars produce ultraviolet radation which splits hydrogen molecules, thus destroying the only coolant which could collapse smaller stars. The collapsing gas might produce enough radiation pressure to prevent fragmentation (Baumgarte & Shapiro, 1999). Another possibility is that the first stars formed dense clusters which merged into one or more SMSs. However, whether or not supermassive stars ended the cosmological dark ages remains a matter of speculation.

4. If *dense clusters* of stars emerged from early star forming regions, then the cluster stars may have merged with each other to create massive stars (say, several hundred solar masses), which then evolved to massive seed black holes. This idea has been mapped out recently by Ebisuzaki et al. (2001), and we will not repeat it here. However, for the sake of timescales, which are the focus of this paper, this scenario is nearly the same as the population III star scenario discussed above, because the merging occurs during the lifetime of the cluster’s massive stars.

Each of the paths in figure 1 sees a phase with an IMBH imbedded in a halo, at which point accretion and halo mergers take over to achieve heavier black holes.

The merger scenario merges central black holes hierarchically along with the host halos, growing both together. It is clear that merging IMBHs and SMBHs are real, interesting phenomena deserving study. For example, Milosavljević & Merritt (2001) have simulated the whole process, and shown that coalescing SMBHs remove the steep cusps that CDM simulations generate, thus matching observed galaxy rotation curves. But it remained to be shown here whether or not mergers are prevalent enough to sufficiently increase black hole mass.

Gas accretion is another major way to grow IMBHs into SMBHs. The belief that quasars are powered by accretion at the Eddington luminosity gives this idea credibility. But over what fraction of cosmic history did these holes accrete? Some authors suppose that Eddington limited accretion is continuous since the black hole first formed (Haiman & Loeb, 2001), and thus its mass grows rapidly. Others have proposed scenarios in which accretion is negligible over the history of the BH (Madau & Rees, 2001), such that $`150\mathrm{M}_{\mathrm{}}`$ holes that formed in the first episode of star formation still populate our Galactic bulge today.

Another possibility is that the BH feeds on dark matter, in which case the halo and the IMBH evolve together. MacMillan & Henriksen (2002) proposed a model in which the halo and hole grow self-similarly, with both analytical and numerical results (MacMillan & Henriksen, 2003) supporting their idea. This model is particularly interesting because it also simultaneously explains both the spectral index for the galaxy power spectrum, $`n2`$, and the $`M_{}\sigma _{\mathrm{bulge}}^4`$ relationship observed to be universal in galaxy centers (Ferrarese & Merritt, 2000; Gebhardt et al., 2000; Tremaine et al., 2002).

All of these methods certainly contribute to the growth of massive black holes. In the next section we will use timescales to help constrain which methods could realistically dominate the process of forming SMBHs in less than a gigayear.

## 3 Timescales

First, we examine the evolution to the IMBH phase, as depicted in Figure 1. If the BH is primordial, then it can be up to $`10^5\mathrm{M}_{\mathrm{}}`$ in the very early universe. Its growth by merger or accretion then awaits the formation of matter clumps to surround the PBH. These could begin to arrive as early as photon decoupling, so timescales alone will not strongly constrain this possibility unless the PBH in question started at $`<1\mathrm{M}_{\mathrm{}}`$, as it would if it formed during the quark-hadron phase transition. In this case the BH seed has a very long way to grow, which is barely possible by accretion, as will be shown in subsection 3.3.

The simulations of population III star formation of Yoshida et al. (2003) produce some very useful information. The earliest simulated star forming region occurs at $`z=32`$, with stars forming at $`10^3\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$ per comoving $`\mathrm{Mpc}^3`$, increasing to $`0.1\mathrm{M}_{\mathrm{}}\mathrm{yr}^1\mathrm{Mpc}^3`$ at $`z=7`$, and holding constant thereafter (Hernquist & Springel, 2003).

Molecular hydrogen coolant forms in these clouds in $`30\mathrm{Myrs}`$. By $`z=23`$, the gas pressure has diminished and gas clumps begin to contract. Halos are $`5\times 10^5`$ to $`2\times 10^6\mathrm{M}_{\mathrm{}}`$ by this time, and the largest virialized star forming clouds are $`10^5\mathrm{M}_{\mathrm{}}`$. The gas takes an approximately isothermal density profile, with constant density inside a core radius of $`10\mathrm{pc}`$.

The virial temperature is a few thousand Kelvins. A simple application of the virial theorem gives a virial radius of $`45\mathrm{pc}`$, and the corresponding circular orbital speed is $`10\mathrm{km}\mathrm{s}^1`$.

This, then, is the starting point for population III stars leading to IMBHs, with star formation beginning in earnest at $`z=23`$. To convert redshifts to times in a $`\mathrm{\Lambda }`$CDM cosmology with $`\mathrm{\Omega }_0=1`$ (Peacock, 1999),

$$t_{\mathrm{age}}=\frac{2}{3}H_0^1\frac{\mathrm{sinh}^1\sqrt{|\mathrm{\Omega }_\mathrm{m}1|(1+z)^3\mathrm{\Omega }_\mathrm{m}^1}}{\sqrt{|\mathrm{\Omega }_\mathrm{m}1|}}$$
(2)
and with WMAP values inserted, this is

$$t_{\mathrm{age}}=10.8\mathrm{sinh}^1\left(\sqrt{2.7(1+z)^3}\right)\mathrm{Gyr}.$$
(3)
So gas collapses into the first stars at $`t_{\mathrm{age}}0.15\mathrm{Gyr}`$.

These stars collapse in the free fall time,

$$t_{\mathrm{ff}}=\sqrt{\frac{3\pi }{32G\rho }},$$
(4)
which is about $`30\mathrm{Myr}`$ at the densities observed in simulations by Abel, Bryan, & Norman (2000). They evolve to black holes (or disruptive supernovae) in a few million years, leaving $`0.7\mathrm{Gyr}`$ remaining to grow into a quasar power source.

Using the dynamical friction formula from Chandrasekhar (1943), quoted here from Binney & Tremaine (1987), we can find the time needed for these population III remnant BHs to sink to the center of their host halo. Treating the BHs as point particles,

$$t_{\mathrm{df}}=\frac{1.17}{G\mathrm{ln}\mathrm{\Lambda }}\frac{vr_0^2}{M_{}},$$
(5)
where $`\mathrm{ln}\mathrm{\Lambda }`$ is called the Coulomb logarithm, which characterizes the distance range perpendicular to the path of motion over which the frictional encounters occur (for a recent study of this parameter, see Spinnato, Fellhauer, & Portegies Zwart (2003)). Using the values mentioned above for the parameters in this equation, and choosing $`\mathrm{\Lambda }20`$ (this choice is common and doesn’t affect the outcome much because of the logarithm), the dynamical friction timescale relevant here is at least $`0.1\mathrm{Gyr}`$, bringing the total down to $`0.6\mathrm{Gyr}`$ available for BH growth. The coalescence time of these holes is not well known, and could be long.

If the star formation yields one or more supermassive stars, rather than multiple normal population III stars, then the evolution to a BH is faster than for normal stars (Shapiro & Teukolsky, 1983). The dynamical friction time to merge smaller stellar BHs is also avoided with SMSs, so in the SMS case, we have $`0.7\mathrm{Gyr}`$ remaining to grow a SMBH.

In the case of an early-forming dense star cluster in which runaway collisions merge the stars, the timing remains essentially the same. One waits until the cluster stars form, merge, and evolve to one or more BHs. The time available to become supermassive is still at most $`0.7\mathrm{Gyr}`$.

### 3.1 Dynamical Friction

We consider the case of a satellite subhalo halo falling into a larger halo. Each halo is modeled as a singular isothermal sphere, so that $`\rho r^2`$ and orbital velocities within a halo are constant. We seek a formulation that tracks the friction until the satellite’s central BH reaches the new center.

Smaller halos formed earlier in the CDM hierarchy, when the universe was denser, and therefore they are more tightly bound. So as a subhalo sinks, we assume it holds itself together, until it gets close to the center of the larger halo. There, which halo an orbiting particle belongs to becomes ambiguous. This occurs roughly when a particle’s distance from its original host subhalo is on order of the separation between the two halo centers. This separation distance is $`r`$, a radial coordinate extending from the center of the larger halo.

The internal speed within an isothermal halo is

$$v=\sqrt{\frac{GM}{R}}.$$
(6)
The square root in this equation is helpful because it reduces the dependence that velocities have on a halo size $`R`$, which is generally not well known.

The mass which can be associated with the sinking satellite is then

$$M_\mathrm{s}M_r+M_{}=\frac{v_\mathrm{s}^2r}{G}+M_{},$$
(7)
where $`M_r`$ is the mass of the satellite halo residing inside a distance equal to $`r`$ around it, $`M_{}`$ is the satellite subhalo’s central BH mass, and $`v_\mathrm{s}`$ is the orbital velocity inside the satellite. Choosing to use $`r`$ in this way assumes that the isothermal satellite halo extends to distances larger than some initial infall radius: $`R_\mathrm{s}r_i`$.

The angular momentum of the satellite in the larger halo is $`\stackrel{}{L}=\stackrel{}{r}\times M_\mathrm{s}\stackrel{}{v}`$, where $`\stackrel{}{v}`$ is the orbital velocity within the larger halo. Dynamical friction induces a torque $`\stackrel{}{\tau }=\stackrel{}{r}\times \stackrel{}{F}=\stackrel{}{r}\times M_\mathrm{s}(d\stackrel{}{v}/dt)`$.

Following the standard derivation presented in Binney & Tremaine (1987), the statement of $`\stackrel{}{\tau }=d\stackrel{}{L}/dt`$ leads to the equation

$$r\frac{dr}{dt}=0.428\frac{GM_\mathrm{s}}{v}\mathrm{ln}\mathrm{\Lambda }.$$
(8)
We next define 3 constants as follows: $`C_1=(0.428/v)\mathrm{ln}\mathrm{\Lambda }`$, $`C_2=C_1GM_{}`$, and $`C_3=C_1v_\mathrm{s}^2`$. Inserting the satellite mass from eq. 7, we have

$$\frac{dr}{dt}+C_2r^1=C_3.$$
(9)
Eq. 9 is a first-order, linear, ordinary differential equation, but the closed-form solution is prohibitively long and complicated due to the $`r^1`$ factor in the second term. A numerical solution readily parameterizes the sinking orbit $`r(t)`$, and a root can be found for the time $`t_{\mathrm{df}}`$ at which $`r0`$. One needs only to specify a starting radius for the satellite halo, $`r_ir(0)`$. An example orbit decay curve computed in this way is given in figure 2.

This dynamical friction model approximately doubles the time spent spiralling inward, relative to the Chandrasekhar formula for an infalling point particle, eq. 5. Colpi, Mayer, & Governato (1999) found the same result using both $`N`$-body simulations and the theory of linear response. A subsequent paper by Taffoni et al. (2003) gives an analytical approximation for NFW halos (Navarro, Frenk, & White, 1996), also finding that the decay time is increased by a factor of a few if the satellite has less than $`10\%`$ of the larger halo’s mass. In the mass range of interest for the present work, all these models essentially agree; however, since we wish to follow length and mass scales ranging down to a central BH, we have chosen to use the method described above.

### 3.2 Mergers

We first consider the viability of black hole merging as a means to grow SMBHs. As halos merge over cosmic history, their central BHs sink to the new center by dynamical friction, and subsequently coalesce into a larger BH.

The usual way to track mergers is by $`N`$-body simulation. We begin with such a simulation, and calculate the total mass of subhalos that can merge into one large halo between the epoch of the first massive BHs and the early quasar epoch. We subject this calculation to the constraint that the merging subhalos must sink to the center in the time allowed between the two epochs.

This can become a very complicated calculation, because the halo population is continuously changing over a long period of cosmic time. To make the problem tractable, we make a series of major approximations. We describe them along the way, and discuss them afterward. We demonstrate first that with our assumptions, mergers are *not* sufficient to grow SMBHs, regardless of how the BHs originally formed. Second, we show that a more accurate calculation to replace our simplifying assumptions would only strengthen that conclusion.

The argument begins with the numerical simulations of Volonteri, Haardt, & Madau (2003). They follow a $`4\times 10^{13}\mathrm{M}_{\mathrm{}}`$ halo today back in time to all its subhalos, presenting halo mass spectra at multiple redshifts. For our purposes, the results they present at $`z=5`$ and $`z=20`$ are most important, approximating the time of the first quasars ($`z=6.41`$) and the first IMBHs in halos (after $`z=23`$). We’ll address the impact of these approximations in redshift and total mass after the calculation.

The halo functions are nearly power law in the mass range of interest. For $`M`$ in $`\mathrm{M}_{\mathrm{}}`$ units, at $`z=5`$

$$\frac{dN}{dM}=10^{14}M^{\frac{7}{3}},$$
(10)
and at $`z=20`$

$$\frac{dN}{dM}=1.6\times 10^{16}M^{2.9}.$$
(11)
Consider the largest halo in existence as the period of interest begins. Imagine (conservatively) that all other subhalos are available to merge into this one, and (also conservatively) that all of these subhalos contain central BHs of $`0.1\%`$ of the parent galaxy’s luminous mass, as observed in large galaxies today (Ferrarese & Merritt, 2002). According to WMAP, the ratio of baryonic matter to total matter is $`0.17`$ (Bennett et al., 2003), so conservatively we may say that a central BH has $`0.02\%`$ of the host subhalo’s mass. Then every subhalo massive enough to dynamically sink to the center in the available time does sink, leading to a more massive combined halo.

If all the IMBHs from all the combined subhalos coalesce quickly (again, highly conservative), then the new SMBH mass is $`0.02\%`$ of the total mass of all subhalos with mass $`MM_{\mathrm{s},\mathrm{min}}`$, the minimum satellite mass needed for dynamical friction to sink a subhalo by the early quasar era.

However, during the merging process, the largest halo is gaining mass and the number of massive halos is increasing. To evade this problem, we suppose that the halo mass spectrum at the *end* of the merging process (at $`z=5`$) applies throughout the merging process. Our calculations show that in most cases, this assumtion conservatively exaggerates the final SMBH mass by more than an order of magnitude.

Using eq. 10, the total mass of the quasar’s SMBH is

$$M_{}=0.0002_{M_{\mathrm{s},\mathrm{min}}}^{\mathrm{}}10^{14}M^{\frac{4}{3}}𝑑M.$$
(12)
The lower limit, $`M_{\mathrm{s},\mathrm{min}}`$, is determined by dynamical friction as described in section 3.1. We follow the infalling halo mass, and transition to the (smaller) infalling BH mass as the two holes close in, achieving results closely consistent with other models that include mass loss but not a central BH (Colpi, Mayer, & Governato, 1999).

In the dynamical friction calculation, three parameters must be chosen. We take the internal velocity of the larger halo to be $`v100\mathrm{km}\mathrm{s}^1`$, noting that for a fixed dynamical sinking time (fixed by the time available for BH growth), the final SMBH mass obeys $`M_{}v^{1/3}`$. So choosing a larger velocity would only very slightly reduce the final SMBH mass by reducing the number of subhalos which sink in the allowed time.

The other needed quantities are radii: $`r_i`$, the initial radius from which the satellite halo begins its descent, and $`R_\mathrm{s}`$, the size of the satellite subhalo. The former is roughly the impact parameter of the colliding subhalos, and the latter indicates the subhalo’s concentration factor. Both radii will clearly vary from one encounter to the next, so we choose several representative values. The results demonstrate that any reasonable (in fact, conservative) choice of these numbers leads to the same conclusion, that merging black holes take too long.

As shown earlier, the timescale between the first IMBH and a quasar’s SMBH depends on how the BH formed. A primordial BH has the age of the universe at the time of the most distant quasar’s emission known today, or $`0.88\mathrm{Gyr}`$. Population III stars and supermassive stars have $`0.7\mathrm{Gyr}`$. If the population III stars merged within their original halo first, then they have $`0.6\mathrm{Gyr}`$ to grow.

We fix the sinking time at the appropriate value and compute the total BH mass that can be assembled by halo mergers, neglecting the time needed for coalescence. Our results are listed in Table 1, which lists a maximum SMBH mass for various choices of $`R_\mathrm{s}`$ and $`r_i`$. No choice brings $`M_{}`$ up near $`3\times 10^9\mathrm{M}_{\mathrm{}}`$, the value needed to explain the most distant known quasars.

One is led to conclude that mergers alone are not sufficient to grow the largest SMBHs in the oldest quasars. Note that although the argument above relies on many simplifying assumptions, all but one tend to overestimate the resulting SMBH mass, so the SMBH masses in table 1 should be considered extreme upper limits. Our assumptions are summarized as follows.

1. All halos contain central BHs at $`0.02\%`$ of the halo mass (this may seriously overestimate the number of BHs).

2. All relevant merging is completed by $`z=5`$, instead of today, so the number of infalling BHs is again exaggerated.

3. In the $`\mathrm{\Lambda }`$CDM world, going from redshift $`20`$ to $`5`$ allows over $`1\mathrm{Gyr}`$ for merging, so the number of more massive subhalos available to merge is artificially large. This overestimates the final SMBH mass, because larger infalling subhalos will dynamically sink faster.

4. Infall radii begin at $`r_iR_\mathrm{s}`$, which for some mergers represents a head start of multiple orders of magnitude.

5. The choice of an isothermal halo model probably exaggerates the central density, artificially speeding up the action of dynamical friction near the center.

6. SMBH coalescence is assumed to occur arbitrarily quickly, but in reality takes at least $`10^8`$ years, and probably much longer (Milosavljević & Merritt, 2001).

7. The choice to trace a $`4\times 10^{13}\mathrm{M}_{\mathrm{}}`$ halo today back in time tends to *under*estimate the number of progenitor subhalos. Extrapolating this mass linearly from Volonteri, Haardt, & Madau (2003), we find that increasing it by $`100`$ times would increase $`dN/dM`$ at $`z=5`$ by a factor of $`10`$ while approximately maintaining its slope.

In this case, the largest halo at $`z=5`$ would be $`10^{13}\mathrm{M}_{\mathrm{}}`$, so halo speeds should be comparable to those in M87 today, the giant elliptical at the center of the Virgo galaxy cluster. Radial velocities for globular clusters around M87 (Cohen, 2000) range from $`1000`$ to $`2000\mathrm{km}\mathrm{s}^1`$. These changes, using $`v1000\mathrm{km}\mathrm{s}^1`$, allow a SMBH to become almost $`10^9\mathrm{M}_{\mathrm{}}`$ in $`0.7\mathrm{Gyr}`$ from an initial radius $`r_i=R_\mathrm{s}=1\mathrm{kpc}`$, which is a factor of $`4`$ greater than our original result.

So if we attempt to correct for this seventh assumption, we approach (but do not reach) the $`3\times 10^9\mathrm{M}_{\mathrm{}}`$ hole that quasar observations require. But this is only true if we take literally the six highly conservative assumptions above.

### 3.3 Accretion

Consider first gas accretion, followed by feeding on dark matter.

*Gas accretion* is probably the simplest and most successful way to explain the growth of an IMBH into a SMBH, or even a large stellar mass BH into a SMBH. Haiman & Loeb (2001) show that if the black hole accretes at the Eddington limit continuously, from birth through detection as a quasar, it achieves the needed growth. They give the Eddingtion luminosity,

$$L_\mathrm{E}=4\pi GM_{}c\mu _\mathrm{e}m_\mathrm{p}\sigma _\mathrm{T}^1,$$
(13)
and introduce a radiative efficiency $`ϵL/\dot{M_{}}c^2`$ and a fraction of Eddington output $`\eta L/L_\mathrm{E}`$. The resulting $`e`$-folding time is

$$t_{\mathrm{acc},\mathrm{e}}=\frac{M_{}}{\dot{M_{}}}=4\times 10^7\left(\frac{ϵ}{0.1}\right)\eta ^1\mathrm{yr},$$
(14)
which equals a $`10`$-folding timescale

$$t_{\mathrm{acc},10}=9.2\times 10^7\left(\frac{ϵ}{0.1}\right)\eta ^1\mathrm{yr}.$$
(15)
For a fiducial value of $`ϵ0.1`$ and continuous Eddington-limited accretion ($`\eta 1`$), one can calculate the initial BH mass needed to grow to $`3\times 10^9\mathrm{M}_{\mathrm{}}`$ in the alloted time.

If the seed BH is primordial, it might accrete steadily from approximately the beginning of time. In this case, it has $`0.88\mathrm{Gyr}`$ to grow; that’s $`9.6`$ orders of magnitude in mass by eq. 15. So the initial seed mass must have been $`M_{\mathrm{seed}}0.8\mathrm{M}_{\mathrm{}}`$, which, interestingly, could have been formed during the cosmological quark-hadron phase transition (Jedamzik & Niemeyer, 1999).

If the seed BH is a normal population III stellar remnant, or a supermassive star’s remnant, then it had $`0.7\mathrm{Gyr}`$. It’s initial mass was then $`M_{\mathrm{seed}}70\mathrm{M}_{\mathrm{}}`$, which is completely plausible for the first generation of stars. One star of $`260\mathrm{M}_{\mathrm{}}`$ would suffice (Heger & Woosley, 2002).

In $`0.6\mathrm{Gyr}`$, the seed BH mass needs to be $`M_{\mathrm{seed}}900\mathrm{M}_{\mathrm{}}`$. This applies for the case where population III stars must sink to the center of their parent halos and merge before beginning to grow significantly by gas accretion.

In the case where the seed BH formed from the evolution of normal population III stars, the mass of the seed star is suggestively close to what one would expect. Thus this route is possible, but only if Eddington limited accretion is maintained during the entire IMBH growth process. If the seed hole evolved from a more massive SMS, then there is time to spare, and the accretion rate could have dropped for some of the growth period.

The success of this quick exponential growth calculation seems to indicate that baryon accretion is the most realistic mechanism for IMBH growth. When the BH is small, the Eddington mass accretion rate is not prohibitively large, and by $`z=6.41`$ we have very strong evidence that the quasar is accreting at its Eddington rate (Willott, McLure, & Jarvis, 2003). So it is reasonable to speculate that it has been accreting at this rate all along.

However it is worth noting that the central SMBHs in modern galaxies appear to grow by a factor of 10 on timescales of $`8\mathrm{Gyr}`$ (Merrifield, Forbes, & Terlevich, 2000). This is very slow compared to Eddington limited growth, requiring $`\eta 0.01`$. This dramatic drop in accretion rate needs explaining, especially if we are to claim that no similar drop has ever happened before $`z6`$ in all of cosmic history, as the notion of SMBH growth by gas accretion onto population III stellar corpses would mandate. This problem is only alleviated if the first stars were supermassive, or the seed BHs were primordial and significantly larger than $`1\mathrm{M}_{\mathrm{}}`$.

*Accreting dark matter* is another possibility. MacMillan & Henriksen (2002) have proposed an interesting way to account for IMBH growth in which the hole grows proportionally to the dark matter halo as matter falls in. As mentioned earlier in section 2, this assumption leads to the observed galaxy power spectrum and the observed $`M_{}\sigma `$ relation.

In numerical simulations (MacMillan & Henriksen, 2003), dark matter was found to form a self-similar region surrounding the BH, which does grow in proportion to the BH growth at the center. The formation of this region proceeds by gravitational interactions between clumps of dark matter (“particles” in the simulations) which must have formed earlier in the CDM universe. The BH grows whenever the particles cross inside the Schwarzschild radius. (The BH growth rate is thus understated in the sense that $`R_{\mathrm{sch}}`$ was used instead of the larger capture radius.)

In the absence of angular momentum, the black hole grows rapidly: $`M_{}t^4`$. But for realistic departures from pure spherically symmetric infall, in which the dark matter has angular momentum to help it resist falling into the hole, the growth of the BH is only linear in time: $`M_{}t`$ (although the $`M_{}\sigma `$ relation is recovered either way).

The simulations were not performed in standard units, and therefore cannot be immediately applied to relevant halo sizes and cosmic timescales. Here we simply note that the simulated BH growth is slow (e.g., linear in time, as compared with exponential in time for gas accretion). One therefore needs an unusually large IMBH to start from, which requires an unusually large SMS as its origin.

## 4 Conclusions

In this paper we have surveyed the state of knowledge of SMBH formation and growth in light of a difficult new data point to satisfy: a quasar at $`z=6.41`$ whose SMBH was $`3\times 10^9\mathrm{M}_{\mathrm{}}`$ when the universe was only $`880\mathrm{Myr}`$ old. We have extended the calculations of others on SMBH growth via mergers and accretion.

In the case of SMBH growth by mergers, we follow the output from a CDM simulation in order to track progenitor subhalos. We stipulate that a subhalo must sink to the center by dynamical friction after a merger within the time available for SMBH growth. This limits the mass of subhalos and therefore BHs which could have coalesced in time, thus constraining the final SMBH mass. We find that in any realistic case, mergers are incapable of growing a BH with the needed speed.

Accreting baryonic matter represents a viable growth process. For this mechanism, we impose time limits sensitive to the manner in which the seed BH first formed. The following possibilities emerge:

1. A primordial black hole that formed during the quark-hadron phase transition has accreted at the Eddington rate continuously ever since baryons have been able to cluster until at least the quasar epoch.

2. A large population III seed star with $`M260\mathrm{M}_{\mathrm{}}`$ evolved to a BH, and accreted at the maximum (Eddington) rate continuously between its formation and the quasar epoch.

3. A collection of merged population III stellar black holes formed an IMBH with $`M900\mathrm{M}_{\mathrm{}}`$, and accreted at the Eddington rate continuously through the quasar epoch.

Each of these options involves seed stars at the upper end of their predicted mass range. So although the Eddington limited gas accretion scenario is adequate to grow the holes in the time available, it is not entirely satisfying, because it requires heavy gas infall which is never significantly diminished before the quasar era, despite the relative gas suffocation which has evidently happened since. One can avoid imposing this requirement by starting with a larger BH:

4. A single SMS, or the merged remains of a dense cluster of stars, or a PBH formed well after the quark-hadron phase transition, yielded an IMBH much larger than $`100\mathrm{M}_{\mathrm{}}`$. It accreted through the quasar epoch, but did *not* need to maintain the maximum accretion rate during that time.

One other option, wherein the BH grows in proportion to its host halo by eating dark matter, is expected to yield slower BH growth. As such, normal population III stellar seeds would be inadequate, although a particularly massive SMS seed remains a possibility.

We conclude by suggesting that *some* unconventional mechanism is needed to realistically beat the clock for the early formation of SMBHs. BH growth through mergers is too slow, and growth by accreting nonstop at the Eddington limit is probably too contrived. The remedy appears to require one of the following: (a) PBHs from unusually large overdensity pockets, (b) population III starbursts well in excess of those predicted by current simulations, or (c) supermassive stars.

This work was supported by NASA under the Colorado Space Grant Consortium. We are grateful to R. N. Henriksen and J. D. MacMillan for helpful correspondence.


format: text

text:
# Condition of the occurrence of phase slip centers in superconducting nanowires under applied current or voltage

## I Introduction

Since the discovery of superconductivity it was expected that if the superconductor was subjected to a constant electric field superconductivity will inevitably be destroyed. The reason is that the superconducting electrons will be accelerated by the electric field and will reach a velocity above its critical velocity. However if the sample is short enough or if the electric field exists only in a small part of the sample such that the path over which the Cooper pairs are accelerated are sufficiently short an electric field can exist in the sample in the presence of superconductivity. Another example is the presence of an electric field in the superconducting sample which is attached to a normal metal. In this geometry the injected current from the normal metal will be converted into a superconducting current on a distance of about the charge imbalance distance ($`\mathrm{\Lambda }_Q`$) and on that scale an electric field will exist in the sample Tinkham1 . In this case the electric field is compensated by the gradient of the chemical potential $`\mu _s`$ of superconducting electrons and it does not lead to an acceleration of the superconducting condensate (see for example the book of Schmidt Schmidt ).

But there is another mechanism which allows superconductivity to survive in the presence of an electric field deep in the superconducting sample of arbitrary length. This is the phase slip mechanism. Initially this phenomena was used in order to estimate the relaxation time of superconducting current in a superconducting wire Langer . If the order parameter vanishes in one point of the wire the phase of the superconducting order parameter exhibits a jump of $`2\pi `$ at that point Langer  and as a result the momentum $`p\varphi `$ decreases by $`2\pi /L`$. So even if the electric field accelerates the superconducting electrons it does not lead to a destruction of supercondictivity because the momentum is able to relax through the phase slip mechanism.

This simple idea was understood already a long time ago. Therefore, it is very surprising that there are practically no experimental nor theoretical works studying in detail what will happen if a voltage (i.e. electric field) is applied to the superconductor. In previous works mainly the situation with applied current (i.e. the I=const regime) was studied. In the latter case the phase slip process was studied theoretically in detail (see for example the review of Ivlev and Kopnin Ivlev  and the books of Tinkham Tinkham1  and Tidecks Tidecks ). On the basis of a numerical solution of the extended time-dependent Ginzburg-Landau equations it was found that the phase slip (PS) phenomena exists in some region of currents. The lowest critical current at which a PS solution first appears in the system may be smaller than the depairing Ginzurg-Landau current density and as a result it leads to hysteresis in the I-V characteristics in the I=const regime Ivlev ; Tidecks .

The I=const regime was also studied experimentally in a number of works Ivlev ; Tidecks  starting from the paper of Meyer and Minnigerode Meyer1 ; Meyer2 . The most characteristic effect of the phase slip mechanism is the appearance of a stair like structure in the current-voltage characteristics. In Ref. Skocpol  a simple phenomenological model was built in order to quantitatively describe this feature. It was proposed that every step in the I-V characteristic is connected with the appearance of a new phase sip center (PSC) that increases the resistivity of the sample by a finite value. In this model, it was conjectured that in a region of the sample with size of about the coherence length $`\xi `$ fast oscillations of the order parameter occurs in time which produces normal quasi-particles. Because of the finite time needed to convert normal electrons into superconducting electrons Tinkham1 ; Schmidt  there is a region of size of about $`\mathrm{\Lambda }_Q`$ near the phase slip center where the electric field and the normal current are different from zero. It means that the chemical potential of superconducting $`\mu _s`$ and normal electrons $`\mu _n`$ are different near PSC and their difference is proportional to the charge imbalance Q in a given point of the superconductor Tinkham1 ; Skocpol . Because usually $`\mathrm{\Lambda }_Q\xi `$ it is possible to neglect the $`\xi `$-sized region with oscillating $`|\psi |`$ and consider that phenomena as a time-independent process. Such a model gives a finite voltage, which is connected with every phase slip center, which is equal to $`V_{PSC}=2\mathrm{\Lambda }_Q\rho _n(I\beta I_c)/S`$ with $`\rho _n`$ the normal resistivity, S the cross-section of the current carrying region, $`I_c`$ is the critical current and $`\beta <1`$ is a phenomenological parameter. In the experiment of Dolan and Jackel Dolan  the distribution of $`\mu _s`$ and $`\mu _n`$ near a PSC were measured which fully supported the idea of Ref. Skocpol  that a difference exist between $`\mu _s`$ and $`\mu _n`$ near the phase slip center.

Despite numerous theoretical and experimental works the physical conditions under which phase slip phenomena can exist is still not clear. In a review Ivlev  on this subject it was claimed that phase slip phenomena are connected with the presence of a limiting cycle in the system which leads to such type of oscillations. But it was not explained how and why it leads to the phase slip phenomena.

Based on our previous work Vodolazov2  on the time-dependent Ginzburg-Landau equations (the investigated systems were superconducting rings in the presence of an external magnetic field) we know that systems which are governed by such equations exhibit two relaxation times. One is the relaxation time of the phase of the order parameter $`\tau _\varphi `$ and the other is the relaxation time of the absolute value of the order parameter $`\tau _{|\psi |}`$. In Ref. Vodolazov2  it was established that phase slip processes can occur in such systems when roughly speaking $`\tau _\varphi <\tau _{|\psi |}`$. In the present paper we discuss this question in the context of superconducting wires in the presence of an applied current or voltage.

To our knowledge there exists only a single theoretical work in which a superconducting wire in the presence of an applied voltage was studied Malomed . The authors used the simple time-dependent Ginzburg-Landau equations and found that the behavior of the system is very complicated and strongly depends on the length of the wire and the applied voltage. Neither a detailed analysis nor any physical interpretation of their results was presented. In a recent letter Vodolazov1  we presented our preliminary theoretical and experimental results on the dynamics of the superconducting condensate in wires under an applied voltage. It turned out that in this case the I-V characteristics exhibit a S-shape which was explained by the appearance of phase slip centers in the wire and their rearrangement in time. In the present paper we will present more details and extend our previous work to the situation in which defects are present, and we investigate the effects of boundary conditions and an applied magnetic field.

The paper is organized as following. In Sec. II we present our theoretical results and give the conditions for the existence of phase slip centers when current (Sec. IIa) or voltage (Sec. IIb) is applied to the superconducting wire. In Sec. III we show our experimental results and in Sec. IV we compare theory and experiment.

## II Theory

We study the current-voltage characteristics of quasi-one-dimensional superconductors using the generalized time-dependent Ginzburg-Landau (TDGL) equation. The latter was first written down in the work of Ref. Kramer1

$`{\displaystyle \frac{u}{\sqrt{1+\gamma ^2|\psi |^2}}}\left({\displaystyle \frac{}{t}}+i\phi +{\displaystyle \frac{\gamma ^2}{2}}{\displaystyle \frac{|\psi |^2}{t}}\right)\psi =`$
$`=(\mathrm{i}𝐀)^2\psi +(1|\psi |^2)\psi .`$ (1)
In comparison with the simple time-dependent Ginzburg-Landau equation where $`\gamma =0`$ it allows us to describe a wider current region (with a proper choice of the parameters $`u`$ and $`\gamma `$) where a superconducting resistive state exists and gives us a wider temperature region in which Eq. (1) is applicable Kramer1 ; Kramer2 ; Watts-Tobin . The inelastic collision time $`\tau _E`$ for electron-phonon scattering is taking into account in the above equation through the temperature dependent parameter $`\gamma =2\tau _E\mathrm{\Delta }_0(T)/\mathrm{}`$ ($`\mathrm{\Delta }_0=4k_BT_cu^{1/2}(1T/T_c)^{1/2}/\pi `$ is the equilibrium value of the order parameter).

Eq. (1) should be supplemented with the equation for the electrostatic potential

$`\mathrm{\Delta }\phi `$ $`=`$ $`\mathrm{div}\left(\mathrm{Im}(\psi ^{}(\mathrm{i}𝐀)\psi )\right),`$ (2)
which is nothing else than the condition for the conservation of the total current in the wire, i.e. $`\mathrm{div}𝐣=0`$. In Eqs. (1,2) all the physical quantities (order parameter $`\psi =|\psi |e^{i\varphi }`$, vector potential $`𝐀`$ and electrostatical potential $`\phi `$) are measured in dimensionless units: the vector potential $`𝐀`$ and momentum of superconducting condensate $`𝐩=\varphi 𝐀`$ is scaled by the unit $`\mathrm{\Phi }_0/(2\pi \xi )`$ (where $`\mathrm{\Phi }_0`$ is the quantum of magnetic flux), the order parameter is in units of $`\mathrm{\Delta }_0`$ and the coordinates are in units of the coherence length $`\xi (T)`$. In these units the magnetic field is scaled by $`H_{c2}`$ and the current density by $`j_0=c\mathrm{\Phi }_0/8\pi ^2\mathrm{\Lambda }^2\xi `$. Time is scaled in units of the Ginzburg-Landau relaxation time $`\tau _{GL}=4\pi \sigma _n\lambda ^2/c^2=2T\mathrm{}/\pi \mathrm{\Delta }_0^2`$, the electrostatic potential ($`\phi `$), is in units of $`\phi _0=c\mathrm{\Phi }_0/8\pi ^2\xi \lambda \sigma _n=\mathrm{}/2e\tau _{GL}`$ ($`\sigma _n`$ is the normal-state conductivity). In our calculations we mainly made use of the bridge geometry boundary conditions: $`|\psi (L/2)|=|\psi (L/2)|=1`$, $`\phi (L/2)=0`$, $`\phi (L/2)=V`$ and $`\psi (L/2,t+dt)=\psi (L/2,t)e^{i\phi (L/2)dt}`$. We chose these boundary conditions because at low temperatures the normal current is converted to a superconducting one due to the Andreev reflection on a distance of about $`\xi _00.18\mathrm{}v_F/k_BT_c`$ near the $`SN`$ boundary Tinkham1 ; Hsiang . It means that practically there is no injection of quasi-particles from the normal material to the superconductor and hence we can neglect the effect of charge imbalance near the S-N boundary. The bridge geometry boundary conditions models this situation. The parameter $`u`$ is about $`5.79`$ according to Ref. Kramer1 . We also put $`A=0`$ in Eq. (1,2) because we considered the one-dimensional model in which the effect of the self-induced magnetic field is negligible and we assume that no external magnetic field is applied.

### II.1 Constant current regime

Lets consider first the more simple case when a constant external current is applied to the sample. In such a case it was theoretically found Kramer1 ; Kramer2  that the system exhibits hysteretic behavior. If one starts from the superconducting state and increases the current the superconducting state switches to the resistive superconducting or normal state at the upper critical current density $`j_{c2}`$ which, in a defectless sample, is equal to the Ginzburg-Landau depairing current density $`j_{GL}=\sqrt{4/27}j_0`$. Starting from the resistive state and decreasing the current it is possible to keep the sample in the resistive state even for currents up to $`j_{c1}<j_{c2}`$ (which we call the low critical current). For $`j_{c1}<j<j_{c2}`$ such a state is realized as a periodic oscillation of the order parameter in time at one point of the superconductor Kramer1 ; Kramer2 . When the order parameter reaches zero in this point a phase slip of $`2\pi `$ occurs. This is the reason why such a state is now called a phase slip state and this point a phase slip center (PSC). Using results obtained in our earlier work Vodolazov2  we claim that the value of $`j_{c1}`$ depends on the ratio between the two characteristic times in the sample: the phase relaxation time of the order parameter $`\tau _\varphi `$ and the relaxation time of the absolute value of the order parameter $`\tau _\psi `$ in the region (with size of about $`\xi `$) where the oscillations of the order parameter occurs.

Having written Eq. (1) for the dynamics of the phase and the absolute value of the order parameter in a quasi-one dimensional wire of length L (-L/2$`<`$s$`<`$L/2)

$`u\sqrt{1+\gamma ^2|\psi |^2}{\displaystyle \frac{|\psi |}{t}}`$ $`=`$ $`{\displaystyle \frac{^2|\psi |}{s^2}}+|\psi |(1|\psi |^2p^2),`$ (3a)
$`{\displaystyle \frac{\varphi }{t}}`$ $`=`$ $`\phi {\displaystyle \frac{\sqrt{1+\gamma ^2|\psi |^2}}{u|\psi |^2}}{\displaystyle \frac{j_n}{s}},`$ (3b)
it is easy to estimate both these times. Indeed, from Eq. (3a) it directly follows that

$$\tau _{|\psi |}u\sqrt{1+\gamma ^2|\psi |^2}u\gamma (\mathrm{for}\gamma 1).$$
(4)
To determine $`\tau _\varphi `$ we need to know how fast the phase changes over the region (with size of about $`\xi `$) at which the order parameter oscillates. Our numerical calculation shows that the time-average of the derivative is $`\varphi /t=C\theta (s)`$ (with C(s)=const and $`\theta (s)`$ the theta function). Because the time-averaged electrostatic potential $`\phi `$ changes over a distance $`\mathrm{\Lambda }_Q1`$ for $`\gamma 1`$ (where $`\mathrm{\Lambda }_Q^2=\sqrt{1+\gamma ^2|\psi |^2}/u|\psi |^2\gamma /u`$ is the decay length of the normal current density and the charge imbalance $`Q`$) and that the order parameter is about unity already at $`s=\pm \xi `$, we can estimate this constant as $`C\mathrm{\Lambda }_Q^2j_n/sj/\mathrm{\Lambda }_Q`$. Consequently we find that

$$\tau _\varphi \frac{1}{\mathrm{\Lambda }_Qj}.$$
(5)
The product $`\mathrm{\Lambda }_Qj`$ is roughly the voltage drop $`V_{PSC}`$ ’produced’ by the phase slip center Skocpol . Thus the time change of the phase of the order parameter in the phase slip center is inversely proportional to the voltage drop over the whole structure (and $`C=V_{PSC}/2`$). This is a consequence of the fact that the time-averaged electro-chemical potential of the superconducting electrons $`\mu _s^e=\varphi /t`$ does not change in space and its jump near the phase slip center is equal to the voltage drop of this structure in the sample. In terms of the language used by Schmid and G. Schön Schmid , $`\tau _{|\psi |}`$ is the ’longitudinal’ time and $`\tau _\varphi `$ decreases with increasing ’transverse’ time because $`\tau _\varphi 1/\mathrm{\Lambda }_Q1/\sqrt{\tau _Q}`$. In Ref. Vodolazov2  it was found that the phase slip events occur periodically in time when $`\tau _\varphi \tau _{|\psi |}`$. This gives us an estimation for the lower critical current

$$j_{c1}\frac{1}{\tau _{|\psi |}\mathrm{\Lambda }_Q}.$$
(6)
By varying the parameter $`\gamma `$ we can change both $`\tau _{|\psi |}`$ and $`\mathrm{\Lambda }_Q`$ and hence we can vary $`j_{c1}`$. In Fig. 1 we show the voltage in the sample, the normal current density and the absolute value of the order parameter in the phase slip center averaged over the time as a function of the external current for two different values of $`\gamma `$. In our simulations we started from the superconducting state with $`j<j_{GL}`$. At $`j>j_{GL}`$ the system jumps instantaneous to the resistive state with finite voltage self1 . Than we decrease the external current and at $`j<j_{c1}`$ the system transits to the purely superconducting state with $`V=0`$. We should emphasize that at $`jj_{c1}`$ the voltage jumps by a finite value (this result is qualitatively different to the results of Ref. Kramer1 ; Watts-Tobin  where the authors found that $`V0`$ at $`jj_{c1}`$.). It means that there is a non-infinite maximal oscillation period for the order parameter in the phase slip center. We believe that the finite voltage jump $`\mathrm{\Delta }V`$ or finite period of oscillations is directly connected with the threshold condition $`\tau _\varphi \tau _{|\psi |}`$ for activation of regular phase slip processes in the constant current regime and it means that $`\mathrm{\Delta }V1/\tau _{|\psi |}`$.

Before going further we should stress here that the above condition for the existence of a phase slip process $`\tau _\varphi \tau _{|\psi |}`$ is a rather rough estimate. Indeed, Eqs. 3(a,b) are a coupled system of equations, and besides $`\tau _{|\psi |}`$, it explicitly (see Eq. (4)) depends on the value of $`|\psi |`$. However, the above condition allows us to explain the general qualitative properties of the phase slip processes (including the existence of $`\mathrm{\Delta }V`$, $`j_{c1}`$ and its dependence on $`\gamma `$ and $`u`$) and predict new features which will be discussed below.

Far enough from $`j_{c1}`$ the dependence of $`V(j)`$ on the current is close to linear. When the current increases, $`\tau _\varphi `$ decreases (see Eq. (5)) and hence the time $`\tau _{PSC}`$ between two phase slips also decreases and the order parameter has less time for recovering at the phase slip center. That is the reason why the time averaged voltage increases (as $`V=2\pi /\tau _{PSC}`$), the averaged order parameter $`|\psi (0)|`$ decreases and the fraction of the normal current $`j_n(0)/j`$ increases with increasing external current. It is interesting to note that from the phenomenological Skocpol-Beasley-Tinkham (SBT) Skocpol  model it follows that $`j_n(0)/j=1\beta j_{c1}/j`$ which qualitatively resembles the dependence shown in Fig. 1(b).

Let us now discuss the effect of defects on the I-V characteristic. This question was considered previously in Ref. Kramer3  for two different models of defects: local variation of the critical temperature and a local variation of the mean free path. We will repeat these calculations partially and interpret it in terms of a competition between $`\tau _\varphi `$ and $`\tau _{|\psi |}`$. In addition to the first type of defect we will also study the effect of the variation of the cross-section of the wire.

The first type of defect is the inclusion of a region in the superconductor which suppresses $`T_c`$ and the order parameter becomes lower than the equilibrium value $`\mathrm{\Delta }_0`$ even in the absence of any external current. This is modelled Kramer3  by introducing the term $`\rho (s)\psi `$ in the right hand side (RHS) of Eq. (1). In the present calculation we choose $`\rho (s)=\rho _0\theta (0.5|s|)`$. The larger $`\rho _0`$ the more the order parameter is suppressed in the center of the wire. Firstly, such a defect leads to a decrease of the upper critical current $`j_{c2}`$. Secondly, it decreases the lower critical current $`j_{c1}`$. Indeed, when we introduce a defect, the RHS of Eq. (3a) decreases at the spatial position where the order parameter oscillations and hence $`|\psi |`$ needs more time to change in that spatial position. Thus this type of defect leads to an increase of the relaxation time of the order parameter in the region of the defect (as an indirect prove of this we found a decreasing $`\mathrm{\Delta }V`$ with increasing ’strength’ of the defect). If the size of the defect is smaller than $`\mathrm{\Lambda }_Q`$ we can neglect its effect on the relaxation time of the phase of the order parameter. Finally, we can conclude that the ’stronger’ the defect the smaller the value of the lower critical current $`j_{c1}`$.

In Fig. 2(a) we present our numerical results for two different values of $`\rho _0`$. With increasing ’strength’ of the defect the upper and lower critical currents decrease and start to merge. As a consequence the hysteresis in the I-V characteristics will disappear when the defect is sufficiently strong.

The second type of defect is one for which we have a local decrease of the cross-section of the wire. In this case we expect that $`\tau _{|\psi |}`$ will not be influenced. The situation with $`\tau _\varphi `$ is more complicated because a large part of the superconductor with size of about $`\mathrm{\Lambda }_Q`$ participates in the formation of this time. Two limiting cases can be distinguished. If the defect is much smaller than $`\mathrm{\Lambda }_Q`$ we can neglect its effect on $`\tau _\varphi `$ and hence we will have the same lower critical current as for the case of an ideal wire. In the opposite case of a defect with size much larger than $`\mathrm{\Lambda }_Q`$ we only have to take into account the increased current density in the part of the wire where we have a smaller cross-section and apply Eq. (5). In this case the current $`j_{c1}`$ is decreased by a factor $`D_{av}/D_l`$, where $`D_{av}`$ is the average cross-section and $`D_l<D_{av}`$ is the local reduced cross-section.

Another important question is the value of the upper critical current $`j_{c2}`$. If the size of the defect is larger than $`\xi `$ then the proximity effect from adjacent parts near the defect will be small and $`j_{c2}`$ will decrease by a factor $`D_{av}/D_l`$ (in the opposite case the current $`j_{c2}`$ is almost defect independent). Therefore, for a defect with length $`\xi <l\mathrm{\Lambda }_Q`$ the I-V characteristic may be reversible with a proper choice of the parameters (like in the case of a local suppression of $`T_c`$). In Fig. 2(b) we show the results of our numerical calculations and we obtain a qualitative agreement with the above physical arguments. When the length of the defect is smaller than $`\mathrm{\Lambda }_Q`$ then the lowest critical current density practically does not change. In the opposite case $`j_{c1}`$ decreases by a factor of $`D_{av}/D_l`$. The upper critical current density changes considerably only if the length of the defect exceeds $`\xi `$. Unfortunately, we are not able to consider the case for which the length of the defect is simultaneously much smaller than $`\mathrm{\Lambda }_Q`$ and much larger than $`\xi `$ because of computational restrictions. For example, when we increase $`\gamma `$ by a factor of two the time of calculations increases also by a factor of two but the ratio $`\mathrm{\Lambda }_Q/\xi `$ increases roughly only by a factor of $`\sqrt{2}`$ for $`\gamma 1`$. The reason is that when we increase $`\gamma `$ we need to take a smaller time step (see Eq. (1)) which increases the computation time.

When the wire has a length less than $`\mathrm{\Lambda }_Q`$ then this will affect the distribution of the normal current density in the wire and hence the time $`\tau _\varphi `$ and the lower critical current $`j_{c1}`$. Indeed from our bridge boundary conditions it follows that $`j_n/s(\pm L/2)=0`$ (see Eq. 3(b)). Then taking into account that $`\varphi /t=V\theta (s)/2`$ we can easily obtain, in the limit $`\mathrm{\Lambda }_Q\xi `$, that

$$j_n(0)=\frac{V(j)}{2\mathrm{\Lambda }_Q}\frac{1}{\mathrm{tanh}(L/2\mathrm{\Lambda }_Q)},$$
(7)
where the current $`j_n(0)`$ is the average normal current in the phase slip center.

The voltage jump $`\mathrm{\Delta }V2\pi /\tau _{|\psi |}`$ at $`j_{c1}`$ should not change with varying $`L`$ (at least for $`L\xi `$ when proximity effect does not have an effect on $`\tau _{|\psi |}`$). From Eq. (7) it follows then that the normal current density increases with increasing length. But $`j_n(0)`$ should always be smaller than the full current $`j`$. It implies that the critical current density $`j_{c1}`$ will increase with increasing length of the wire as $`j_{c1}1/\mathrm{tanh}(\mathrm{L}/2\mathrm{\Lambda }_\mathrm{Q})`$ in order to keep the ratio $`j_n(0)/j`$ constant. Fig. 3 illustrates the dependence of $`j_{c1}`$ on the length of the wire which we obtained on the basis of a numerical solution of Eqs. (1,2). Unfortunately, in our calculations we are not able to use very large values of $`\gamma `$ and the maximal value of $`\mathrm{\Lambda }_Q`$ was $`4.1\xi `$ for $`\gamma =40`$. But nevertheless, we found that $`j_{c1}`$ increases with decreasing wire length. We should add here that we also found that $`\mathrm{\Delta }V`$ also decreases a little bit. We connect it with a small change in $`\tau _{|\psi |}`$ due to the small length of the wire and hence the increased effect of the boundaries.

The more pronounced effect of the finite length of the wire on the value of $`j_{c1}`$ is found for the case that we used N-S boundary conditions: $`\psi (\pm L/2)=0`$ and $`\phi /s(\pm L/2)=j`$. These boundary conditions are approximately valid for samples at temperatures close to $`T_c`$. It is easy to show that in this case

$$j_n(0)=\frac{V(j)}{2\mathrm{\Lambda }_Q}\frac{\mathrm{tanh}(L/2\mathrm{\Lambda }_Q)}{11/\alpha \mathrm{cosh}(L/2\mathrm{\Lambda }_Q)},$$
(8)
where $`\alpha =j_n(0)/j<1`$ (see Fig. 1(b)). In this case $`j_{c1}`$ also increases with decreasing wire length (see inset in Fig. 3). But besides there is a finite length $`L_0`$ for which $`j_n(0)\mathrm{}`$. It implies that the phase slip process is not possible in wires with length $`L<L_0`$. In such wires the system goes from the superconducting state directly to the normal state.

Not only the finite length of the sample is able to change the value of $`j_{c1}`$. If we apply a magnetic field parallel to the length of the wire it will suppress the order parameter in the sample. Our analysis shows that if the diameter of the wire $`d`$ is less than $`2\xi `$ and if we can neglect screening effects ($`\lambda >\xi `$) then the distribution of the order parameter will be uniform along the cross-section of the wire. The order parameter depends on H as

$$|\psi |^2=1(H/H_c)^2$$

with $`H_c2.9\mathrm{\Phi }_0/\pi \xi d`$. This behavior is very similar to the behavior of a thin plate in a parallel magnetic field Ginzburg ; Douglass  or a thin and narrow ring in a perpendicular magnetic field Vodolazov3 . In all cases the transition to the normal state is of second order and the vorticity in the wire will be equal to zero due to the small cross-section of the sample.

Because the order parameter practically does not depend on the radial coordinate of the wire we can use the one-dimensional model in order to study the response of the system on the applied current. In order to take into account the suppression of $`|\psi |`$ by magnetic field ($`H`$) we add to the RHS of Eq. (1) the term $`(H/H_c)^2\psi `$. In some respect it is similar to the way we introduced the first type of defect in our wire. We can expect that $`\tau _{|\psi |}`$ will increase with increasing $`H`$ (because the ’strength’ of the ’defect’ increases). But because the magnetic field suppresses the order parameter everywhere in the sample it also leads to an increase of $`\mathrm{\Lambda }_Q`$, because in the TDGL model $`\mathrm{\Lambda }_Q1/\sqrt{|\psi |}`$. It is clear that both these processes should decrease $`j_{c1}`$. The strongest mechanism is connected with the change in $`\tau _{|\psi |}`$. Indeed it is easy to estimate that $`\tau _{|\psi |}1/(1(H/H_c)^2)`$ and $`\tau _\varphi (1(H/H_c)^2)^{1/4}`$. Even if we take into account that the parameter $`\gamma `$ may decrease with increasing H (because in non-zero magnetic field there is another pair-breaking mechanism and instead of $`\tau _E`$ we should useSchmid ; Kadin1  $`\tau _E/\sqrt{1+2\tau _E\tau _s}`$ with $`\tau _s=\mathrm{\Delta }(T=0,H=0)/\mathrm{}(H/H_c(T=0))^2`$ ) it does not lead to an increase of $`j_{c1}`$ with an increase of $`H`$ because $`\tau _{|\psi |}`$ changes faster than $`\tau _\varphi `$ even in this case.

In the above model it is easy to show that $`j_{c2}(H)=\sqrt{4/27}(1(H/H_c)^2)^{1.5}`$ for the case of an uniform wire. In Fig. 4 we present the results of our numerical calculations. We found that both $`j_{c1}`$ and $`j_{c2}`$ decreases with increasing magnetic field and at some $`H^{}`$ they practically merge. Unfortunately, it is quite difficult to find an analytical expression for the dependence of $`j_{c1}(H)`$ like we had for $`j_{c2}`$. The reason is the complicated behavior of the dynamics of $`\psi `$ in the phase slip center.

The variation $`\mathrm{\Lambda }_Q`$ with increasing $`H`$ was obtained experimentally in Ref. Kadin1 . To compare with the theory the authors of Ref. Kadin1  used the expressions found in the work of Schmid and Schön Schmid  which are valid in the limit $`TT_c`$. It is interesting that in Ref. Kadin1  it was found that expressions of Schmid and Schön are quantitatively valid even far from $`T_c`$. From this observation we may hope that the present theoretical results are also valid over a wider temperature range than only near $`T_c`$.

The main conclusion which follows from the change of $`j_{c1}`$ with decreasing wire length and/or applying magnetic field is that there exist a critical length $`L^{}`$ (or critical field $`H^{}`$) below(above) which the current $`j_{c1}`$ becomes equal to $`j_{c2}`$. It implies that for wires with lengths $`L<L^{}`$ and/or fields $`H>H^{}`$ there will be no jump in the voltage in the current-voltage characteristic and the I-V curve will be reversible. It will also result in the absence of a S-behavior in the V=const regime (see section below).

### II.2 Constant voltage regime

In our earlier work we found that the I-V in the $`V=const`$ regime exhibits an S-like behavior Vodolazov1  and furthermore for low voltages there is an oscillatory dependence of the current on the applied voltage (see Fig. 5). As was shown in Ref. Vodolazov1  these properties are connected with the existence of two critical currents $`j_{c1}`$ and $`j_{c2}`$. Here, we will discuss the characteristic voltages $`V_1`$, $`V_2`$ (see Fig. 5) and their dependence on the length of the sample.

If we apply a voltage $`V`$ to the wire of length $`L`$ then in the sample an electric field $`E=V/L`$ exist which will accelerate the superconducting electrons. When the current density approaches $`j_{c2}`$ phase slip centers will spontaneously appear in the sample. As a result the momentum of the superconducting electrons (and hence the current) will decrease by $`2\pi /L`$ after each phase slip event. When the current density decreases below $`j_{c1}`$ the phase slip process will no longer be active and the applied electric field will be able to accelerate the superconducting condensate. This process leads to periodic oscillations in time of the current in the sample.

We can divide the period of oscillations (at least at voltages $`V<V_1`$) in two parts. First, the longest part is the one during which the condensate is accelerated by the electric field till the moment reaches $`jj_{c2}`$. The second part we call the transition period $`T_{tr}`$. The transition period also consists of two parts: the time needed for the phase slip processes (which is proportional to the number of phase slip events and hence the length of the wire) and the time $`T_0`$ for the decay of the order parameter from $`\sqrt{2/3}`$ (when $`jj_{c2}`$) till the first phase slip event and for the recovering of $`|\psi |`$ back to $`\sqrt{2/3}`$ from zero after the last phase slip event (see Fig. 6). The minimal number of phase slip events which occur during the transition time is determined by the internal parameters of the superconductor ($`j_{c1}`$) and its length Vodolazov1  which is given by

$$\mathrm{N}_{\mathrm{min}}=\mathrm{Nint}\left((\mathrm{p}_\mathrm{c}\mathrm{p}_{\mathrm{c1}})(\mathrm{L}/2\pi +1)\right),$$
(9)
where $`p_{c1}`$ is the smallest real root of the equation $`j_{c1}=p_{c1}(1p_{c1}^2)`$ and $`\mathrm{Nint}(\mathrm{x})`$ returns the nearest integer value.

At $`V<V_1`$ the time-averaged current $`j`$ increases (with oscillations) and in the range $`V_1<V<V_2`$ it decreases with increasing voltage. The lowest minimal current in the latter region is $`j_{c1}^V`$ which depends on the length of the system (see Fig. 5 and Ref. Vodolazov1 ). In Fig. 7 we present the dependence of the above voltages on the length of the wire. The explanation for their different behavior is the following. At a voltage $`V_1`$ the period of the oscillation $`T=2\pi N/V_1`$ (N is the number of phase slip events during the transition time) of the current becomes of order $`2T_{tr}`$. The time $`T_0`$ does not depend on the length and the voltage at $`VV_1`$. So, we can estimate $`V_1`$ as

$$V_1\frac{2\pi N}{2T_{tr}}=\frac{\pi N}{T_0+\tau _{PSC}N}=\frac{\pi \alpha L}{T_0+\tau _{PSC}\alpha L},$$
(10)
where $`\tau _{PSC}`$ is the average time between two phase slip events and $`\alpha `$ is the coefficient which depends on $`p_{c1}`$ (see Eq. (9)) and hence on the specific superconductor. When $`T_0\tau _{PSC}\alpha L`$ the voltage $`V_1`$ becomes practically independent of the length. Because the time $`T_0`$ depends on the relaxation time of the absolute value of the order parameter the length of the wire $`L_{sat}`$ at which $`V_1`$ saturates depends on the internal parameters of the superconductor.

The voltage $`V_2`$ decreases with increasing length of the sample because the lower critical current $`j_{c1}^V`$ decreases. It is interesting to note that with accuracy of finding $`V_1`$ the saturated value of the lower critical voltage coincides with the voltage jump $`\mathrm{\Delta }V`$ at $`j=j_{c1}`$ in the constant current regime. Because the minimal value of $`V_2`$ is equal to $`\mathrm{\Delta }V`$ for an infinitely long wire we may conclude that with increasing wire length the I-V characteristic in the range of voltages $`(V_1,V_2)`$ approaches the horizontal line.

Defects, magnetic field, short length of the sample etc. will decrease the hysteresis and consequently the currents $`j_{c1}`$ and $`j_{c2}`$ approaches each other. If the difference between them will be small enough the S-shape of the I-V characteristic at the constant voltage regime changes to the usual monotonic behavior (see Fig. 8). This happens because when the voltage approaches $`V_1`$, the maximal current during the period of the current oscillation $`T`$ will be about $`j_{c1}`$ (because $`j_{max}`$ cannot substantially exceed $`j_{c2}`$ but $`j_{c2}`$ is already close to $`j_{c1}`$ in this case).

We also would like to discuss how the change of the boundary conditions will change the shape of the I-V characteristic. If we apply the N-S boundary conditions the I-V curve in the voltage driven regime also exhibits an S-behavior but without oscillations in the current at small voltages (see Fig. 9). In this case there will always be an inevitable voltage drop near the boundaries connected with the current in the wire through the relation $`V_{NS}j\mathrm{\Lambda }_Q`$ (because the electric field and the normal current decays on a scale of the charge imbalance length $`\mathrm{\Lambda }_Q`$ near the boundaries Tinkham1 ). Near the boundary the nonzero electric field is compensated by the term $`\mathrm{\Lambda }_Q^2E/s`$ (see Eq. 3(b) - and $`E=j_n`$ in our units) and the superconducting electrons are not accelerated by this field. The situation is very similar to the case when we inject a current in the wire. When the current generated by this voltage reaches $`j_{c1}`$ the phase slip process becomes possible in the system. But like the case of the current driven regime the superconducting state can be stable(metastable) until the current reaches $`j_{c2}`$ at $`Vj_{c2}\mathrm{\Lambda }_Q`$. From our estimates it follows (see section below) that the fluctuations of the order parameter are very important in our samples. Therefore, in our calculations we introduce at some moment of time (at fixed voltage) a phase slip center in the center of the wire and checked if it will survive or not. As a result we obtained the I-V characteristics as presented in Fig. 9. When the voltage is less than some critical value ($`V^{}`$) the phase slip process decays in time and the resulting current in the wire is time-independent. The whole voltage drop occurs near the boundaries. At $`V>V^{}`$ (and $`j>j_{c1}`$) the dynamics of the condensate in the wire will be similar to the case considered above for $`V>V_1`$.

The reason for this is as follows. The current density in the wire at small voltages will always be less than $`j_{c1}`$ and $`j_{c2}`$ due to the relation $`jV_{SN}/\mathrm{\Lambda }_Q=V/\mathrm{\Lambda }_Q`$. When the current density reaches $`j_{c1}`$ the phase slip process becomes possible in the sample. But such a phase slip process leads to a finite voltage. As a consequence the voltage drop at the boundaries will sharply decrease when a phase slip center is created. But it implies that the full current density will also decrease and consequently will become less than $`j_{c1}`$. We can conclude that the phase slip process may survive in such a type of superconductor only if the applied voltage will be roughly equal to $`\mathrm{\Delta }V(j)`$ plus the voltage drop near the boundaries $`V_{SN}`$ necessary for the creation of a current larger than $`j_{c1}^V`$.

Concluding the theoretical part, we present in Fig. 10 the I-V characteristics in the I=const and V=const regimes at different temperatures close to $`T_c`$. With decreasing temperature the range of currents, where there is a S-behavior increases and the voltages $`V_1`$ and $`V_2`$ increases. It resembles the experimental results presented in Vodolazov1 ; Michotte  and in the subsequent Section III.

## III Experiment

Measurements of the I-V characteristics were done on single Pb nanowires Vodolazov1 ; Michotte  with typical diameter about $`50nm`$ and length $`2250\mu m`$ (see Table I). At first we would like to discuss the rather wide resistive transitions and strong magneto-resistance effect (see Figs. 11-12) in our samples. We may claim, taking into account the small diameter of our nanowires, that the effect of the thermo-activited Langer  and the quantum-activated Giordano ; Tinkham2  phase slip phenomena is very strong in our samples. It is easy to see that with increase of the diameter of the nanowire the width of the resistive transition decreases. Our estimations, based on Giordano expressions, Giordano  showed that for sample A even at $`T=0`$ the number of quantum phase slip events should be about $`10^5`$ per second. That is the reason why we did not observe in our experiment any hysteresis in the current driven regime (see below). But this rate is not large enough for destroying the S-behavior in the voltage driven regime. Indeed, the period of oscillations of the current for Pb is about $`10^9`$s at $`VV_1`$ and $`T=0`$. It means that only at temperatures close to $`T_c`$ the fluctuating PSC’s will interrupt the internal temporal oscillations in the order parameter and ruin the S-behavior. In other words the system can be in a state with $`j>j_{c1}`$ only during a time which is less than the time between two phase slips. This results in the coincidence of the I-V characteristics both in the current and the voltage driven regimes at temperatures close to $`T_c`$ (see Fig. 3 in Ref. Michotte ). We would like also to mention the small difference in the critical temperatures which implies that all our samples have practically the same value of the superconducting gap.

As can be seen in Fig. 12, the wider the transition width in $`R(T)`$, the smaller the magnetic fields where the resistance starts to increase with increasing magnetic field.

Of course, our nanowires are not free from imperfections. Although our fabrication technique produce samples which are quite regular (see Ref. Michotte  for Scanning Electron Microscopy (SEM) picture), we may distinguish two kinds of imperfections. Firstly, there can be deviations from an ideal cylindrical shape. Indeed, although the nanopores of the membrane in which the nanowires are grown are designed to be as regular as possible, there may exist a smooth variation in diameter along the length of the nanowire, which we estimated to be not more than 10 $`\%`$. However, in addition to these smooth variations, it may happen that a defect in the membrane leads to a constriction in the measured nanowire. In this case, the diameter can be locally reduced quite importantly and this results in a lower critical current and a larger critical magnetic field. This is in fact what we observed in sample C. In addition to such shape imperfections, we note also that the parameters of our nanowires strongly vary from sample to sample - see table I (we estimated the coherence length using the expression $`H_c=2.9\mathrm{\Phi }_0/(\pi \xi d)`$ ($`\xi _1`$) and $`H_c=\mathrm{\Phi }_0/(2\pi \xi ^2)`$ ($`\xi _2`$) for the critical field). The reason for this difference in resistivity is related to the second kind of imperfections, which is structural disorder that is formed during the electrodeposition of these nanowires inside the nanopores. Indeed, as shown in Ref. Dubois , these nanowires are polycrystalline and inevitably contain structural defects such as dislocations, twins, etc. These two kinds of imperfections results in differences in the I-V characteristics (Fig. 13). For example for samples A and B we observed two jumps in the voltage which we explain by the appearance of two successive phase slip centers in the nanowire but for sample C we found only one jump in the voltage. Indeed, due to the presence of a constriction in this sample, heating may drastically affect the behavior of this sample, precipitating the return to the normal state. Unfortunately sample D was broken during the measurements and we could not measure its I-V characteristic.

We found that jumps in the voltage $`\mathrm{\Delta }V`$ are practically the same for all three samples (see Figs. 13-14). It agrees with the theory presented in Sec. II. Indeed all the three samples have almost the same $`T_c`$ and hence the same superconducting gap. It is naturally also to suppose that the time $`\tau _E`$ is almost the same for all our samples. We can expect that the parameters $`\gamma `$ and $`u`$ are the same for all our samples and consequently the relaxation time $`\tau _{|\psi |}`$ and the ratio $`\mathrm{\Lambda }_Q/\xi `$ are the same. This automatically leads to the invariance of $`\mathrm{\Delta }V`$ on $`\rho _n`$. From the results of Sec. IIb follows that the voltage $`V_2`$ decreases with increasing (in $`\xi `$) nanowire length (see Fig. 6). For our shortest sample A: L$`6001400\xi `$ (see Table I). It means that in our nanowires the voltage $`V_2`$ already reaches the minimal value $`\mathrm{\Delta }V`$ and hence $`V_2`$ is independent of $`\rho _n`$ for our nanowires. The voltage $`V_1`$ depends not only on the length of the nanowire but also on the time change of the order parameter during the transition period $`T_0`$ (see Eq. (9)) and hence the voltage $`V_2`$ may reach $`\mathrm{\Delta }V`$ at shorter lengths than $`V_1`$ if the time $`T_0`$ is large enough (in Fig. 6 the opposite situation is presented - at first $`V_1`$ reaches the saturated value).

If the current density is uniformly distributed over the cross-section of the sample then the oscillations of the order parameter will be in phase along the cross-section and in this case the $`\mathrm{\Delta }V`$ will not depend on the size of the sample. It is a direct consequence of the fact that the time-averaged chemical potential of the superconducting electrons should be equal and constant on both sides of the phase slip center or phase slip line/surface.

Our results show if one uses the Skocpol-Beasley-Tinkham (SBT) Skocpol  model for the estimation of $`\mathrm{\Lambda }_Q`$ one should be very careful. Indeed, those authors replaced the normal current density in the PS center by the expression $`I_n(0)=(I\beta I_{c1})`$. As a result the derivative $`dI_n(0)/dI=1`$ becomes current independent. But in general the time derivative may be large than unity (see Figs. 1(a,b)). Secondly, if the length of the sample is comparable with $`\mathrm{\Lambda }_Q`$ we should replace $`2\mathrm{\Lambda }_Q`$ by $`2\mathrm{\Lambda }_Q\mathrm{tanh}(L/2\mathrm{\Lambda }_Q)`$ (see Eq. (7)). As a result the differential resistance in samples with $`L\mathrm{\Lambda }_Q`$ may be larger, in general, than the normal one ours2 . This case corresponds to our samples ($`R_{dif}320`$ Ohm, $`R_{dif}351`$ Ohm and $`R_{dif}456`$ Ohm for samples A,B,C respectively after the first jump in voltage). But from the SBT model follows that $`R_{dif}R_n`$ (the equality sign holds for samples with $`L\mathrm{\Lambda }_Q`$).

Unfortunately we do not know the actual dependence $`I_n(0)(I)`$ for our samples. If we use the values obtained from the SBT model ($`\mathrm{\Lambda }_Q=12.6\mu m`$, $`48.5\mu m`$ and $`30.3\mu m`$ for samples A,B and C, respectively) which is qualitatively understandable for longer samples, we have a too small number of PSC (compare samples A and C). And in this case, the question arises why $`\mathrm{\Lambda }_Q`$ changes so much. Probably, the SBT model gives us only the correct order of magnitude for $`\mathrm{\Lambda }_Q`$ which is only useful as an estimate.

Finally, we present our results on the influence of an applied magnetic field on the I-V characteristics. We limit ourselves to data for sample A in the current driven regime (see Fig. 15). These results already support our theoretical predictions of Section IIb. It is evident that the lower critical current density decreases with increasing applied magnetic field. From some range of H-values only one jump in the voltage exists. We explain it by an increase of $`\mathrm{\Lambda }_Q`$ at relatively large magnetic fields and hence there is a lack in space for the coexistence of two phase slip centers in the nanowire. At high magnetic fields the order parameter is strongly suppressed by H and the effect of quantum phase-slip fluctuations becomes more pronounced. This is the reason for a smoothing of the I-V characteristics at high magnetic fields.

## IV Discussion

In conventional superconductor near T<sub>c</sub> the N-S boundary conditions (in the sense mentioned in Sec. IIa) are valid Hsiang . For $`T0`$ the bridge boundary conditions are more applicable due to Andreev reflections at the border between the normal metal and the superconductor Andreev ; Hsiang . At intermediate temperatures we expect a mixture of the N-S and bridge boundary conditions. It means that part of the voltage will drop at the boundaries and part in the sample. The situation in our experiment is even more complicated because in our two-point measurements there is also a voltage drop at the contacts. Probably, this will not allow us to observe the oscillations of the current in the voltage driven regime. Another reason is that our samples are very long compared to $`\xi `$ (even for our 22 $`\mu m`$ sample $`L1000\xi `$) and consequently the amplitude of the oscillations is very small. But nevertheless our calculations in both limiting cases of bridge and N-S boundary conditions predicts an S-shape of the I-V characteristic.

We found that the type of boundary conditions are not so important for the process of nucleation of phase slip centers if the length of the sample is much larger than $`\mathrm{\Lambda }_Q`$. However, the difference in the process of the conversion of superconducting electrons to normal ones and vice versa at the N-S boundary at various temperatures plays a crucial role for the creation of phase slip centers in shorter nanowires. It turned out that for similar parameters (nanowire’s length, coherence length, superconducting gap, $`\tau _E`$) phase slip centers appear in the superconducting wire at smaller currents for the case of the bridge geometry boundary conditions.

All our theoretical results are strictly speaking only valid near $`T_c`$, which is the temperature region where Eqs. (1,2) are quantitatively correct. Nevertheless, experimental results supports our prediction on the competition of two relaxation times in the creation of a phase slip center even far from $`T_c`$. Indeed, an applied parallel magnetic field decreases the lower critical current (compare Fig. 15 and Fig. 4). Jumps in the voltage $`\mathrm{\Delta }V`$ turned out to be almost an independent function of the disorder (of the resistance of the sample) as it follows from theory. Unfortunately, it is rather difficult to check the dependence of $`j_{c1}`$ on the length of the nanowire using our technique because every preparation of a new sample leads to a different level of disorder and hence different values for $`\xi `$ and $`\mathrm{\Lambda }_Q`$.

Therefore, we expect only quantitative differences in the dependence of $`j_{c1}(T,H)`$, $`j_{c2}(T,H)`$, as compared with our theoretical results. Some qualitative differences (see Ref. Baratoff ) in the dynamic of the order parameter at the phase slip center or the creation of charge imbalance waves Kadin2  cannot affect the main properties of our theoretical results because it does not influence the existence of the two different critical currents: $`j_{c1}`$ and $`j_{c2}`$. It may lead to quantitative differences in the dependence of $`\tau _\varphi `$ and $`\tau _{|\psi |}`$ on the microscopic parameters of the superconductor.

Finally, we would like to discuss other mechanisms which, for our geometry, may lead to an S-behavior of the I-V in the V=const regime. The first is heating Volkov . In order to explain the double S-structure for our samples A, B by this mechanism we have to assume that heat dissipation and heat evacuation have a very complicated and non trivial dependence on temperature. We do not know any mechanisms which can lead to such adependence in our case. Besides we did not any observe hysteresis in the current driven regime which is an inevitable property of that mechanism if the I-V characteristic would have a S-shape in the voltage driven regime. For these reasons we believe that heating is not responsible for the observed behavior.

Secondly, if the nanowire contains a normal region the I-V characteristic will exhibit an S-behavior due to multiple Andreev reflection in the SNS structure Nicolsky . We do not have any indication for this process from our R(T) and R(H) measurements which shows that our samples are quite homogeneous. Furthermore the structure would occur at voltages much smaller than $`\mathrm{\Delta }/e`$. In our case the S-behavior is seen for $`V>\mathrm{\Delta }/e`$ \- see Fig. 13 (for Pb, $`\mathrm{\Delta }(0)1.4meV`$). Andreev reflection on the N-S boundaries may give rise to a zigzag shape of the I-V in the V-const regime Sols  but this effect is negligible for our samples with $`L\xi _0`$.

###### Acknowledgements.

This work was supported by IUAP(P5/1/1), GOA (University of Antwerp), ESF on ”Vortex matter”, the Flemish Science Foundation (FWO-Vl) and the ”Communauté Française de Belgique” through the program ”Actions de Recherches Concertées”. D.V. is supported by DWTC to promote S $`\&`$ T collaboration between Central and Eastern Europe. We thank the POLY lab. at UCL for the polymer membranes. S. Michotte is a Research Fellow of the FNRS.