text:
# Vortex dynamics and states of artificially layered superconducting films with correlated defects

## I Introduction

The early recognition that the most interesting cuprate superconductors consist of alternating layers of superconducting and insulating material has encouraged an enormous effort to model the complex vortex structure and dynamics in such layered systems. Blatter et al. (1994) Difficulties achieving strong pinning of vortices in these layered systems, especially those in which the Josephson coupling between adjacent superconducting layers is weak, lends special urgency to the problem, for weakly pinned vortices lead to resistive losses. The research in the cuprates has been carried out on several fronts, with a search for enhanced pinning in these layered systems running parallel to advances in the modelling and understanding of the thermodynamic vortex states in material with various pinning properties. Thus the former work has concentrated on the introduction of correlated defects which traverse many layers, while more fundamental work has demonstrated a first-order solid-liquid vortex transition in samples which are largely free of pinning defects. Schilling et al. (1996); Roulin et al. (1996) The less dramatic phase transition in a pinned vortex fluid is usually regarded and treated as second order, Fisher (1989); Fisher et al. (1991); Nelson and Vinokur (1993) though there continue to be doubts as to whether there is a thermodynamic phase transition at all in this situation. Strachan et al. (2001); Landau and Ott (2002)

It is fully recognized that the state and dynamics of vortices in these material depend strongly on the degree of coupling between the superconducting layers, so that, for example, vortex motion is much more easily controlled in YBa2Cu3O7-δ (Y123) than in the less strongly coupled layers in Bi2Sr2CaCu2O8+δ (Bi2212). There are other examples of cuprates in which the anisotropy and coupling is intermediate between these two, but the lack of good single crystals spanning a full range of anisotropies has encouraged a focus almost exclusively on these two systems. The range of correlated defects is similarly limited, and has depended largely on heavy ion damage, and to a lesser extent on grain boundaries. Although we now understand much of the vortices in the materials, it is important to test that understanding in a greater range of layering characteristics and defect geometries.

One option available for a fuller study is based on the preparation of artificially layered superconductors. There have indeed been a number of studies on layered films, including multilayers of Nb/Ge, Ruggiero et al. (1980, 1982); Koorevaar et al. (1994) Mo/Ge, White et al. (1991) and Y123/PrBa2Cu3O7. Brunner et al. (1991); Triscone et al. (1994); Radovan and Ziemann (1999); Yang et al. (1999) For some years we have contributed to the work with experiments performed on Ta0.3Ge0.7/Ge multilayers. The system is attractive on two important counts: (i) the constituents form a consistent amorphous structure across the entire range of thicknesses required for the study, and (ii) all of the important measurements are at temperatures below $`2.17`$ K. The first of these assures that the layers do not change their properties as their thickness is varied, and furthermore encourages the formation of smooth layers. The second allows measurements to be performed while the films are immersed in superfluid He; the excellent heat extraction ensures a temperature very close to the bath temperature even at the high currents required to drive the material into the normal state.

In a number of earlier papers we have described the behavior of homogeneous TaGe films, Ruck et al. (1997) of multilayers with varying coupling strengths, Trodahl et al. (1996); Ruck et al. (1997) and of non-layered films with correlated defects. Abele et al. (1999) More recently we reported some aspects of the behavior of a film in which correlated defects were introduced into strongly coupled layered films. Engel et al. (2001a) The present paper completes the picture as regards the effects of correlated defects in a layered environment, and seeks to answer specific questions in relation to the existence of a second order phase transition for varying interlayer coupling strengths and the vortex dynamics in the solid vortex state in the presence of correlated disorder.

The study we report is based on thin films prepared by vacuum deposition. Alternating layers are formed by simply alternating the composition of the vapor stream, and correlated defects threading the films at about $`20`$$`30^{}`$ to the normal are introduced by depositing onto a substrate held at an oblique angle to the vapor stream. The vortex dynamics are then investigated by $`I`$/$`V`$ measurements, for vortices aligned both along the extended defects and at a large angle to them. The measurements are performed at many magnetic fields and temperatures in order to sample the full phase diagram.

The remainder of this paper is organized as follows: We start in Sec. II with experimental housekeeping followed in Sec. III with an investigation of $`R`$ vs. $`\theta `$ to demonstrate the presence of correlated defects and to establish their orientation. In Sec. IV we present measurements of $`R`$ vs. $`T`$ (below $`H_{c2}`$ but far from the solid vortex state) to estimate activation energies, followed by $`IV`$-characteristics (V) to establish phase diagrams (VI), and detailed investigations of the dynamics in the glass phase (VII). In each of the last four of these there will be two subheadings contrasting strongly to weakly coupled samples. Finally (VIII) the conclusions.

## II Samples and Experimental Details

The amorphous multilayered films in this study were prepared by vapor deposition onto glass substrates in a vacuum chamber with a base pressure of less than $`10^8`$ mbar. Ta was evaporated from an electron-gun heated source, and Ge from a tungsten boat. The two evaporation rates were independently controlled to within $`5`$% of their average with quartz crystal deposition monitors. A Ta concentration of about $`35`$% was chosen to maximize $`T_c`$. Trodahl et al. (1996) The layering characteristics were achieved by periodically placing a shutter across the Ta vapor stream. Tilting the substrates with respect to the incoming vapor flux resulted in the formation of columns with regions of lower density in between them. van Kranenburg and Lodder (1994) In previous studies we have demonstrated the effectiveness of these regions as extended pinning sites. Abele et al. (1999); Engel et al. (2001a) The layering and defect parameters were established by transmission electron microscopy (TEM), and their chemical composition has been measured by Rutherford backscattering spectroscopy (RBS). Engel et al. (2001b) Table 1 lists these details for the films contributing to this report.

All three films for which we report results here were deposited at an angle of $`40^{}`$ to the substrate normal, leading to columns oriented near $`25^{}`$ to the normal. We have previously demonstrated that films deposited at this angle show significantly enhanced pinning, especially for magnetic fields well aligned with the columnar orientation. Engel et al. (2001a) Two of the films, labelled as C40 and J40, show relatively strong Josephson coupling between adjacent layers, due to their Ge layers being well below the $`5`$ nm thickness above which that coupling becomes weak. Note that J40 has superconducting TaGe layers that are about $`30`$% thinner than C40, which we show below does not affect the qualitative behavior of their superconductivity. The third film, P40, has Ge layers well above $`5`$ nm, ensuring that the coupling is weak. It will be demonstrated below that the strength of this coupling has a profound effect on the vortex hopping energy, the phase diagram and the vortex dynamics in the vortex glass phase.

The films were prepared for the four-terminal $`I/V`$ measurements by scribing a path $`2`$ mm wide and $`5`$ mm long, defining the current to flow perpendicular to the columnar microstructure (Fig. 1).

They were inserted in the cryostat in an orientation which allowed the magnetic field to be applied at any arbitrary direction in the plane perpendicular to the current, including along the defect structure. Below we signify that geometry, with the field parallel to the extended coplanar defects, by the letter C. As discussed above the columns are directed approximately $`25^{}`$ from the normal, so that the symmetric orientation relative to the normal, labelled here the AC direction, placed the field about $`50^{}`$ to the columns. We will compare measurements made in these two field directions below.

The samples were immersed directly in 4He, with a magnetic field applied by either a Varian electromagnet or a superconducting solenoid capable of reaching $`7`$ T. The software-controlled current was supplied by a Keithley model 224 precision current source and the voltage measured with a Keithley model 182 nanovoltmeter. More complete details of the measurements systems are found elsewhere. Engel et al. (2001a); Engel (2001)

## III Angular Dependence of Resistance

Careful determination of the C-direction is very important for the interpretation of the subsequent conductivity measurements. Since the samples had to be broken to obtain cross-sectional views from the TEM investigations which were available only after the completion of the conductivity measurements, the columnar orientation was determined from resistance measurements as a function of $`\theta `$, where $`\theta `$ is the angle enclosed between the magnetic field direction and the film normal. For the right combination of temperature and applied field $`R(\theta )`$-curves show a pronounced dip near the columnar angle in addition to the anisotropy due to the layered structure of the films. Figure 2

shows a set of $`R(\theta )`$-curves for sample J40 taken at constant temperature and varying magnetic fields. The C-directions for all the samples were determined from such $`R(\theta )`$-curves as the angle of the local minimum in resistance; results are given in Table 1. Comparing these angles to those derived from TEM images, a systematic deviation towards larger angles is observed. This trend is also seen in other samples not presented here. The deviation is most likely due to the development of vortex staircases because of the competing pinning by the coplanar defects and the intrinsic pinning by the insulating germanium layers. Studies on Y123 with competing extended pinning sites Silhanek et al. (1999) have shown that such staircase structures can experience strong pinning.

Zero field transition temperatures $`T_c(0)`$ and upper critical fields $`H_{c2}(T)`$ were determined from measurements of the resistance $`R`$ versus temperature $`T`$ in applied fields ranging from zero to $`2.5`$ T. The zero-field data could be well described using the expression for fluctuation conductivity in three dimensional (3D) superconductors. Schmidt (1968); Maki (1968); Thompson (1970) In-field $`R(T)`$ data were analyzed according to the scaling approach of Ullah and Dorsey Ullah and Dorsey (1990, 1991) and assuming the upper critical field to be linear in $`T`$. The upper critical field was determined for fields in the C and AC-direction, and the extrapolated values at zero temperature were identical in the two orientations within the accuracy of the scaling analysis. The extrapolated $`H_{c2}(0)`$ values were used to estimate the coherence length $`\xi _{ab}(0)`$ according to the following relation from Ginzburg-Landau (GL) theory: Ginzburg and Landau (1950)

$$H_{c2}=\frac{\mathrm{\Phi }_0}{2\pi \xi _{ab}^2},$$
(1)
where $`H_{c2}`$ is the upper critical field perpendicular to the $`ab`$-planes, $`\mathrm{\Phi }_0`$ is the flux quantum $`h/2e`$, and $`\xi _{ab}`$ is the coherence length within the $`ab`$-planes. For moderate anisotropy111The anisotropy ratio $`\gamma =H_{c2}/H_{c2}`$ was not determined for these films, but estimated from previous studies to be less than $`1/10`$. the systematic error in $`\xi _{ab}`$ as a result of using $`H_{c2}`$ in the C-direction amounts to $`<5`$%. Blatter et al. (1992) The penetration length $`\lambda _{ab}(0)`$ at zero temperature can be calculated using the relation of Kes and Tsuei:Kes and Tsuei (1983)

$$\lambda _{ab}(0)=1.05\times 10^3\left[\rho _N(0)/T_c(0)\right]^{1/2}$$
(2)
with $`\rho _N(0)`$ the extrapolated normal state resistivity at zero temperature. Values for these important parameters are given in Table 2. It is worth noting the exceptionally large GL parameters $`\kappa \lambda /\xi 200`$ for these multilayered films.

## IV Activation Energy for TAFF

It is generally accepted that immediately below $`T_c(H)`$ the vortex system is in a liquid state with no long range order and vanishing shear modulus. The liquid state can be further divided into a free flux flow region, which is extremely narrow for the superconductors under investigation, and a thermally activated flux flow (TAFF) region. The pinning potential which hinders the free flow of vortices is caused by a combination of point-like defects222caused mostly by oxygen impurities and the correlated defects of the columnar microstructure. In general, the temperature dependence of the resistance can be described by the following exponential for activated behavior

$$R\mathrm{exp}\left(\frac{U_A(H,T)}{k_BT}\right),$$
(3)
with $`k_B`$ the Boltzmann constant and $`U_A(H,T)`$ the field and temperature dependent activation energy.

The temperature and field dependence of the activation energy for plastic motion of the flux lines has been calculated within the limits of weak collective pinning Geshkenbein et al. (1989); Vinokur et al. (1990, 1991)

$$U_{\text{pl}}\frac{T_cT}{\sqrt{H}}.$$
(4)
Replacing $`U_A`$ in Eq. 3 with $`U_{\text{pl}}`$ the $`T`$ and $`H`$-dependence can be separated and the activation energy as determined from Arrhenius-plots of $`R(T)`$-curves taken at constant $`H`$ should follow a simple $`H^{1/2}`$ power law. In two-dimensional (2D) or quasi-2D superconductors another mechanism for TAFF is possible, the movement of unbound dislocations. The expression for the activation energy is given as Feigel’man et al. (1990)

$$U_{\text{dis}}\left(T_cT\right)\mathrm{ln}\left(\frac{H_0}{H}\right),$$
(5)
with $`H_0`$ a characteristic field of the order of $`H_{c2}(0)`$. This kind of dissipation mechanism is topologically forbidden in three-dimensional superconductors since only dislocations with zero net flux can be realized. The $`T`$ and $`H`$-dependence can be separated again when inserting expression (5) for $`U_A`$ in Eq. 3 and in thin films or strongly layered superconductors the activation energy for TAFF should obey a logarithmic field dependence in case $`U_{\text{dis}}<U_{\text{pl}}`$.

### IV.1 Strongly Coupled Layers

Figure 3(a) shows the activation energies in the two strongly coupled samples in a double-logarithmic plot.

As a guideline, the solid line represents the results for a homogeneous, unstructured Ta$`_\text{x}`$Ge$`_{\mathrm{\hspace{0.17em}1}\text{x}}`$-film fitted to a field-dependence according to Eq. 4. Ruck (1998) For comparison the $`\mathrm{ln}H`$ least-square fit in the lower graph (b) is shown as the dotted curve. For C-aligned fields in the range from $`0.1`$ up to about $`2.0`$ T both samples represented in (a) have a significantly increased activation energy for TAFF due to increased pinning by extended defects. The enhanced pinning is seen to slowly diminish for fields exceeding $`0.2`$ to $`0.3`$ T. If this is identified with the matching field, for which the density of vortices is equal to the number of pinning defects, this corresponds to an average spacing of strong pinning defects of approximately $`100`$ nm. The $`1/\sqrt{H}`$-dependence is approximately preserved for fields in the AC-direction, although the activation energies for $`H>0.2`$ T are consistently higher than for the reference film. We understand this as a combined effect of the layering and the coplanar defect structure in the present films. Engel et al. (2001a); Engel (2001) Overall, the field-dependence of the activation energies of the strongly-coupled films are consistent with anisotropic, three dimensional superconductors.

### IV.2 Weakly Coupled Layers

In contrast with the results for the strongly coupled films described in the previous paragraph, the weakly coupled film P40 is more complex. The data for the AC field direction, aligned far from the columns, follows the 2D prediction of Eq. 5. To emphasize this behavior the activation energies for this film are displayed in a linear-log plot in Fig. 3(b). For comparison the $`1/\sqrt{H}`$-relation from graph (a) is plotted as the curved, dotted line in this semi-logarithmic plot. Especially at lower field values the measured activation energies deviate from the power-law behavior, but can be fitted very successfully using an expression according to Eq. 5 with a logarithmic dependence on the applied magnetic field (solid line). The fitting parameter, the characteristic field $`H_0=10.0`$ T, is comparable to the upper critical field $`H_{c2}(0)=7.3`$ T as required by the theory. Thus, the most likely process leading to dissipation is the movement of flux-carrying dislocations, possible only in 2D or quasi-2D, layered superconductors.

The interpretation of the activation energies for fields aligned with the extended defect structure in sample P40 is more difficult. Firstly, it has to be noted that activation energies for fields smaller than $`0.8`$ T are significantly higher compared to AC-oriented fields. At higher fields on the other hand, the activation energies are equal for the two field orientations within experimental errors, in contrast to the strongly coupled samples where the difference between the orientations remains significant up to at least $`2.0`$ T. Although the field-dependence can be roughly described by the power-law behavior of Eq. 4, with increased $`U_A`$ over a certain field-range similar to C40 and J40, the values obtained below $`0.1`$ T are significantly below the expectations of Eq. 4 and approach the results from AC-oriented fields, just as at high fields.

It appears unlikely that the differences between the C and AC aligned fields can be understood based on strictly uncoupled quasi-2D layers, for in this picture the extended defects act as simple point defects centered on the superconducting layers. We propose instead a picture involving dimensional transitions, suggested by theoretical descriptions of pancake liquids Bulaevskii et al. (1998) and by experiments on Bi2212 containing columnar defects. Sato et al. (1997); Kosugi et al. (1997); Morozov et al. (1999) At small magnetic fields and corresponding high temperatures pancake vortices are trapped by the extended defects, but independently in each layer. As the field is increased the correlations in the direction of the applied field increase and the pancakes recouple to form vortex lines. In this intermediate field range the activation energies for TAFF are much higher for C than for AC-aligned fields. At even higher magnetic fields the magnetic intra-layer interactions begin to dominate and the pancakes decouple again.

## V Current-Voltage Characteristics

The state and dynamics of the vortex fluid can be studied directly by measuring the current-voltage relation in samples exposed to a magnetic field. The interpretation of these data are based on the recognition that the product of current and field determines the force applied to the vortices, while the measured voltage is proportional to the product of the field and the vortex flow velocity. Thus the usual presentation of the data, as a plot of voltage versus current, is equivalent to a plot of the speed of the response (flow velocity) vs. force for the vortex fluid. In this section we describe $`I`$-$`V`$ results for a range of field and temperatures, probing the vortex dynamics over the portion of the phase diagram accessible within the limits of our experimental arrangements.

One of the most fundamental issues, whether a vortex state with true zero resistance exists for $`T>0`$ in the mixed phase of dirty type-II superconductors, is still unresolved. Strachan et al. (2001); Landau and Ott (2002) The simple question whether the barriers for flux creep diverge in the limit of zero applied current below a field dependent transition temperature is inherently difficult to answer experimentally. The most popular theoretical model predicting a second order phase transition from a vortex liquid into a low-temperature solid vortex phase with zero resistance is based on an analogy to magnetic spin glass systems. Huse et al. (1992) Depending on the anisotropy of the superconductor and the character of the defects (point-like or extended) one differentiates between 2D and 3D vortex glasses (VG) Fisher (1989); Fisher et al. (1991) and Bose-glasses (BG). Nelson and Vinokur (1992)

At the second order phase transition a characteristic length and corresponding time diverge algebraically and these define a set of critical glass exponents. A scaling ansatz between electric field $`E`$ and current density $`J`$ is obtained from which a number of relations between applied current $`I`$ and voltage $`V`$ can be derived for both the vortex liquid and solid state. Those relations which were important for the analysis of the $`IV`$-characteristics presented here will be repeated. Neglecting some specific differences between the VG and BG theories these relations can be summarized conveniently using a parameter $`D`$, with $`D=2,3`$ for quasi-2D and 3D VG and $`D=4`$ for BG. Right at the glass transition temperature $`T=T_g`$ the voltage increases as a simple power-law of the applied current

$$VI^{\left(z+1\right)/\left(D1\right)}.$$
(6)
As the transition is approached from above ($`T>T_g`$) the resistivity $`\rho `$ vanishes as

$$\rho |TT_g|^{\nu \left(z+2D\right)},$$
(7)
and below the transition temperature $`IV`$-curves show an exponential relation with zero linear resistance

$$V\mathrm{exp}\left[c\left(I_c/I\right)^\mu \right].$$
(8)
Relations (7) and (8) cross over to critical power-law behavior above a crossover current given by

$$I_\text{x}|TT_g|^{\nu \left(D1\right)}.$$
(9)
$`\nu `$ and $`z`$ are the static and dynamic critical exponents, respectively; they should be field-independent and equal for films which belong to the same universality class. $`\mu `$ is the glass exponent and will be discussed in detail below, $`c`$ and $`I_c`$ are material-dependent parameters. Appropriately rescaled sets of $`IV`$-curves should collapse onto just two universal functions above and below the transition temperature using the critical exponents $`\nu `$ and $`z`$ as obtained from above relations. In order to use this condition as an additional check of the analysis $`(V/I)|TT_g|^{\nu (zD+2)}`$ is plotted versus $`I|TT_g|^{\nu (1D)}`$.

However, the collapse of $`IV`$-curves onto an upper and lower branch as the only criterion for a vortex glass transition is not suitable. Due to the experimentally limited current and voltage resolution a collapse may be achieved for a range of glass exponents. Voss-de Haan *et al.* Voss-de Haan et al. (1999, 2000) have demonstrated the collapse of a set of $`IV`$-curves for deliberately limited current ranges but varying parameters. It is therefore mandatory to achieve a collapse of the $`IV`$-curves and to fulfill Eq. 6 \- 9 simultaneously for the same set of parameters.

### V.1 $`IV`$-Curves in Strongly Coupled Layers

As an example, Fig. 4 shows the analysis of $`IV`$-curves taken for sample J40 and nine magnetic field magnitudes in the range $`0.1`$ to $`1.3`$ T applied in the AC direction at constant temperature.

The temperature range spans typically about $`300`$ mK on either side of the glass transition temperature $`T_g`$ in intervals of $`50`$ mK. The main graph is a double-logarithmic plot of the linear resistance $`R`$ at low currents above $`T_g`$ (left ordinate) and the crossover current $`I_\text{x}`$ (right ordinate) marking the onset of nonlinear behavior in the vortex liquid phase above $`T_g`$ versus $`TT_g`$. $`I_\text{x}`$ was defined here as the current for which the resistance exceeded twice the linear resistance. To construct the plot the value of $`T_g`$ was first estimated as the temperature which showed power-law behavior (Eq. 6). Subsequently these $`T_g`$ values were adjusted until the data in Fig. 4 followed straight lines. A least-squares fit to the data allows one to determine the exponents in Eqs. 7 and 9.

Knowing these exponents one can check the collapse of the $`IV`$-curves as shown in the inset of Fig. 4 in a double-logarithmic plot of the rescaled current $`I^{}=I|TT_g|^{\nu (1D)}`$ and voltage $`V^{}=(V/I)|TT_g|^{\nu (zD+2)}`$. Often this procedure already leads to a very satisfactory collapse. If not, the parameters are slightly adjusted and then the scaling of the linear resistance and crossover current checked and if necessary corrected, and so forth until all conditions for a glass transition are fulfilled within the experimental limits. To extract the glass exponents $`\nu `$ and $`z`$ one needs to know which of the different glass models is appropriate, i.e. $`D=2,3`$ or $`4`$. As discussed in an earlier paper Engel et al. (2001a) we use $`D=3`$ for the strongly coupled samples C40 and J40 for AC-aligned fields and the BG model, $`D=4`$, when the external fields were applied in the C-direction. The resulting critical glass exponents listed in Table 3 are very similar to those reported for vortex glasses in HTSCs. Grigera et al. (1998); Voss-de Haan et al. (1999); Ammor et al. (2000); Klein et al. (2000) We stress that our values are field-independent over at least one order in magnitude of the applied field.

Strachan *et al.* argue that the above analysis is still not sufficient evidence for a transition to a vortex state with vanishing linear resistance. Strachan et al. (2001) According to Eqs. 6 \- 9 $`IV`$-curves plotted in double-logarithmic fashion must change their curvature from positive above $`T_g`$ to negative below the transition temperature. This property is preserved in the universal functions, which are specified by the scaled data sets in the inset of Fig. 4. Whereas individual data sets above $`T_g`$ describe the upper branch including the characteristic kink to the horizontal part, rescaled current-voltage curves below $`T_g`$ describe only very short sections and only a set of $`IV`$-curves for different $`TT_g`$ give a full representation of the lower branch as in the inset. This problem arises because of the extremely rapidly decreasing voltage with decreasing current and the experimental voltage resolution allows $`IV`$-curves to be taken over only a relatively small current range. As a consequence an apparent collapse of $`IV`$-curves can be achieved for a range of scaling parameters as reported by Voss-de Haan Voss-de Haan et al. (1999, 2000) and Strachan *et al.*Strachan et al. (2001) demand the detection of positive curvatures above the glass transition and negative curvatures below as a necessary condition for a vortex glass transition.

Due to the extremely rapidly vanishing voltage level at low temperatures for decreasing currents, our data do not fulfill this experimentally extremely demanding condition for a second order phase transition. However, data taken for the weakly coupled film includes $`IV`$-characteristics at constant field which are obviously lacking any sign change in curvature as will be demonstrated below. Based on these observations and thermodynamic measurements on Y123 \[Ref. (Roulin et al., 1998)\] we believe that an interpretation of our data as a second order phase transition into a vortex glass at low temperatures is justified.

### V.2 $`IV`$-Curves in Weakly Coupled Layers

The $`IV`$-curves measured for the weakly coupled film P40 differ qualitatively from those of the strongly coupled films. The differences are most easily seen by plotting the logarithmic gradient $`\text{d}(\mathrm{ln}V)/\text{d}(\mathrm{ln}I)`$ versus the applied current, as displayed in Fig. 5 for a selection of situations.

For reference, Fig. 5(a) includes data from a strongly coupled film showing the high-temperature nearly-ohmic behavior with a logarithmic derivative very close to unity for all currents. As the temperature is lowered the derivative is close to unity (ohmic) at low and high currents, passing through a maximum (super-ohmic) at intermediate currents. At even lower temperatures, below the glass transition temperature, the low-current ohmic behavior is missing; the continuing rise of the logarithmic derivative with decreasing current signals a voltage which falls more rapidly than any power law with decreasing current. Note that this presentation identifies immediately the problem of using $`IV`$-curves to establish the presence of super-ohmic behavior at vanishing currents as a signature of a solid vortex phase. At low currents, where the higher temperature curves show an unambiguous tendency toward ohmic behavior, the voltage level at lower temperatures is too small to be measured.

We now compare this behavior to that of film P40, shown for C-oriented fields in Figs. 5(b) and (c). At the lowest applied field, Fig. 5(b), the logarithmic derivatives are qualitatively similar to those in a strongly coupled film. Their behavior at lower temperatures continues to rise with decreasing current, with no suggestion of a low-current ohmic dependence in the accessible current range. In a higher field, however, that behavior is entirely different, with even the low-temperature curves showing a tendency to fall at the lowest currents. There is clearly no pinned solid vortex state at all in the weakly coupled film at these higher C-oriented fields. Interestingly, similar data for a field of $`0.3`$ T shows nearly horizontal lines at all temperatures, corresponding to the power-law behavior expected at a phase transition.

Further insight into the vortex states in this film is gained by application of the $`IV`$ scaling analysis. At the low fields of Fig. 5(b) the analysis is successful, but at higher fields only the upper branch ($`T>T_g`$) can be formed. That high-field scaling analysis returns a transition temperature within our measurement range. The lower temperatures cannot be accommodated as in the glass phase, nor can they be accommodated within a model in which the glass transition is assumed to be at temperatures below our lower limit. It appears that a *would-be* transition is sensed in the vortex liquid phase, but that the transition never, in fact, takes place.

The scaling analysis suggested by Eqs. 7 and 9 is still available for these results, for it relies on data only above the glass transition. The resulting scaling of $`R`$ and $`I_\text{x}`$ is shown in Fig. 6 analogous to the main graph in Fig. 4.

But contrarily to the scaling of J40 and the other strongly coupled samples for which all data points followed the same power-law behavior irrespective of the applied field strength, the data for P40 clearly separate into a low field ($`H=0.1`$ and $`0.3`$ T, filled symbols) and a high field group ($`H0.6`$ T, open symbols). Assuming BG behavior at low fields, the $`IV`$-curves collapse very satisfactorily onto the two universal functions above and below the transition temperature and with comparable glass exponents as the BG phases in C40 and J40. The $`IV`$-curves at higher fields and temperatures above the estimated transition temperature collapse as well. Not surprisingly, the collapse at lower temperatures ($`T<T_g(H)`$) is not possible, because the $`IV`$-curves have still positive curvature. It is significant that a lower assumed transition temperature, placing the $`IV`$-curves at the lowest temperatures above the transition, leads to a violation of the collapse and scaling of $`R`$ and $`I_\text{x}`$. Thus in a C-oriented field of $`0.6`$ T the behavior above the transition temperature determined from scaling of $`R`$ and $`I_\text{x}`$ is exactly as expected for a liquid, but below that temperature the vortex fluid resembles neither a liquid nor a glass. Similar behavior is observed also for all investigated fields higher than 0.6 T.

Turning the magnetic field in the AC-direction so that field and defect structure are misaligned leads to the disappearance of even the low-field glass phase; at the lowest field $`H=0.1`$ T the $`IV`$-curves show only positive curvatures. Consequently, an attempt to analyze the data according to glass theories leads to the same difficulties as for the high field $`IV`$-curves and C-aligned fields.

The TAFF activation energies in Section IV suggest an interpretation for the lack of a clear second-order liquid-glass transition in this film. Recall that those energies had a field dependence characteristic of moving dislocations in two-dimensional superconductors for AC alignment, and that the C-aligned energies converged with those 2D values at fields above $`0.6`$ T. These observations suggest that the low-temperature vortex fluid state is essentially 2D for AC-aligned fields, while in C-aligned fields there is a crossover from a 3D glass phase at low fields to a mobile 2D phase at higher fields.

This conclusion relates directly to the qualitative change seen in the $`IV`$-characteristics. At low fields, the pancake vortices are highly correlated along the defect structure and at low temperatures they freeze into a Bose-glass phase. With increasing magnetic field the interlayer correlations are lost and the superconducting layers become effectively decoupled, preventing the formation of a BG phase as has been suggested also to occur in ion-irradiated Bi2212. Morozov et al. (1999) For misaligned magnetic fields (AC), on the other hand, the interlayer correlations between pancake vortices are not strong enough to form a 3D-VG for any applied field.

The true nature of the low-temperature, high-field phase remains unclear, however. For a truly 2D superconductor Eqs. 69 are not valid, instead it is expected that the glass phase exists at $`T=0`$ K, only. Fisher et al. (1991) At finite temperatures $`IV`$-curves should still be divided into a low-current linear and a high-current nonlinear part. The crossover current relation Eq. 9 is replaced by

$$I_\text{x}^{2D}T^{1+\nu _{2D}},$$
(10)
with $`\nu _{2D}=2`$. Analyzing the $`IV`$-characteristics of P40 according to Eq. 10 over the accessible temperature range leads to a non-constant exponent with $`\nu _{2D}30`$, most certainly ruling out a truly 2D-VG phase.

## VI Vortex Phase Diagrams

### VI.1 Strongly Coupled Layers

The analysis in the previous section puts us into a position to plot vortex phase diagrams for each of the three samples. Turning first to the strongly-coupled samples Fig. 7 shows the reduced $`H`$-$`T`$ phase diagram for J40.

The pattern is almost identical to that in C40. The uppermost line is the linearized upper critical field $`H_{c2}(T)`$ separating the normal from the superconducting phases. The vortex glass melting temperatures as derived from $`IV`$-curves are plotted as squares (C-) and circles (AC-alignment). At lower fields and temperatures the vortices form the respective Bose-glass and 3D-VG for C and AC-aligned fields. In the 3D-VG model the glass melting line can be described as a simple power-law Fisher et al. (1991)

$$[T_c(0)T_g]H^{1/2\nu _0}$$
(11)
with an expected $`\nu _02/3`$. The experimentally determined $`T_g(H)`$ do follow a power-law behavior very nicely, but with an exponent $`\nu _0=1`$.

The field dependence of the BG melting line for C-aligned magnetic fields is more complicated. Blatter *et al.*Blatter et al. (1992) give an expression derived for superconductors containing columnar defects and applied magnetic fields smaller than the matching field, but large enough that the vortices interact, $`B>\mathrm{\Phi }_0/\lambda ^2`$:

$$T_g(B)\chi T_m(B)+(1\chi )T_c\left(1\frac{B}{\mu _0H_{c2}(0)}\right),$$
(12)
with $`T_m(B)`$ the melting temperature of the flux line lattice (FLL) for the material in the clean limit without columnar or point defects. $`\chi `$ is a disorder parameter which equals unity in the clean limit and approaches zero for very strong disorder. The clean-limit melting temperature is unknown for these Ta$`_\text{x}`$Ge$`_{\mathrm{\hspace{0.17em}1}\text{x}}`$/Ge -films. Instead, we used the glass melting line for a homogeneous Ta0.27Ge0.73-film without extended defects but with weak-pinning point defects Ruck (1998) (dashed curve in Fig. 7) to approximate $`T_m`$ in the above equation. A fit of Eq. 12 then provides a satisfactory description of the field-dependence of the Bose-glass melting line and the disorder parameter $`\chi `$ a quantitative measure of the shift of the glass transition to higher temperatures.

This phenomenological fitting procedure yields the disorder parameters listed in Table 3, from which it can be seen that the disorder appears marginally stronger in film J40 than in C40. We have not been able to correlate this difference with any measured structural parameters, though it is worth noting that the superconducting layers are thinner in the more strongly disordered film.

### VI.2 Weakly Coupled Layers

The phase diagram for the weakly coupled film P40 is dominated by its quasi-2D behavior and a solid vortex glass phase is only realized for small, C-oriented fields. In Fig. 8 the reduced $`H`$-$`T`$ phase diagram is shown analogous to Fig. 7.

The uppermost solid line is the upper critical field $`H_{c2}(T)`$ and the lowest dotted curve is the glass melting line of the reference Ta0.27Ge0.73-film. Ruck (1998)

As outlined in Sec. V.2 Bose-glass scaling procedures could only be successfully applied to $`IV`$-curves taken at the two lowest fields applied parallel to the defect structure. The derived glass melting temperatures $`T_g`$ are given as the two lowest squares in Fig. 8. A glass melting line according to Eq. 12 has been successfully fitted to these data points (solid curve); the resulting disorder parameter indicates reduced disorder compared to the strongly coupled samples (see Table 3). At temperatures below $`T_g`$ the vortex system forms a solid BG (hatched area), which melts into a 3D vortex liquid at higher temperatures. The three dimensional character of the vortex liquid in this field range is supported by the field-dependence of the activation energies for TAFF (Sec. IV.2). Both 3D-phases are limited by an upper magnetic field, signaled in the case of the BG by the failing low temperature scaling and in the case of the vortex liquid deduced from activation energies. The approximate crossover to 2D behavior is given by the horizontal dotted line. It is worth mentioning that we have found signs of an reentrant 2D vortex liquid in Sec. IV.2 at very low magnetic fields, not shown in Fig. 8.

For higher magnetic fields and all field values applied in the AC-direction the $`IV`$-characteristics gave no evidence of a pinned solid vortex phase. The would-be transition temperatures sensed by the high-temperature vortex liquid are given as squares (C) and circles (AC). The dashed and dotted lines are guides to the eye and may mark the crossover to a highly viscous vortex liquid. This remains however an open question which demands further investigations.

A comparison of these phase diagrams with vortex phase diagrams for HTSCs is also interesting. The strongly coupled samples resemble those for weakly anisotropic Y123 \[Ref. (Krusin-Elbaum et al., 1994; Nakielski et al., 1996)\] and isotropic (K,Ba)BiO3. Klein et al. (2000) P40 on the other hand shows remarkable similarities with the moderately anisotropic Bi2Sr2Ca2Cu3O10 \[Ref. (Li et al., 1994)\] and especially Tl2Ba2Ca2Cu3O10. Yu and Yeh (1997)

## VII Dynamics in the Bose-glass Phase

The exponential $`IV`$-characteristics in the glass phase, Eq. 8, are caused by vortex creep due to finite applied currents and finite probabilities for vortex excitations at non-zero temperatures. These vortex dynamics in the glass phase can be characterized by the glass exponent $`\mu `$, with $`\mu >0`$ for a truly glassy phase. Various models for vortex creep in the presence of defects allow the calculation of the field and temperature-dependence of $`\mu `$. For weak pinning by point defects the weak collective pinning (WCP) theory offers one way to calculate the glass exponent $`\mu `$. However, $`\mu `$ depends non-trivially on the various length scales determining the size of the vortex bundles, which are the smallest independently moving objects in WCP theory. Feigel’man et al. (1989); Nattermann (1990); Blatter et al. (1994)

Within the Bose-glass theory for superconductors containing extended defects vortex creep is explained by the excitation of vortex loops, half-loops and kinks. Nelson and Vinokur (1993) Under the acting Lorentz-force due to the finite probing current, these excitations lead to the characteristic creep of the vortices through the pinning potential. These ideas can be developed for the case of pinning by columnar defects Nelson and Vinokur (1993) as well as for pinning by coplanar defects. Marchetti and Vinokur (1994, 1995) The problem can be mapped onto electron transport in doped semiconductors Mott (1969); Shklovskii and Efros (1984) wherefore it is also known as variable range hopping (VRH) of vortices. The following results are obtained for VRH by the creation of double superkinks

$$\mu =1/3\text{(columnar defects)},$$
(13)
$$\mu =1/2\text{(coplanar defects)}.$$
(14)
And for larger currents half-loop excitations can extend under the acting Lorentz-force until they reach the next favorable pinning site leading to a glass exponent

$$\mu =1.$$
(15)
Several groups have found some evidence for VRH in Y123 \[Ref. (Dekker et al., 1992; Thompson et al., 1997)\] and Bi2212. Soret et al. (1999, 2000); Ammor et al. (2000)

To check the applicability of these results to our films it is necessary to extract the exponent $`\mu `$ from $`IV`$-curves measured at temperatures below the glass melting temperature $`T_g`$. After recasting Eq. 8 into

$$\frac{\text{d}\left(\mathrm{ln}V/I\right)}{\text{d}I}I^{\left(\mu +1\right)}$$
(16)
this is most easily done by plotting the left side of Eq. 16 versus the current $`I`$ on a double-logarithmic scale. The data should then follow a straight line with gradient $`(\mu +1)`$ from which $`\mu `$ can be easily calculated. We apply this analysis in the following two sections.

### VII.1 Weakly Coupled Layers

We start the discussion of the experimental results with the weakly coupled sample P40. From the phase diagram Fig. 8 it can be seen that a solid glass phase exists only for small, C-oriented fields. In Fig. 9 the glass exponent $`\mu `$ is plotted as a function of the temperature below the glass melting temperature, $`T_gT`$, for $`H=0.1`$ T.

The inset shows the respective rescaled $`IV`$-curves and least square fits according to Eq. 16, demonstrating the applicability of above model to our data. Within the estimated errors the glass exponent $`\mu `$ is temperature independent with an averaged $`\overline{\mu }=0.22`$. This is approximately half the theoretically calculated value of $`1/2`$ for VRH in the presence of coplanar pinning sites (Eq. 14) and still significantly less than the other predicted value for pinning by columnar defects (Eq. 13).

Although above $`IV`$-characteristics can be described successfully with an exponential dependence according to Eq. 8, existing models describing the dynamics in the Bose-glass phase fail to explain the magnitude of the glass exponent. The available data did not allow for a conclusive determination of a field-dependence of $`\mu `$. The only other set of $`IV`$-curves taken in the Bose-glass phase at $`H=0.3`$ T is very close to the upper field limit of the glass phase. These $`IV`$-curves are very close to power-law behavior resulting in glass exponents $`\mu 0`$.

### VII.2 Strongly Coupled Layers

The glass exponent $`\mu `$ has a markedly different temperature-dependence in the strongly coupled samples C40 and J40, and the extended field-range of the glass phases additionally enabled us to study the field-dependence. In Fig. 10 $`\mu `$ is plotted versus $`T_gT`$ for a wide range of applied fields and temperatures below $`T_g(H)`$ in the case of C40 and C-oriented fields.

Very similar behavior was observed for magnetic fields applied in the AC-direction and the other strongly coupled film J40.

Looking at the field-dependence for fixed $`T_gT`$ the glass exponent $`\mu `$ decreases monotonically with increasing field. This is true for the complete observed temperature and field range; the latter spans more than one order in magnitude. The temperature-dependence shows even more interesting characteristics. As the temperature is lowered below the glass transition temperature $`T_g`$ the glass exponent increases down to temperatures approximately $`0.45`$ K below $`T_g`$. This increase may well be described by a linear dependence on temperature. Extrapolating $`\mu `$ to $`T_g`$ one obtains non-zero glass exponents at the glass transition, which decrease with increasing field as at lower temperatures. The rate at which $`\mu `$ increases also decreases with higher magnetic field. It is important to note that this happens in a temperature interval below $`T_g`$ for which we observed a constant $`\mu `$ in the weakly coupled sample P40 (compare Fig. 9).

At a temperature roughly $`0.45`$ K below $`T_g`$ (marked by the dashed vertical line in Fig. 10) a crossover to an approximately constant glass exponent is observed. 333For the two highest magnetic fields the temperature range was limited by our $`1.2`$ K limit. Therefore, the crossover could not be observed. The level is again field-dependent and ranges from about $`1.75`$ at low fields to estimated $`0.5`$ at the highest fields. The correlation between the crossover to a constant $`\mu `$ and $`T_gT`$ links the glass exponent to the formation of the glass phase, but is not explained by the simple VRH picture presented above.

## VIII Conclusions

We have performed extensive conductance measurements on artificially layered Ta$`_\text{x}`$Ge$`_{\mathrm{\hspace{0.17em}1}\text{x}}`$/Ge multilayers and compared the results for films with strong and weak interlayer coupling and with and without extended strongly pinning defects, respectively. Due to their large GL parameter $`\kappa `$ conduction below $`T_c(H)`$ is dominated by the dynamic behavior of vortices in response to an applied transport current and the interactions between vortices and the samples’ microstructure.

The strongly coupled films are best modelled as weakly anisotropic three dimensional superconductors. The activation energies for TAFF show a significant increase above the expected $`1/\sqrt{H}`$-dependence when the magnetic field is aligned with the extended pinning sites. The $`IV`$-characteristics were interpreted within the existing vortex glass models. A thorough analysis results in a very consistent second order phase transition from a vortex liquid at high temperatures to a pinned vortex glass phase at a field-dependent transition temperature $`T_g(H)`$. The critical glass exponents $`\nu `$ and $`z`$ are field-independent and all $`IV`$-curves taken at different fields and temperatures could be collapsed onto two scaling functions. The correlated defects lead to a large shift of the glass transition line to higher reduced temperatures when compared to an unstructured homogeneous film. This shift is maximized for fields well aligned with the correlated defect structure.

In the case of the more weakly coupled sample P40 one has to distinguish between fields parallel to the pinning sites (C-oriented) and magnetic fields at large angles to the columnar structure (AC-oriented). For the latter situation the film behaves like an anisotropic quasi-2D superconductor. The activation energies for TAFF are best described by a $`\mathrm{ln}H`$-dependence over the complete magnetic field range investigated, typical for 2D superconductors. Additionally, we have seen no real signs of a second order phase transition into a glass phase at low temperatures, in line with predictions that a glass phase in two dimensions should at best exist at $`T=0`$. On the other hand, turning the magnetic field into the C-direction leads to a reoccurrence of 3D vortex behavior in a narrow field range. The strongest argument for this interpretation is the formation of a BG phase at low magnetic fields. This three dimensional vortex phase is supported by the extended defects which assist the formation of well-aligned stacks of pancake vortices leading to increased correlations across the superconducting layers as suggested by results on Bi2212. Sato et al. (1997); Kosugi et al. (1997); Morozov et al. (1999) This picture is further supported by the qualitative change observed in the $`IV`$-characteristics between high and low fields.

Although the low-temperature $`IV`$-curves can be described very consistently within the BG model there are significant differences in the vortex creep between the strongly and weakly coupled samples and neither of them is fully explained by the simple VRH mechanism. This is particularly true for vortex creep in the strongly coupled samples. Since the glass exponent $`\mu `$ showed very similar $`T`$ and $`H`$-dependence for C and AC-aligned fields, details of the pinning sites, such as point-like or extended, seem to be less important and the magnitude of $`\mu `$ varied over a wide range with temperature and field. However, the strong correlation between the crossover from increasing $`\mu `$ to a constant value and the temperature difference $`T_gT`$ suggests a connection between the formation of the glass phase and the dynamics of the vortices. More specific, two fundamental variables of the glass theories, the correlation length $`\xi `$ and the relaxation time $`\tau `$, are proportional to $`|T_gT|^x`$, where $`x`$ is determined by the critical exponents $`\nu `$ and $`z`$. To elucidate a possible connection between the $`T`$ and $`H`$-dependence of $`\mu `$ on the one hand and $`\nu `$ and $`z`$ on the other, more theoretical work in this area is called upon.

The authors acknowledge support by the New Zealand Marsden Fund and thank J. C. Abele for assistance in gathering experimental data.


format: text

text:
# AN ALGORITHM GENERATING SCALE FREE GRAPHS

## 1 Introduction

A wide variety of systems in biology, communication technology, sociology and economics are best described as complex networks of agents linked by various physical and informational links. Clearly, the architecture of such networks is crucially important for the vitality and performance of these systems. It is generally recognized that a thorough investigation of possible mechanisms which determine the topology of complex ”real-world” networks is necessary for understanding their behavior. Numerous studies of the Internet - and World Wide Web (WWW) -, the web of human sexual contacts  and movie actors collaboration networks -, national power-grids  and phone-call networks -, protein folding  and cellular networks -, etc. have revealed, among other facts, that the probability that a node of these networks has $`k`$ connections follows a power law

$$P\left(k\right)k^\gamma $$
(1)
over a large range of $`k`$, with an exponent $`\gamma `$ that ranges between 1 and 3 depending on the system . Such networks are called scale-free . Many other ”real-world” networks display an inherently uniform exponential topology that still deviates from a Poisson distribution expected for usual random graphs $`𝔾(N,p)`$.

Recently, a model generating a scale-free network has been suggested in . It establishes a random graph process in which vertices are added to the graph one at a time and connected to a fixed number of earlier vertices, selected with probabilities proportional to their degrees. This preferential attachment assumption is based on the idea that a new site is more likely to link to existing sites which are ”popular” at time the site is added. Numerical simulations for the above process reported in  show that after many iterations the proportion of vertices with degree $`k`$ indeed obeys the power law (1) with the exponent $`\gamma =3.`$ In , the power law statistics for $`P\left(k\right)`$ in this model has been proven analytically, and it has been justified that $`\gamma =3`$.

Probably, a complex network formed by the citation patterns of scientific publications where the nodes standing for published articles and a directed edge representing a reference to a previously published article exhibiting the power law with the exponent $`\gamma =3`$ , can be well fitted by the above process as discussed in . However, in other ”real-world” networks having degree sequences matching a power law these exponents usually deviate from $`3`$. For example, $`\gamma =2.3\pm 0.1`$ for the movie actors collaboration network , $`\gamma =2.1`$ for a graph of long distance telephone calls made during a single day , etc. To fit these deviations, one assumes that the preferential attachment for a real network is nonlinear , i.e. the probability that a node $`i`$ attaches to node $`j`$ can be characterized by some function $`\mathrm{\Pi }\left(k\right)`$ also following a power law with some index $`\alpha .`$

The effect of nonlinear preferential attachment $`\mathrm{\Pi }\left(k\right)`$ on the network dynamics and topology was studied in . They demonstrated that the scale free network topology holds only if $`\mathrm{\Pi }\left(k\right)`$ is asymptotically linear and is destroyed otherwise. This way the exponent of the degree distribution $`P\left(k\right)`$ can be tuned to any value between 2 and $`\mathrm{}`$ .

Nevertheless, the detailed studies of various ”real-world” networks still challenge to the preferential attachment approach . Some ecological networks such as food webs quantifying the interactions between species , some networks in linguistics , a network representing an aggregation of the WWW domains , etc. exhibit a scale free topology characterized by the indices $`1<\gamma <2`$. Moreover, for some networks, the in-degree and out-degree probability distributions are dissimilar: the WWW-topology where the nodes are the webpages and the edges are the hyperlinks that point from one document to another meets the power laws with $`\gamma _{\mathrm{in}}=2.1`$ and $`\gamma _{\mathrm{out}}=2.72`$ , the distributions of the total number of sexual partners in a single year indicate the power law decay with $`\gamma =2.54\pm 0.2`$ for females and $`\gamma =2.31\pm 0.2`$ for males . A recent study  extends the results on the citation patterns (incoming links) of scientific publications ($`\gamma _{\mathrm{in}}=3`$, ) to the outgoing degree distribution obtaining that it has an exponential tail (not a Poisson distribution). This challenge stimulates a search for the new models inducing graphs of the scale free topology.

A random procedure which we use to generate scale free graphs draws back to the ”toy” model for a system being at a threshold of stability reported in  and to various coherent-noise models - discussed in connection to the standard sandpile model  developed in self-organized criticality, where the statistics of avalanche sizes and durations also take a power law form. The random process in our model begins with no edges at time $`0,`$ at a chosen vertex $`i,`$ and adds new edges linked to $`i`$, one at a time; each new edge is selected at random, uniformly among all outgoing edges not already present. When at some time a performance of the system occasionally exceeds a threshold value, the system looses its stability, so that the process stops at the given vertex and then starts again at some other vertex lacks of outgoing edges. This way the out-degree distribution of the resulting graph equals to the distribution of residence times below the threshold which is closely related to the ergodic properties of the system and its complexity, .

Although it has no a definite relation to the preferential attachment principle , probably, both algorithms are similar: we also monitor a global state of the system (while inspecting its stability) and introduce an auxiliary function (which is, in some sense, analogous to $`\mathrm{\Pi }\left(k\right)`$) to characterize the ergodic properties of the system. Nevertheless, the model we propose in the present paper has also a distinguishing feature: we allow for the threshold value to quench with a given probability $`\eta `$, so that a natural control parameter in our model is the relative frequency of threshold changes. Varying this frequency, one can control the statistics of switching events (when the system performance parameter passes the stability threshold) tuning it from the truly exponential decay distribution (when the threshold representing the environmental stress changes each time a new edge added to the graph) to a power law with some exponent $`\gamma `$ when the threshold is frozen.

Among clear advantages of the proposed algorithm, one can mention its flexibility (the possible value of index $`\gamma `$ varies in the wide range $`(1,\mathrm{}),`$ the in-degree and out-degree statistics can be tuned independently, the variation of frequency $`\eta `$ switches $`P\left(k\right)`$ from an exponential decay to a power law) as well as its relative simplicity for both analytical and numerical studies.

The paper is organized as follows. In Sec. 2, we describe the model we introduce. In Sec. 3, we investigate the degree statistics for the random graphs generated in accordance with the introduced random procedure. Section 6 is devoted to conclusions.

## 2 Motivation and description of the model

Probably, the most fascinating property of living systems is a surprising degree of tolerance against environmental interventions revealing the robustness of the underlying metabolic and genetic networks . A typical question which is usually posed in the literature in concern to this feature is on the role of the network topology in the error tolerance (see  and references therein). However, in the present work, we bring an opposite question to a focus: what is the network topology (say, of the protein-protein interactions map) resulted from the very durable evolution process selecting from among species enjoying sufficiently good adaptivity? A principle of evolutionary selection of a common large-scale structure of biological networks put forward in  is closely related to this question.

The most of biological systems exhibiting a high degree of adaptivity can be referred as systems being close to a stability threshold, so that a slight modification of the environment would prompt a drastic change in the system performance increasing or reducing the stamina of biological spices. From the dynamical systems point of view, such a behavior can be interpreted as a forcing of bifurcation parameter through a bifurcation point triggering the system into a new state. In our study, we do not refer to any definite physical system displaying such a complex behavior. Instead, following , we use a ”toy” model, in which the bifurcation parameter (the current state of the system characterized by its net performance) and the bifurcation point (the threshold of stability characterizing a challenge of the environments) are considered as random independent variables. It is supposed that the system is subjected to a bifurcation, when its net performance takes on values above the threshold. It is also important that all network’s changes acquired in the previous stages of evolution of the graph are not discarded, when the system suffers a bifurcation, but still into effect for the forthcoming evolution.

We shall characterize the current performance of the whole system by the real number $`x[0,1].`$ Another real number $`y[0,1]`$ plays the role of the stability threshold. The network is supposed to be stable until $`x<y`$ and is condemned otherwise. We consider $`x`$ as a random variable distributed with respect to some given probability distribution function (pdf) $`f\left(u\right):[0,1]^+`$. In a deterministic version of the model, $`f`$ would be an ergodic invariant density function on the phase space which has an essential influence on the system dynamics close to the stability threshold value.

The value of threshold $`y`$ also fluctuates reflecting the changes of environments and is also treated as a random variable distributed over the interval $`[0,1]`$ with respect to some other pdf $`g\left(u\right):[0,1]^+`$. In general, $`f`$ and $`g`$ are two arbitrary positive integrable functions satisfying the normalization condition

$$_0^1𝑑uf\left(u\right)=_0^1𝑑ug\left(u\right)=1.$$
(2)
We define the random procedure on the fixed set of $`N`$ vertices describing the evolution of network as a random graph evolving in time. The process begins with no edges at time $`0,`$ at a chosen vertex $`i`$. Given a fixed real number $`\eta [0,1]`$, we define a discrete time random process in the following way. At time $`t=0,`$ the variable $`x`$ is chosen with respect to pdf $`f`$, and $`y`$ is chosen with pdf $`g`$. If $`x<y`$, we add an edge selected at random, uniformly among all edges outgoing from $`i`$, the process continues and goes to time $`t=1`$. At time $`t1,`$ one of the following events happens:

i) with probability $`\eta `$, the random variable $`x`$ is chosen with pdf $`f`$ but the threshold $`y`$ keeps the value it had at time $`t1`$. Otherwise,

ii) with probability $`1\eta ,`$ the random variable $`x`$ is chosen with pdf $`f,`$ and the threshold $`y`$ is chosen with pdf $`g`$.

If $`xy,`$ the process at the given vertex ends and then resumes at some other vertex chosen randomly from the vertices having no outgoing edges yet; if $`x<y,`$ the process continues at the same vertex and goes to time $`t+1.`$

Eventually, when at some time $`T`$ limiting the duration of the ”monotonous” evolution of the graph, the process moves to some other vertex, and $`i`$ has exactly $`k=T`$ outgoing edges.

For example, in Fig. 1 we have presented the results of numerical simulation for the above algorithm for a network consisting of $`2\times 10^3`$ vertices, for $`f=g=1`$ and $`\eta =1`$. It displays the proportion of vertices having degree $`k`$ versus $`k`$ in the log-log scale. Circles stay for the out-degree $`k_{\mathrm{out}}`$ and diamonds are for the in-degree $`k_{\mathrm{in}}.`$ One can see that the out-degree probability distribution is fitted well by a power law, while the in-degree probability distribution forms a Gaussian profile centered at the most probable number of incoming links $`k_{\mathrm{in}}^{}=11.`$

The same algorithm can be extended readily to moderate the incoming edges also. Let us consider three random variables $`x,y,`$ and $`z`$ that are the real numbers distributed in accordance to the distributions $`f,`$ $`g`$, and $`v`$ within the unit interval $`[0,1].`$ We assume that $`x`$ represents the current performance of the system, $`y`$ and $`z`$ are the thresholds for emission and absorption of edges respectively.

The generalized process begins on the set of $`N`$ vertices with no edges at time $`0,`$ at a chosen vertex $`i`$. Given two fixed numbers $`\eta [0,1]`$ and $`\nu [0,1]`$, the variable $`x`$ is chosen with respect to pdf $`f`$, $`y`$ is chosen with pdf $`g`$, and $`z`$ is chosen with pdf $`v`$, we draw $`e_{ij}`$ edge outgoing from $`i`$ vertex and entering $`j`$ vertex if $`x<y`$ and $`x<z`$ and continue the process to time $`t=1.`$ Otherwise, if $`xy`$ ($`xz`$), the process moves to other vertices having no outgoing (incoming) links yet.

At time $`t1,`$ one of the three events happens:

i) with probability $`\eta `$, the random variable $`x`$ is chosen with pdf $`f`$ but the thresholds $`y`$ and $`z`$ keep their values they had at time $`t1`$.

ii) with probability $`1\eta ,`$ the random variable $`x`$ is chosen with pdf $`f,`$ and the thresholds $`y`$ and $`z`$ are chosen with pdf $`g`$ and $`v`$ respectively.

iii) with probability $`\nu ,`$ the random variable $`x`$ is chosen with pdf $`f`$, and the threshold $`z`$ is chosen with pdf $`v`$ but the threshold $`y`$ keeps the value it had at time $`t1`$.

If $`xy`$, the process stops at $`i`$ vertex and then starts at some other vertex having no outgoing edges yet. If $`xz,`$ the accepting vertex $`j`$ is blocked and does not admit any more incoming link (provided it has any). If $`x<y`$ and $`x<z`$, the process continues at the same vertex $`i`$ and goes to time $`t+1.`$

Figure 2 displays an example of a graph generated in accordance this algorithm on the set of 100 vertices for the case of uniform densities $`f=g=v=1`$ for $`\eta =\nu =1.`$ In Fig. 3, we have presented the distribution $`P\left(k\right)`$ vs. $`k`$ in the log-log scale over $`10^5`$ vertices. Here, the circles stay for outgoing degrees, and diamonds are for incoming degrees. For $`k1,`$ both profiles enjoy a power law decay with $`\gamma _{\mathrm{in}}=\gamma _{\mathrm{out}}=2.`$

In the forthcoming section, we demonstrate that for some probability densities $`f,g`$, and $`v`$ the resulting probability distribution $`P\left(k\right)`$ meets asymptotically a power law (1) as $`k1`$ where the index $`\gamma `$ can vary in the wide range $`\gamma (1,\mathrm{})`$. Moreover, by means of adjusting probability densities $`f,g,`$ and $`v`$ from a proper class of functions, one can tune $`\gamma _{\mathrm{in}}`$ and $`\gamma _{\mathrm{out}}`$ to different values.

If $`\nu =0`$, both thresholds $`y`$ and $`z`$ have synchronized dynamics, and sliding the value of $`\eta `$ form 0 to 1, one can tune the statistics of out-degrees and in-degrees simultaneously out from the pure exponential decay to the power laws provided $`f,`$ $`g`$, and $`v`$ belong to the proper class of functions. When $`\nu =1`$, the threshold $`z`$ changes each time but the threshold $`y`$ is frozen. In this case, $`P\left(k_{\mathrm{in}}\right)`$ decays exponentially for any choice of $`v`$, but $`P\left(k_{\mathrm{out}}\right)`$ still meets a power law for $`k_{\mathrm{out}}1`$.

## 3 Probability degree distribution $`P(k)`$

Below, we compute the probability $`P\left(k\right)`$ that a given vertex has precisely $`k`$ outgoing edges as a result of the random procedure defined in the previous section. For sake of clarity, we investigate the statistics of switching events quenching the monotonous evolution of the graph in the simplest case when only one threshold ($`y`$) is present. All computations can be extended readily to the general case when both thresholds are present.

Directly from the definitions of Sec. 2, one can show that

$$P\left(0\right)=_0^1𝑑yg\left(y\right)_y^1𝑑xf\left(x\right).$$

In $`t1`$ time steps, the system can either ”survive” (S) or ”die” (D). Both issues can take place either in the ”correlated” way (with probability $`\eta `$, see i) (we denote these scenarios as $`S_c`$ and $`D_c`$ consequently), or in the ”uncorrelated” way (with probability $`1\eta `$, see ii) ($`S_u`$ and $`D_u`$, consequently). For example, for $`k=1`$, we have

$$P\left(1\right)=P\left[SD_c\right]+P\left[SD_u\right]=\eta B\left(1\right)+\left(1\eta \right)A^2\left(1\right)$$

where we have defined, for $`n=0,1,2,\mathrm{},`$

$$A\left(n\right)=_0^1𝑑yg\left(y\right)\left[_0^y𝑑xf\left(x\right)\right]^n$$
(3)
and

$$B\left(n\right)=_0^1𝑑yg\left(y\right)\left[_0^y𝑑xf\left(x\right)\right]^n_y^1𝑑uf\left(u\right)=A\left(n\right)A\left(n+1\right).$$
(4)
Similarly,

$$P\left(2\right)=P\left[SS_cD_c\right]+P\left[SS_cD_u\right]+P\left[SS_uD_c\right]+P\left[SS_uD_u\right],$$

where

$$P\left[SS_cD_c\right]=\eta ^2B\left(2\right),P\left[SS_cD_u\right]=\eta \left(1\eta \right)A\left(2\right)B\left(0\right),$$

$$P\left[SS_uD_c\right]=\eta \left(1\eta \right)A\left(1\right)B\left(1\right),P\left[SS_uD_u\right]=\left(1\eta \right)^2A^2\left(1\right)B\left(0\right).$$

The general formula for $`k3`$ is tedious and can be found in .

Computations can be essentially simplified with the use of generating functions. Let us define the generating function for $`P\left(k\right)`$ by

$$\widehat{P}\left(s\right)=\underset{k=0}{\overset{\mathrm{}}{}}s^kP\left(k\right),P\left(k\right)=\frac{1}{k!}\frac{d^k\widehat{P}\left(s\right)}{ds^k}|_{s=0},$$
(5)
and introduce the new auxiliary functions,

$$p\left(l\right)=\eta ^lA\left(l+1\right),\mathrm{for}l1,p\left(0\right)=0,$$

$$q\left(l\right)=\left(1\eta \right)^lA^{l1}\left(1\right),\mathrm{for}l1,q\left(0\right)=0,$$

$$r\left(l\right)=\eta ^l\left[\eta B\left(l+1\right)+\left(1\eta \right)A\left(l+1\right)B\left(0\right)\right],\mathrm{for}l1,r\left(0\right)=0,$$

$$\rho =\eta B\left(1\right)+\left(1\eta \right)A\left(1\right)B\left(0\right).$$
(6)
Then, the use of convolutions property  yields

$$\widehat{P}\left(s\right)=B\left(0\right)+\rho s+\frac{s}{1\widehat{p}\left(s\right)\widehat{q}\left(s\right)}\left[\widehat{r}\left(s\right)+\rho \widehat{p}\left(s\right)\widehat{q}\left(s\right)+\rho A\left(1\right)\widehat{q}\left(s\right)+A\left(1\right)\widehat{q}\left(s\right)\widehat{r}\left(s\right)\right],$$
(7)
where $`\widehat{p}\left(s\right),\widehat{q}\left(s\right),\widehat{r}\left(s\right)`$ are the generating functions of $`p\left(l\right),q\left(l\right),r\left(l\right),`$ respectively.

In the marginal cases of $`\eta =0`$ and $`\eta =1`$, the probability $`P\left(k\right)`$ can be readily calculated.

### 3.1 The uncorrelated case, $`\eta =0`$

This case corresponds to the external stress being sensitive to any modification of a network and changing its figure as soon as a new edge appears in the graph. For $`\eta =0,`$ one has $`\widehat{p}\left(s\right)=\widehat{r}\left(s\right)=0,`$ and $`\widehat{q}\left(s\right)=s\left(1sA\left(1\right)\right)^1,`$ $`\rho =A\left(1\right)B\left(0\right)`$, then equation (7) gives

$$\widehat{P}_{\eta =0}\left(s\right)=\frac{B\left(0\right)}{1sA\left(1\right)}.$$
(8)
Applying the inverse formula in (5) to equation (8), one arrives at

$$P_{\eta =0}\left(k\right)=A^k\left(1\right)B\left(0\right)=\left[_0^1𝑑yg\left(y\right)_0^y𝑑xf\left(x\right)\right]^k_0^1𝑑yg\left(y\right)_y^1𝑑xf\left(x\right).$$

Therefore, in the ”uncorrelated” case ($`\eta =0`$), for any choice of the pdf $`f\left(u\right)`$ and $`g\left(u\right),`$ the probability $`P\left(k\right)`$ decays exponentially.

### 3.2 The correlated case, $`\eta =1`$

Another marginal case corresponds to the external stress remaining unchanged during all time of network evolution. For $`\eta =1`$, one has $`\widehat{q}\left(s\right)=0,`$ $`\rho =B\left(1\right),`$ and

$$s\widehat{r}\left(s\right)=\underset{n=1}{\overset{\mathrm{}}{}}s^{n+1}B\left(n+1\right)=\underset{n=2}{\overset{\mathrm{}}{}}s^nB\left(n\right),$$

so that formula (7) yields

$$\widehat{P}_{\eta =1}\left(s\right)=\widehat{B}\left(s\right).$$
(9)
Consequently,

$$P_{\eta =1}\left(k\right)=B\left(k\right)=_0^1𝑑yg\left(y\right)\left[_0^y𝑑xf\left(x\right)\right]^k_y^1𝑑uf\left(u\right).$$
(10)
In the ”correlated” case, $`\eta =1`$, many different types of behaviour are possible, depending on the form of the pdf $`f\left(u\right)`$ and $`g\left(u\right)`$. We will examine an important class of $`f\left(u\right)`$ and $`g\left(u\right)`$, for which $`P_{\eta =1}\left(k\right)`$ can be explicitly computed from equation (10). We will take

$$\begin{array}{cc}f\left(u\right)=\left(1+\alpha \right)u^\alpha ,\hfill & \alpha >1,\\ g\left(u\right)=\left(1+\beta \right)\left(1u\right)^\beta ,\hfill & \beta >1.\end{array}$$
(11)
Equation (10) gives in this case:

$$P_{\eta =1}\left(k\right)=\frac{\mathrm{\Gamma }\left(2+\beta \right)\mathrm{\Gamma }\left(1+k\left(1+\alpha \right)\right)}{\mathrm{\Gamma }\left(2+\beta +k\left(1+\alpha \right)\right)}\frac{\mathrm{\Gamma }\left(2+\beta \right)\mathrm{\Gamma }\left(1+\left(k+1\right)\left(1+\alpha \right)\right)}{\mathrm{\Gamma }\left(2+\beta +\left(k+1\right)\left(1+\alpha \right)\right)}.$$

Using Stirling’s approximation, we get for large $`k>>1`$:

$$P_{\eta =1}\left(k\right)=\frac{\left(1+\beta \right)\mathrm{\Gamma }\left(2+\beta \right)\left(1+\alpha \right)^{1\beta }}{k^{2+\beta }}\left(1+0\left(\frac{1}{k}\right)\right).$$

For different values of $`\beta `$, the exponent of the threshold distribution, we get all possible power law decays of $`P_{\eta =1}\left(k\right)`$. Notice that the exponent $`\left(2\beta \right)`$ characterizing the decay of $`P_{\eta =1}\left(k\right)`$ is independent of the distribution $`f\left(u\right)`$ of the state variable.

Therefore, for $`\eta =1`$, one has to provide just one power law function (which defines the statistics of the threshold) to achieve a power law for $`P\left(k\right)`$ characterized by the exponent $`\gamma (1,\mathrm{}).`$

### 3.3 The case of uniform densities $`f=g=1`$

For the case of uniform densities $`f=g=1,`$ and after some tedious but trivial computation, we get from equation (7):

$$\widehat{P}\left(s\right)=\frac{1}{1+\left(1\eta \right)\mu \left(s\right)}\left[\frac{1+\mu \left(s\right)}{s}\eta \mu \left(s\right)\right],$$
(12)
where $`\mu \left(s\right)`$ is defined by

$$\mu \left(s\right)=\frac{\mathrm{ln}\left(1\eta s\right)}{\eta s}.$$

The asymptotic behavior of $`P\left(k\right)`$ for $`k`$ large is determined by the singularity of the generating function $`\widehat{P}\left(s\right)`$ that is closest to the origin.

For $`\eta =0,`$ the generating function has a simple pole $`\widehat{P}\left(s\right)=\left(2s\right)^1,`$ and therefore $`P\left(k\right)`$ decays exponentially, which agrees with the general conclusion of Sec. 3.1.

For intermediate values $`1>\eta >0`$, the generating function $`\widehat{P}\left(s\right)`$ has two singularities. One pole, $`s=s_0,`$ corresponds to the vanishing denominator $`1+\left(1\eta \right)\mu \left(s\right),`$ where $`s_0`$ is the unique nontrivial solution of the equation

$$\mathrm{ln}\left(1\eta s\right)=s\frac{\eta }{1\eta }.$$
(13)
Another singularity, $`s=s_1,`$ corresponds to the vanishing argument of the logarithm, $`s_1=\eta ^1,`$ such that $`s_1>s_0>1.`$ The dominant singularity of $`\widehat{P}\left(s\right)`$ is of polar type, and the corresponding decay of $`P\left(k\right)`$ is asymptotically exponential for the degrees $`k`$ larger than $`k_c\mathrm{ln}\left(s_0\right)^1`$, with a rate $`\mathrm{ln}\left(s_0\right)`$ that vanishes like $`1\eta `$ as $`\eta `$ tends to 1.

When $`\eta `$ tends to one, two singularities $`s_0`$ and $`s_1`$ merge. More precisely, we have

$$\widehat{P}_{\eta =1}\left(s\right)=\frac{s+\left(1s\right)\mathrm{ln}\left(1s\right)}{s^2}.$$
(14)
The corresponding dominant term in (14) is of order $`\mathrm{O}\left(k^2\right)`$ . This obviously agrees with the exact result one can get from equation (10), with $`f\left(x\right)=g\left(x\right)=1`$:

$$P_{\eta =1}\left(k\right)=\frac{1}{\left(k+1\right)\left(k+2\right)}.$$
(15)
In the case of uniform densities, it is possible to get an expression of $`P\left(k\right)`$ for all $`k`$, and for any value of $`\eta `$, by applying the inversion formula (5) to equation (14) and use of the formulae . Combining all this, we get

$$P\left(k\right)=\frac{\eta ^k}{\left(k+1\right)\left(k+2\right)}+\underset{l=1}{\overset{k}{}}\frac{\eta ^k}{\left(kl+1\right)\left(kl+2\right)l}\underset{m=1}{\overset{l}{}}\left(\frac{1\eta }{\eta }\right)^mc(m,l),$$
(16)
where $`c(m,l)`$ is defined by

$`c(m,l)=m!{\displaystyle \underset{\begin{array}{c}l_1+l_2+\mathrm{}+l_m=l\\ l_i1\end{array}}{}}`$
$`{\displaystyle \frac{l_1l_2\mathrm{}l_{m1}l_m}{\left(l_1+1\right)\left(ll_1\right)\left(l_2+1\right)\left(ll_1l_2\right)\mathrm{}\left(l_{m1}+1\right)\left(ll_1l_2\mathrm{}l_{m1}\right)\left(l_m+1\right)}}.`$
When $`\eta 0`$, there is an alternative way of writing the previous expression:

$$P\left(k\right)=\frac{\eta ^k}{\left(k+1\right)\left(k+2\right)}+\underset{l=1}{\overset{k}{}}\frac{\eta ^{k+1}}{\left(kl+1\right)\left(kl+2\right)}\underset{m=1}{\overset{\mathrm{}}{}}\left(1\eta \right)^mb(m,l),$$
(18)
where $`b(m,l)`$ is defined by

$$b(m,l)=\underset{\begin{array}{c}i_1+i_2+\mathrm{}+i_l=l\\ i_j0\end{array}}{}\frac{1}{\left(i_1+1\right)\left(i_2+1\right)\mathrm{}\left(i_{m1}+1\right)\left(i_m+1\right)}.$$

In Fig. 4, we have plotted these asymptotes in the log-linear scale for the consequent values $`\eta =0,`$ $`\eta =0.5`$, $`\eta =0.7`$, $`\eta =0.9`$, $`\eta =1`$ (bottom to top) together with the curves representing analytical results (15-16).

Notice that for $`\eta 1`$, we only plot distributions $`P\left(k\right)`$ up to relatively small degrees ($`k=16`$, $`k=25`$) correspondent to short quiescent times in the system, since these times are already bigger then the crossover value $`k_c\left(\eta \right)1/\mathrm{ln}\left(s_0\right)`$ to the exponential decay $`\mathrm{exp}\left(\mathrm{ln}\left(s_0\right)k\right)`$ (where $`s_0`$ is defined by the Eq. (13)). For much longer times, very few survivals are observed, and the statistics gets bad. Of course, $`k_c\left(\eta \right)`$ grows as the parameter $`\eta `$ tends to 1, so that we have good statistics for larger $`k`$.

## 4 Discussion and Conclusion

In this paper, we have presented a simple random procedure inducing various types of random graphs and the scale free random graphs among others. The model we use essentially differs from the preferential attachment approach  discussed in the literature before. The proposed approach is simple and more flexible both for numerical simulations and analytical studies. The model draws back to the original ”toy” model for a system being at a threshold of stability .

In the proposed model, the degree statistics of a random graph depends from the value of the control parameter $`\eta `$ stirring the pure exponential statistics of degrees (at $`\eta =0,`$ when the threshold value is sensitive to each new edge risen in the graph) to a power law (at $`\eta =1`$, when the threshold value is frozen). The exponent characterizing the power law depends as $`\gamma =2+\beta `$ from the power law probability density distribution of the threshold $`g\left(u\right)\left(1u\right)^\beta ,`$ $`\beta >1,`$ and can vary in the wide range $`\gamma (1,\mathrm{}).`$ For the intermediate values of $`\eta `$, the decay rate is mixed. In particular, we showed that, even if the asymptotic behavior is exponential, when $`\eta `$ is close to 1 this only holds for very large degrees $`k`$ of order of $`\left(1\eta \right)^1`$. The model can be readily generalized to various particular types of networks. The crucial advantage of the model is that, for most physically meaningful distributions, it can be solved analytically. In other cases, a numerical solution can be readily found.

In the present paper, we have limited the consideration to the special cases of uniform and power law distributions of thresholds. Indeed, one can take other (normalized) distributions for $`f,g,`$ and $`v`$ in this model and obtain different decaying asymptotic profiles for $`P\left(k\right)`$ either for synchronous $`\nu =0`$ or asynchronous $`\nu =1`$ dynamics of thresholds, they have, in general, nothing to do with power laws at $`\eta =1`$, but still uniformly exponential at $`\eta =0.`$

## 5 Acknowledgements

One of the authors (D.V.) benefits from a scholarship of the Alexander von Humboldt Foundation (Germany) that he gratefully acknowledges.