text:
# On the origin of the pseudogap in underdoped cuprates

## Introduction

The existence of a pseudogap, i.e. a large suppression of low-frequency spectral weight in the normal state of underdoped high-temperature superconductors, is now a commonly accepted experimental fact. A multitude of experiments has probed the magnetic, thermodynamic, transport and optical properties of the underdoped cuprates (for an overview see Refs. timusk ; tallon ).

Some of the earliest indications of a pseudogap in the spin-channel were found in NMR-experiments warren , where the spin-susceptibility as measured by the Knight-shift was seen to decrease with decreasing temperature well above the superconducting critical temperature $`T_c`$. Further indirect evidence of a pseudogap was found in specific heat measurements loram . In underdoped samples, the electronic contribution starts to decrease with decreasing temperature in the normal state well above $`T_c`$. In transport measurements the crossover of the linear dependence of the ab-plane resistivity in temperature to a stronger dependence takagi  was taken as evidence for the opening of a pseudogap. However, the most direct and reliable measurements of the normal state pseudogap are angle-resolved photoemission experiments ding . A highly anisotropic ($`d_{x^2y^2}`$-like) suppression of low-energy spectral weight persisting far above $`T_c`$ has been found in underdoped samples. This pseudogap closes at a temperature $`T^{}`$ consistent with the crossover temperatures determined from other measurements.

Theorists have responded with several scenarioskampf ; emery ; anderson ; varma ; anderson2 ; zhang ; rice  for the pseudogap. Scenarios, based on the vicinity of the underdoped system to antiferromagnetic ordering, hold the coupling of quasiparticles to antiferromagnetic spin fluctuations responsible for the pseudogap behavior kampf : As a consequence of short-ranged antiferromagnetic correlations, a shadow-band forms and a pseudogap opens up in the density of states.

Motivated by the fact that the pseudogap has the same symmetry as the superconducting gap, ”preformed pair” scenarios associate the pseudogap with fluctuations that lead to $`d`$-wave superconductivity emery . In these scenarios the transition to the normal state is controlled by the vanishing of the superfluid density and not the closing of the superconducting gap. Due to strong pairing correlations above $`T_c`$, precursors of the superconducting gap are seen as pseudogap in the low-energy excitations.

Ideas involving spin-charge separation are based on the resonating valence bond (RVB) picture anderson . In these scenarios the pseudogap is due to d-wave singlet pairing of spin $`1/2`$, charge neutral fermions, called spinons. At $`T_c`$ the holons, i.e. spin 0 charge excitations, become coherent and recombine with spinons to form electron pairs which renders the system superconducting.

Further scenarios are based on the existence of a quantum critical point (QCP) close to optimal doping that associate the pseudogap with a broken symmetry (see eg. Ref. varma ). Other approaches that ascribe the formation of the pseudogap to the presence of strong charge fluctuations are motivated by the observation of charge stripes in some cuprate superconductors.

Presently, there is no consensus as to the origin of the pseudogap in the underdoped cuprates so numerical investigations based on model calculations are highly desirable. Early in the history of high-$`T_c`$ superconductors it was realized that the two-dimensional (2D) Hubbard model in the intermediate coupling regime, i.e. where the Coulomb interaction is of the order of the bandwidth, or closely related models like the t-J-model, should capture the essential low-energy physics of the cuprates anderson2 .

In this report we study the pseudogap behavior of the 2D Hubbard model using the dynamical cluster approximation DCA\_Hettler1 ; DCA\_Hettler2 ; DCA\_Maier1 ; DCA\_Jarrell1  (DCA). The DCA is a non-perturbative approach for the thermodynamic limit, which systematically incorporates the effects of nonlocal correlations to local approximations like the dynamical mean field approach (DMFA) DMFA , by mapping the lattice onto a self-consistently embedded cluster of size $`N_c`$. We solve the cluster problem using a combination of Quantum Monte Carlo (QMC) and the maximum entropy method to obtain dynamics DCA\_Jarrell3 . Since further insight into the pseudogap phenomenon can be obtained by studying the effects of lattice frustration, we also present results for the 2D Hubbard model with additional next nearest neighbor hopping $`t^{}`$. We show that our results for the pseudogap behavior in the 2D Hubbard model support the RVB picture and exclude other mechanisms for the pseudogap.

## Formalism

A detailed discussion of the DCA formalism was given in previous publications DCA\_Hettler1 ; DCA\_Hettler2 ; DCA\_Maier1 ; DCA\_Jarrell1 ; DCA\_Jarrell3 . The DCA is based on the assumption thurston  that the lattice self-energy is only weakly momentum dependent and can be approximated by a constant within each of a set of cells centered at a corresponding set of cluster $`𝐊`$-points. The Green-functions used to calculate the self-energy are coarse-grained or averaged over these cells. This greatly reduces the complexity of the lattice problem to that of a periodic cluster embedded in a host which has to be determined self-consistently. The DCA reduces to the DMFA for a cluster size of one and becomes exact when the cluster size is equal to the system size.

## Results

We present results of DCA calculations for the conventional 2D Hubbard model characterized by an overlap integral $`t`$ between nearest neighbors and $`t^{}`$ between next nearest neighbors and an on-site Coulomb repulsion $`U`$. We set $`t=0.25`$ and choose the magnitude $`\left|t^{}/t\right|0.5`$, so that the band-width $`W=8t=2`$, and study the intermediate coupling regime $`U=W`$. We study the initial corrections to the DMFA by setting the cluster size to $`N_c=4`$, the smallest cluster size which allows for d-wave pairing. We have previously shown that this cluster size is large enough to capture the qualitative low-energy physics of the cuprate superconductors DCA\_Jarrell2 ; DCA\_Maier2 . Due to numerical restrictions we are currently not able to perform systematic studies with increasing cluster size. The $`N_c=4`$ simulations presented in this paper should therefore be interpreted as qualitative extended mean-field results that describe effects on short-length scales (within the cluster).

Our results are summarized in the temperature-doping ($`T`$-$`\delta `$) phase diagram shown in Fig.1. We find regions of antiferromagnetism, d-wave superconductivity, pseudogap and Fermi-liquid like behaviors. In this report we focus on the pseudogap behavior found at low doping $`\delta `$. The pseudogap region is characterized by a suppression of spin excitations moreo  below a crossover temperature $`T^{}`$, defined via the maximum in the spin-susceptibility (see Fig.2c for $`t^{}=0.0`$ and doping $`\delta =0.05`$) when accompanied by the formation of a pseudogap in the density of states for temperatures $`T<T^{}`$ (see Fig.2b). $`T^{}`$ at zero doping is of the order of the magnetic exchange coupling $`J=4t^2/U=0.125`$. To a very good approximation, $`T^{}`$ is equal to the mean-field Neél temperature. These observations indicate that short-ranged antiferromagnetic correlations are a prerequisite to the pseudogap.

Thus, it is natural to explore the scenario in which the antiferromagnetic short-range order causes the pseudogap. To this end we study the effects of lattice frustration, i.e. of the magnitude of the next nearest neighbor hopping integral $`t^{}`$ on the pseudogap. With increasing $`\left|t^{}\right|`$ the spins on the lattice become frustrated and antiferromagnetic correlations are expected to be suppressed. On the other hand, we expect a possible RVB spin-liquid state to be essentially resistant to, or even be stabilized by finite values of $`t^{}`$. In Ref.luca  it was shown that the ground state of the frustrated antiferromagnetic Heisenberg model is almost exactly reproduced by a RVB wave function. In Fig. 2a we plot our results for the antiferromagnetic susceptibility for different values of $`t^{}`$. Clearly antiferromagnetic correlations are suppressed as a function of increasing magnitude of $`t^{}`$.

On the other hand, the density of states plotted in Fig. 2b remains unchanged near the Fermi surface with increasing magnitude of $`t^{}`$. In addition, the temperature $`T^{}`$ where the anomaly in the uniform magnetic susceptibility shown in Fig. 2c occurs, is essentially unaffected by the value of $`t^{}`$. Lattice frustration has little effect on the pseudogap. Thus we conclude that short-ranged antiferromagnetic correlations alone cannot be responsible for the evolution of the pseudogap, inconsistent with scenarios based on the coupling of quasiparticles to antiferromagnetic spin-fluctuations kampf . Moreover, we infer that these results favor RVB physics over antiferromagnetic short-range order as the origin of the pseudogap.

In the ”preformed pairs” scenarios the pseudogap is ascribed to fluctuations that lead to superconductivity at a lower temperature $`T_c`$. In this case we would expect a corresponding signature in the d-wave pair-field susceptibility, i.e. an enhancement due to pairing correlations in the pseudogap region. In the left panel of Fig. 3 we compare the inverse d-wave pair-field susceptibility $`P_d^1`$ of the clean system ($`x=0`$) as a function of reduced temperature $`TT_c`$ for different dopings $`\delta `$ in the pseudogap ($`\delta =0.025,0.05`$) and overdoped ($`\delta =0.20`$) regions. Here we are interested in temperatures below $`T^{}`$ indicated by the arrow for $`\delta =0.05`$. As the doping increases from the pseudogap to the overdoped region, $`P_d^1`$ decreases at fixed reduced temperature indicating that pairing correlations in the overdoped region are even more pronounced than in the pseudogap region. Our result in the inset of Fig.3 for the local cluster cluster  equal-time d-wave pair-field susceptibility $`\overline{P}_d^{local}=\frac{1}{4N_c}_{ijk}g_{ij}g_{ik}c_ic_jc_k^{}c_i^{}`$, where $`g_{ij}`$ is the Fourier-transform of $`\mathrm{cos}k_x\mathrm{cos}k_y`$, supports this finding. This quantity does not contain phase fluctuations and thus can be taken as direct evidence against strong pairing correlations at underdoping. These results lead us to eliminate the ”preformed-pairs” scenarios as an origin for the pseudogap in our results.

Among other scenarios for the pseudogap include the idea of QCP, stripes, and spin-charge separation. In most scenarios, the QCP exists at the zero temperature terminus of $`T^{}`$ versus doping. However, near this point, we find no evidence for a QCP in our simulations; i.e. none of the many susceptibilities that we measure are enhanced as we approach this point. We do not see any signatures of stripes in the charge susceptibility (not shown). However, stripes might be important for larger clusters.

As already discussed, the results for the frustrated system are consistent with a RVB, i.e. spin-charge separated state. In order to address the question of spin-charge separation unambiguously, we would have to calculate the lattice dynamic spin- and charge- structure factors for a range of momenta and extract the spin and charge dispersions. However, these calculations, although formally possible, are presently numerically too expensive to perform. Nevertheless, we believe that we can gain insight into this question by studying the behavior of the corresponding quantities calculated on the cluster onlycluster . Preliminary results (not shown here) for the dynamic cluster spin- and charge susceptibilities at low doping are consistent with spin-charge separation: The spin-susceptibility, dominated by fluctuations at $`𝐐=(\pi ,\pi )`$, becomes suppressed at low frequencies with decreasing temperature well above $`T_c`$. In the charge-susceptibility, which is strongest at $`𝐪=(0,0)`$, weight builds at low frequencies as the temperature is lowered and no sign of a pseudogap is seen.

The results shown in Fig. 4 for the static uniform lattice spin- (circles) and charge-susceptibilities (squares) versus temperature at underdoping ($`\delta =0.05`$, left) and overdoping ($`\delta =0.30`$, right) support this picture. As already discussed, the spin-susceptibility $`\chi _s`$ in the underdoped system shows an anomaly at a temperature $`T^{}`$ below which spin-excitations are suppressed. The charge-susceptibility $`\chi _c`$ in the underdoped system however displays qualitatively different behavior: When the spin degrees of freedom become suppressed at $`TT^{}`$, $`\chi _c`$ starts to rise with decreasing temperature. This clearly indicates that coherence in the charge channel starts to set in with decreasing temperature whereas the spin excitations become suppressed. These results are consistent with recent angle-resolved photoemission experiments Ding2  where superconducting order in the underdoped system was found to be accompanied by an emerging quasiparticle coherence as seen in the spectra.

Thus, our results for the underdoped region can be interpreted within an RVB, spin-charge separated picture, where the development of the pseudogap originates in the pairing of the spin degrees of freedom and superconductivity is driven by coherence of the charge excitations. In the overdoped region (right panel in Fig. 4) our results show that the system becomes more conventional, i.e. spin and charge excitations behave in a similar way.

## Summary

We have investigated the microscopic origin of the pseudogap in DCA simulations of the underdoped 2D Hubbard model. We find that short-ranged correlations alone, although necessary for the pseudogap to emerge, is not sufficient to describe its origins. Similarly, although stripe fluctuations may emerge for larger clusters, no indications were found in our calculations, and therefore stripes are not the origin of the pseudogap. Likewise, we find no indication for strong pairing fluctuations for temperatures $`T_c<T<T^{}`$, in disagreement with pre-formed pairing scenarios. Finally, although we cannot exclude a quantum critical point at finite doping as an origin, we find no evidence for it. Instead, all of our numerical evidence points to spin-charge separation as described by Anderson’s RVB theory as the microscopic origin of the pseudogap. In particular, in the pseudogap region, spin degrees of freedom are suppressed while charge excitations become more coherent with decreasing temperature. At overdoping where no pseudogap is found, spin and charge behave qualitatively similar.

## Acknowledgments

This work was supported by NSF grant DMR-0073308 and by NSF cooperative agreement ACI-9619020 through computing resources provided by the NPACI at the Pittsburgh Supercomputer Center.


format: text

text:
# Potential and current distribution in strongly anisotropic Bi2Sr2CaCu2O8 single crystals at current breakdown

## I Introduction

Understanding the current and the potential distributions at the onset of dissipation in a superconductor is important; at a fundamental level for interpreting the results of transport measurements in terms of the force-velocity relation for the vortices and the Josephson coupling between layers as well as on a practical level for understanding how to maximize the current carrying capacity of a wire. The family of high-$`T_c`$ cuprate superconductors naturally raises the question of the influence of high anisotropy on current carrying properties. How does one interpolate between the extreme anisotropic limit, where current injected into a superconducting plane has no transfer to other planes, and the isotropic case where the penetration of the current in the superconducting state is determined by the London screening length.London (1960) To give a concrete illustration of the problem, a naive interpretation of the threshold current for dissipation in a strongly anisotropic $`\mathrm{Bi}_2\mathrm{Sr}_2\mathrm{CaCu}_2\mathrm{O}_8`$ (BSCCO) single crystal based on penetration limited to the London screening length or, even more naively, supposed to be uniform over the typical thickness of a few micrometers gives threshold current densities one to two orders of magnitude lower than a standard interpretation of a magnetic hysteresis loop.

Anisotropy is a central feature of the family of the high-$`T_c`$ cuprate superconductors which are usually modeled by discrete superconducting sheets parallel to the $`ab`$ crystallographic plane weakly coupled together in the $`c`$ direction by extended Josephson junctions.Lawrence and Doniach (1971) The weak coupling between planes gives rise to very high conduction anisotropy in both the normal and superconducting states. It is also responsible for the richness of the phase diagram of the vortices created when one applies a magnetic field. A convenient experimental probe for investigating these phases, particularly as regards their pinning to the host lattice disorder, is to look at the voltage response to transport current ($`VI`$ characteristics) for currents which take one into the resistive regime where the vortices are dislodged. But to interpret these measurements correctly, it is important to know about the current distribution in the resistanceless regime, how it is modified as dissipation sets in and of course what the potential distribution then looks like.

It was remarked some time ago,Busch et al. (1992) in connection with resistivity experiments, that if the response of the superconductor is Ohmic, $`𝐄\mathbf{=}𝝆\mathbf{}𝐉`$ (the resistivity tensor $`𝝆`$ is constant and independent of current and position), the problem may be treated as for the normal state where it is implicit that the phase slip giving rise to the voltage also relaxes the magnetic dephasing to allow the current to penetrate beyond the London screening depth. In this regime, at least, the high anisotropy of the (BSCCO) samples studied limited the current penetration into the depth of the sample to a small fraction of its thickness and that it would therefore be quite erroneous to deduce a critical current density on the basis of uniform penetration across the section of a short sample. Experiments to date on the current profile have been able to show a concentration of current towards the sample edges, but have been insensitive to the depth dependenceFuchs et al. (1998) which is the aspect that interests us chiefly here.

The present paper reports an experimental investigation and a few basic reflections on this problem in the non-Ohmic, low-temperature, high magnetic-field vortex solid phase, also in monocrystalline BSCCO in a c-directed magnetic field. The experiment consists of measuring the voltage response to short pulses of current up to and beyond dissipative breakdown. Potentials are measured at contacts placed in the usual way on the ab faces supplemented by in-depth potential contacts to have direct access to depth-dependent features. One of the important points that emerges is that the physically unrelated $`c`$-axis and $`ab`$-plane dissipative breakdowns interact to create a resistive/nonresistive front which has the effect of altering the current distribution and giving rise to a critical current which involves both $`ab`$ and $`c`$ characteristics. Not only does the voltage response along an $`ab`$-plane surface into which the current is injected and withdrawn (the “$`ab`$ configuration” of Fig. 1) show features of the Josephson junctions between planes, but also the response to current injection and withdrawal on opposing $`ab`$-plane faces (the “$`c`$ configuration” of Fig. 1) shows features of breakdown in the $`ab`$ plane. One manifestation of this is that the threshold current in the $`ab`$ configuration on BSCCO single crystalsSas et al. (2000) shows the same temperature and preparation dependence as for the $`c`$ configuration.Rodríguez et al. (1993) This brings up the question as to who measures what and does one measure what one thinks; but beyond that what do the results mean?

## II Experimental technique

Our samples are thin monocrystalline platelets with the $`c`$ axis perpendicular to the face. As all the samples studied gave qualitatively the same results for similar contact configurations, we limit ourselves here to the results on one sample of each different contact configuration. The sizes of samples $`A`$, $`B`$, and $`C`$ are about $`2\times 0.35\times 0.005`$ $`\mathrm{mm}^3`$, $`1.1\times 0.3\times 0.002`$ $`\mathrm{mm}^3`$, and $`0.7\times 0.5\times 0.005`$ $`\mathrm{mm}^3`$. They were all fabricated by a melt cooling technique.Keszei et al. (1989); Cooper et al. (1990) The critical temperatures were between 88 and 90 K with a transition width of about 2 K at zero field. The anisotropy coefficient $`\gamma _n500`$ was estimated from normal-state resistivities with $`\rho _{\text{ab}}`$ 100 $`\mu \mathrm{\Omega }`$ cm at 90 K. In the case of samples $`A`$ and $`B`$, electrical contact was made by bonding 25-$`\mu `$m gold wires with silver epoxy fired at 900 K. For sample $`C`$ the contacts were made by depositing silver on a lithographically defined area using a lift-off technique and heat treating at 700 K for 1000 s. In all cases the contact resistance was less than 3 $`\mathrm{\Omega }`$. The sample was mounted flat against a 7-mm-diameter 0.5-mm-thick sapphire disk with silicone grease to ensure thermal homogeneity and mechanical freedom and placed in the bore of a superconducting magnet with the field along the $`c`$ direction. The sample and disk were surrounded by exchange gas and the temperature was electronically regulated. The longitudinal voltage-current ($`VI`$) characteristics were measured by a symmetrized differential four-point technique using 25-$`\mu `$s triangular pulses (12.5 $`\mu \text{s}`$ from zero to maximum current) of maximum amplitude in the range 10 mA $`<I_m<`$ 350 mA, usually restricted to exceed the threshold by about 30%, with repetition time $`0.1`$ s. As reported earlier,Sas et al. (2000) the observed threshold is independent of pulse duration for times up to about 250 $`\mu \text{s}`$, argued there to be the time required for heat produced in the current contact to diffuse towards the voltage contact, whereas bulk heating in the superconductor was estimated not to be important. Further technical details may be found in Ref. Sas et al., 2000.

## III Experimental results

We present first the results of experiments on the $`VI`$ response with surface contacts in the standard $`ab`$ and $`c`$ configurations on the same sample using contacts common to the two configurations. The results lead rather naturally to a description in terms of resistive/nonresistive fronts moving outwards from the current contacts. Experiments are then presented for other surface contact configurations intended to check this hypothesis. These are followed by experiments with a different type of contact configuration designed to probe the depth of the sample as breakdown progresses.

The first experiment was done on samples contacted to be able to measure in both $`ab`$ and $`c`$ configurations on the same sample at the same time as illustrated in Fig. 2. The “$`ab`$ configuration” measurements refer to the potential drop between contacts ($`A2,A4`$) for current injected at $`A1`$ and withdrawn from $`A5`$. The “$`c`$ configuration” refers to the potential measured across ($`A2,A7`$) for current passed through ($`A1,A6`$). Figure 3 shows the temperature dependence of the threshold current in the $`VI`$ characteristic at $`1.5`$ T in sample $`A`$ for both configurations. The threshold current is defined by the break in slope criterion of Ref. Sas et al., 2000 which gives essentially the same values as the criterion of crossover from Kim-Anderson behavior at low current to power law close to linear behavior at higher current,Portier et al. (2002) illustrated in detail for the first pair of curves of Fig. 2 in the semilog plot in the same figure. The measured thresholds, as represented in Fig. 3, are seen to be independent of the configuration, whether the sample is field cooled (FC) or zero-field cooled (ZFC) prepared, despite the fact that the two preparations give different thresholds for temperatures lower than the peak in the ZFC results.Sas et al. (2000) For the FC preparation, the field is applied above $`T_c`$ and the sample is cooled at constant field to the lowest measuring temperature, subsequent measurements being made on increasing the temperature at the same field. For ZFC preparation, the system is cooled in zero field from above $`T_c`$ to the lowest temperature before applying the field and then making measurements at sequentially higher temperatures at the same field. Figure 4 shows the magnetic-field dependence in both contact configurations for ZFC preparation at $`5`$ K on increasing the field (FC preparation followed by field variation gives the same result, except for the initial FC prepared point Portier et al. (2002)). At fields above about $`0.3`$ T, the threshold currents for the two different configurations again show the same behavior, this time in field, varying approximately as $`H^{1/2}`$. Not only are the temperature and field dependences of the threshold currents in the two configurations similar in form, but the values themselves coincide. Measurements along several other lines in the $`(H,T)`$ plane confirm this indistinguishability of $`ab`$ and $`c`$ configuration thresholds as generic behavior for the low-temperature high-field domain. This fact led us to look in more detail with the following experiment.

Again referring to the inset of Fig. 2, we used three potential contacts ($`A2,A3,A4`$) on the top layer and two ($`A7,A8`$) on the bottom. The contact $`A2`$ is half as far away from the left current contact $`A1`$ as is the contact $`A4`$ from the right-hand current contact $`A5`$, while $`A3`$ is midway between the two current contacts. If we apply the current between contacts $`A1`$ and $`A5`$ ($`ab`$ direction) we measure a much lower threshold current between $`A2`$ and $`A4`$ than between $`A3`$ and $`A4`$. On the other hand, we find virtually identical $`VI`$ curves on the $`ab`$ potential contacts ($`A2,A4`$) whether the current is applied in the $`c`$ direction through contacts ($`A1,A6`$) or along the $`ab`$ direction ($`A1,A5`$); similarly the $`VI`$ across ($`A3,A4`$) is nearly identical for $`ab`$-directed current through ($`A1,A5`$) or for $`c`$-directed current through ($`A5,A9`$).

These results suggest the idea that the sample first becomes resistive around the injection points and that this resistive region progressively invades the sample as the current is increased. This led us to experiment with the contact configuration of sample $`B`$ illustrated in Fig. 5 where we drove the current through the $`c`$ direction and also measured the voltage drop across that direction. On sending the current through ($`B1,B4`$) the $`VI`$ characteristic measured across ($`B2,B5`$) shows a threshold current which is about half that measured across contacts ($`B3,B6`$). On the other hand, if we put the current through the middle contacts ($`B2,B5`$) we find the same threshold at either of the two end potential contacts ($`B1,B4`$) or ($`B3,B6`$) with a value intermediate between the previous two.

Although these experiments lend considerable support to the idea of a resistive/nonresistive front progressing along the surface, they say nothing about its profile with depth since at no point were we able, in this low-temperature high-field phase and despite micron thin samples, to drive the lower face contacts resistive with current sent along the top face. To have more direct information about the shape of the front with depth, we prepared a third sample $`C`$ as in Fig. 6 with a lithographically defined terrace argon ion etched along the length of one $`ab`$ face to a depth of $`220`$ nm. The upper part of the face was also argon ion etched to avoid any difference in surface pinning between it and the terrace, the width of which was restricted to a small fraction ($`20\%`$) of the total sample width to minimize the trivial effect of current spread as the distribution attains the depth of the terrace. The contacts were accurately aligned between the two levels by a scanning electron-beam masking microscope. On a control sample with the same contact configuration but with no terrace, the $`VI`$ response on contacts ($`C2,C3`$) gave the same values of the threshold current and dynamic resistance as contacts ($`C5,C6`$), themselves very similar to those for the top contacts of the terraced sample.

Sending an $`ab`$-directed current through contacts ($`C1,C4`$) on the upper part of the face, we measured the potential drop across ($`C2,C3`$) on the upper part and across ($`C5,C6`$) on the step. The threshold current for the terrace contacts is considerably higher than for the corresponding top contacts: about $`25`$ mA and $`10`$ mA, respectively, much greater than the decrease of $`20\%`$ in the current density to be expected from simple current spread into the wider part of the sample. We interpret the large difference to arise from the current distribution with depth.

## IV Interpretation of experimental results

All these observations are consistent with the idea of a resistive front moving progressively outward from the current contacts as illustrated for sample $`A`$ in Fig. 7. Also, and more surprisingly, it seems that the position of the front depends only on the magnitude of the current injected and seems to be independent of its ultimate destination if one is to account for the nearly identical $`VI`$ response on, for example, the potential contact pair ($`A2,A4`$) independently of whether the current is withdrawn on the same face ($`A1,A5`$) or the opposite face ($`A1,A6`$); similarly for the potential pair ($`A3,A4`$) for current through ($`A1,A5`$) or ($`A5,A9`$). The difference between the values of the threshold current measured on ($`A2,A4`$) and ($`A3,A4`$) contact pairs is accounted for by the distances of the potential contact from the nearest current contact ($`A2`$ from $`A1`$ or $`A4`$ from $`A5`$). In these experiments, the fronts did not attain contact $`A3`$ which remained at the potential of the resistanceless region. The triangular form used for the resistive region in the figure is purely schematic to show the progression of the front. Further consideration of possible forms is given below where it is argued that the form is in fact most probably rectangular.

Figure 8 represents how the resistive front is imagined to progress in the experiment with sample $`B`$ when the current is run through the $`c`$ direction. Here too the lower threshold current was obtained for potential pair ($`B2,B5`$) nearest the current contacts. The resistive front from current contacts ($`B1,B4`$) attains ($`B2,B5`$) before ($`B3,B6`$). When the current was run through the middle of the sample by ($`B2,B5`$) the injected current is divided between two directions with a corresponding reduction in current density so that the voltage drop appears at higher current for the same distance and gives the same response at both ends.

Figure 9 represents the front for the experiment with the terraced sample. Here the form of the front has been taken to be rectangular for the reasons outlined in the following section. As the front is defined by the line at which the current density reaches the threshold for resistive breakdown ($`ab`$ plane and/or $`c`$ direction) the integral of the threshold current density crossing the front is just the total current. Furthermore if the front does not change geometrical form with current amplitude its dimensions increase linearly with the current. Then, knowing the distance from the nearest current contact and the depth of the terrace, the threshold current values for the upper and step contacts give information on the depth and potentially some information on the form of the front as it attains the upper and terrace level contacts. We can expect the aspect ratio of the front to be related to the anisotropy.

## V Elementary Understanding

We shall assume for the present discussion that edge effects are not a dominant factor for our observations. Although we expect an integrable singularity at the edge of the sample, a large part of the middle of the sample has an essentially flat distribution with distance from the center line 222 Fuchs and collaborators (Ref. Fuchs et al., 1998) have done a very interesting experiment to determine the distribution of the $`\widehat{𝐜}`$-directed component of the current-induced magnetic field across and just outside a current carrying sample. In the low-temperature phase which interests us here, they deduce that there is a stronger current density on the edge of the sample. Their observations are compatible with the expected finite width correction $`j_x=[1(2y/w)^2]^{1/2}`$, where $`w`$ represents the full width and $`y`$ is measured from the midline. This distribution is constructed to ensure the zero $`\widehat{𝐜}`$-directed induced field in the sample as required by the superconducting state with fixed vortices. and we deal in the first instance uniquely with the distribution with depth $`\widehat{𝐳}`$ and along the length $`\widehat{𝐱}`$ of the sample. We know that in the normal resistive state the current distribution is governed by the anisotropy of the resistivity. By rescaling the coordinates according to the anisotropy factor $`\gamma _n=\sqrt{\rho _c/\rho _{ab}}`$ the problem can be solved as a Laplace equation in the potential. The necessary coordinate rescaling is determined by the anisotropy factor $`\gamma `$: $`z\gamma z`$. In the normal state, it is clearly given by $`\gamma _n^2=\rho _c/\rho _{ab}(500)^2`$ and would give rise to an $`ab`$ current penetration depth of about $`p_{ab}\mathrm{}_{ab}/\gamma _n`$ $`1`$ $`\mu \mathrm{m}`$ for a distance $`\mathrm{}_{ab}=1`$ $`\mathrm{mm}`$ between $`ab`$ current contacts.Busch et al. (1992)

It is less recognized, perhaps, that the low current distribution in the nonresistive superconducting regime can be governed by the same sort of equation for the superconductor phase in similarly rescaled coordinates if the penetration depth is smaller than the London screening depth. In the nonresistive superconducting phase, the Laplacian in the potential is replaced by a Laplacian in the superconductor phase which has the same form of relationship to the superconducting current as does the potential to the current in the normal state. The anisotropy ratio in the superconducting state is given by $`\gamma _s^2=(n_se\mathrm{}/m_{ab})/sJ_J`$, the ratio of the coefficients of the $`ab`$ and $`c`$ phase gradients in the relation between the current operator and the phase. The $`sJ_J`$ term is the low current expansion of the Josephson relation between successive superconducting planes separated by $`s`$, while $`n_s`$ is the three-dimensional (3D) superconducting electron density and $`1/𝐦`$ the inverse mass tensor. If one were to neglect London screening, $`\gamma _s`$ would govern the current penetration in the same way as does $`\gamma _n`$ in the normal state.

The superconducting anisotropy is usually expressed as a mass ratio or London screening length ratio for $`c`$ and $`ab`$ currents $`\gamma _s^2=m_c/m_{ab}=\lambda _c^2/\lambda _{ab}^2`$ and is customarily expected to have about the same value as the extrapolated normal-state resistivity ratio. In BSCCO, however, experimentsLatyshev and Yamashita (1999) on mesa structures have shown a Josephson critical current reduced by a factor of about $`30`$ with respect to the usual Ambegaokar-BaratoffAmbegaokar and Baratoff (1963) relation to the high current resistance; the low current resistive slope is proportionately altered to $`V_{\mathrm{gap}}/sJ_J30\rho _c`$ and the coefficient in the Josephson relation $`J_z=J_J\mathrm{sin}(\phi _n\phi _{n+1})sJ_J\phi /z`$ at low current is reduced by the same factor with a corresponding enhancement of the anisotropy factor \[$`\phi _n`$ is the phase $`\phi (𝐫)`$ of the superconductor wave function at the $`n\mathrm{th}`$ superconducting plane\]. This has been attributed to the $`d`$-wave nature of the order parameter and already enhances the anisotropy factor by $`\sqrt{30}`$ over the normal state. But if one also considers the effect of phase fluctuations across the Josephson junction caused by misalignment of vortices in magnetic field, another multiplicative factor $`\mathrm{cos}\phi _{n,n+1}^{1/2}`$ should be incorporated into $`\gamma _s`$. The factor $`\mathrm{cos}\phi _{n,n+1}`$ has been variously estimated from Josephson plasma resonance Matsuda et al. (1995) experiments to be between about $`10^3`$ and $`10^1`$, extrapolating to about $`4\times 10^3`$ (Refs. Shibauchi et al., 1999 and Gaifullin et al., 2000) at low temperature and $`\mu _0H=1.5`$ T, which suggests a very strongly enhanced anisotropy parameter in the superconducting phase as compared with the normal state: $`3\times 10^4<\gamma _s<1\times 10^6`$. Current penetration in the superconducting phase is governed by a combination of London screening and anisotropy, both of which act to limit the penetration. But, even before including the effects of London screening, the anisotropy limits the penetration depth $`p_{ab}\mathrm{}_{ab}/\gamma _s`$ in the middle of a 1-mm-long BSCCO sample to $`2<p_{ab}<20\mathrm{nm}\lambda _{ab}200`$ nm. Thus in short samples of BSCCO like ours the current penetration in the superconducting phase is always dominated by anisotropy rather than by London screening.

This suggests that the effective thickness for current flow might be given approximately by $`p_{ab}\mathrm{}_{ab}/\gamma _s<\lambda _{ab}t`$ (the real thickness), and furthermore that the current density is nonuniform along the sample, being higher closer to the contacts, causing breakdown to resistive behavior to be progressive and not to penetrate immediately throughout the bulk. Under these conditions, $`ab`$ breakdown alone is not possible: if only $`ab`$ breakdown were to occur, the front would have to be parallel to the $`c`$ direction because no part of the interface with the nonresistive region can be parallel to the electric field produced along the $`ab`$ breakdown direction. But if breakdown does not occur uniformly throughout the thickness of the sample, the front must somewhere change orientation. The electric field in the resistive portion, which must be oriented perpendicular to the interface, therefore also changes orientation and that requires $`c`$ breakdown. Hence a resistive/nonresistive front within the sample which does not go straight through it along a principal direction necessarily implies simultaneous $`ab`$ and $`c`$ breakdown. The regions of $`ab`$ and $`c`$ breakdown need not be spatially coincident however, but they must be at least contiguous. The same can be argued for the case of pure $`c`$ breakdown.

Another consequence of the nonuniformity of the current distribution in depth is the presence of shear forces on the vortices and the possibility of sliding between planes of (semi)ordered vortex segments.Busch et al. (1992); Khaykovich et al. (2000) This has important consequences on the locality of the $`𝐄(𝐉)`$ relation. If the force required to shear two adjoining planes derives from an effective Josephson coupling $`J_J\mathrm{cos}\phi _{n,n+1}`$, the condition for nonelastic shear (slip) is $`J_{ab}/zJ_J\mathrm{cos}\phi _{n,n+1}a/s^2,`$ where $`a`$ and $`s`$ are the distances between vortices in the plane and the distance between neighboring planes, respectively. Using the expression for the anisotropy penetration depth $`p_{ab}\mathrm{}_{ab}/\gamma _s`$ and identifying the effective Josephson coupling with the $`c`$-axis critical current density $`J_c^{\mathrm{th}}`$, the condition for shear to occur at the onset of dissipation is $`J_{ab}^{\mathrm{th}}/J_c^{\mathrm{th}}\gamma a\mathrm{}_{ab}/s^2`$. At high values of anisotropy, $`\gamma _s10^2`$, however, the shear strength is dominated by magnetic coupling between vortex segments.Buzdin and Feinberg (1990); Koshelev and Kes (1993) Either by relating the critical force per vortex segment for shear slip $`f_c`$ to the tilt modulus $`C_{44}`$ by $`f_cC_{44}as`$ and using the evaluations of the modulus for magnetic coupling,Koshelev and Kes (1993) or using directly the calculation of the shear strength,Pe et al. (1997) it appears that shear slip between planes should occur for $`J_{ab}/z(c\varphi _0/32\pi ^3\lambda _{ab}^4)\mathrm{cos}\phi _{n,n+1}a/s`$. It will be seen that this condition is met with the values obtained from the data analyzed according to the reasoning outlined below and it is a necessary condition to have a local $`𝐄(𝐉)`$ relation.

The above considerations lead to a scenario where the current injected and withdrawn from the same $`ab`$ surface is confined by anisotropy to a penetration depth $`\mathrm{}_{ab}/\gamma _s`$ in the nonresistive regime with consequent enhancement of the current density. When the current density attains the critical value, for either $`J_c`$ or $`J_{ab}`$, the current distribution is modified until both reach their breakdown values and the sample is resistive on one side and nonresistive on the other side of a front which advances out from each current contact point as the current is increased. As the fronts attain successive voltage contacts, these display a potential drop with respect to the portion of the sample which has remained in the nonresistive state. The shape of the resistive region can be expected to be a function of the $`𝐄`$-$`𝐉`$ response around the breakdown values (vortex depinning for $`ab`$ current and Josephson-junction behavior for $`c`$ current) described for the most part by the anisotropy factor $`\gamma `$ of the resistive state and $`\mathrm{\Gamma }=J_{ab}^{\mathrm{th}}/J_c^{\mathrm{th}}`$. If the shape were determined uniquely by the $`𝐄`$-$`𝐉`$ response, its position would depend only on the contact from which it originated and the magnitude of the current injected.

Restrictions are imposed on the shape of the front by the irrotational nature of the electrochemical field and by current conservation. The inductive field from the current is negligible compared to resistive fields for our pulses and sample sizes and we can suppose the field $`𝐄`$ to be irrotational: $`\times 𝐄=𝐁/t=0.`$ We can then, as in London’s treatment of breakdown in a type-I current carrying wire,London (1960) demand that $`𝐄\mathrm{d}`$$`\mathbf{}`$ $`=0`$ be satisfied through the resistive/nonresistive interface (Fig. 10) with the consequence that the interface must be perpendicular to the field. As argued above, our quasipoint contact situation enforces that both $`ab`$ and $`c`$ breakdown appear simultaneously and the $`VI`$ response in either configuration will contain both $`c`$ and $`ab`$ features. Either $`ab`$ and $`c`$ breakdown occur together at the interface or they occur alone in separated but necessarily contiguous or overlapping regions. The first case is considered in part (a) of Fig. 10 where a segment of the interface is locally oriented at an angle $`\mathrm{\Theta }`$ with the $`ab`$ plane. A model with a local (locality implicitly assumes shearing of the vortex segments from plane to plane) $`𝐄`$-$`𝐉`$ relation for $`𝐉>𝐉^{\mathrm{th}}`$ of the form

$$E_x=0,|J_x|J_{ab}^{\mathrm{th}}$$

$$E_x=\pm E_{ab}^{\mathrm{th}}+\rho _{ab}^{}(J_xJ_{ab}^{\mathrm{th}}),J_x\pm J_{ab}^{\mathrm{th}}$$
(1)
$$E_z=0,|J_z|J_c^{\mathrm{th}}$$

$$E_z=\pm E_c^{\mathrm{th}}+\rho _c^{}(J_zJ_c^{\mathrm{th}}),J_z\pm J_c^{\mathrm{th}}$$
(2)
on the resistive side is shown in part (b) of the figure. $`\rho _{ab}^{}`$ and $`\rho _c^{}`$ are the dynamic resistivities in the $`ab`$ and $`c`$ direction in the superconductive state beyond breakdown. The condition that the field be perpendicular to the interface would demand that $`\mathrm{tan}\mathrm{\Theta }=\mathrm{\Delta }z/\mathrm{\Delta }x=\pm E_{ab}^{\mathrm{th}}/E_c^{\mathrm{th}}`$, where the $`E_{ab}^{\mathrm{th}}`$ and $`E_c^{\mathrm{th}}`$ components of the electrochemical field both point towards the interface, in the directions of the $`ab`$ and $`c`$ components of the current flow. However because both components of the current density must diminish on approaching the interface, $`J_x/x`$ and $`J_y/y`$ both have the same sign whereas current conservation $`𝐉=0`$ requires them to have opposite signs. More explicitly, if the interface is imagined to have positive intercepts with the axes of an ($`\widehat{x},\widehat{z}`$) coordinate system whose origin is at $`P`$ in the resistive region as in Fig. 10, both $`J_x>J_{ab}^{\mathrm{th}}`$ and $`J_z>J_c^{\mathrm{th}}`$ (all with positive values, current flowing from left to right and top to bottom). But the fact that the current diminishes on approaching the interface along each of the principal directions to reach $`J_{ab}^{\mathrm{th}}`$ and $`J_c^{\mathrm{th}}`$ which defines the interface imposes that both $`J_x/x<0`$ and $`J_z/z<0`$ at the same point $`P`$, in violation of current conservation which demands that $`J_x/x=J_z/z`$.

We conclude that a region of spatially coincident $`ab`$ and $`c`$ breakdown may not have an interface with the nonresistive region.

Resistive/nonresistive interfaces perpendicular to pure $`ab`$ or pure $`c`$ breakdown, however, do not violate current conservation. For sole $`c`$ breakdown, $`J_z>J_c^{\mathrm{th}}`$, $`J_x<J_{ab}^{\mathrm{th}}`$ at the interface and the positive $`J_x/x`$ which must result from the negative $`J_z/z`$ to conserve current has no consequence for the electric field because the $`x`$ component is subcritical and produces no electric field. A similar argument can be made for pure $`ab`$ breakdown, establishing that nonresistive/resistive interfaces are possible for single-component breakdown whereas they are not for spatially coincident two-component breakdown. Furthermore it is possible to change the orientation of the interface between the two perpendicular principal directions of the superconductor on the condition that the two separated mutually perpendicular $`ab`$ and $`c`$ fronts meet at a common point and that the interface between the two types of breakdown converges to the same point. The latter interface could be either a simple common boundary or more generally a lens shaped region of coincident $`ab`$ and $`c`$ breakdown. Such a situation is illustrated in part (c) of Fig. 10. We conclude that the fronts must be formed of segments of interfaces perpendicular to pure $`ab`$ or pure $`c`$ current breakdown. Numerical solutions of the modelKriza et al. () (unpublished) indicate a simple rectangle with a generic form of the type illustrated in Fig. 11.

The length $`\mathrm{}`$ and depth $`d`$ of the resistive region are related to the current by

$$I=w(dJ_{ab}^{\mathrm{th}}+\mathrm{}J_c^{\mathrm{th}}),$$
(3)
where $`w`$ is the width of the sample. The threshold current measured for the potential between the nonresistive region and a contact on the surface at a distance $`\mathrm{}_1`$ from the current injection contact is given by the condition that the front arrives at the contact, $`\mathrm{}=\mathrm{}_1`$:

$$I_\mathrm{}_1^{\mathrm{th}}=w\mathrm{}_1(J_{ab}^{\mathrm{th}}d/\mathrm{}+J_c^{\mathrm{th}})$$

and for a voltage contact at depth $`d_1`$ by $`d=d_1`$:

$$I_{d_1}^{\mathrm{th}}=wd_1(\mathrm{}/d)(J_{ab}^{\mathrm{th}}d/\mathrm{}+J_c^{\mathrm{th}}),$$

provided that $`\mathrm{}(I)`$ is sufficient that the front attains the contact in the $`\widehat{x}`$ direction ($`I_\mathrm{}_1^{\mathrm{th}}>I_{d_1}^{\mathrm{th}}`$). In general, $`\mathrm{}=\mathrm{}(I)`$ and $`d=d(I)`$, but if there are no other relevant lengths in the problem (distances between $`ab`$ configuration current contacts $`\mathrm{}_{ab}`$ $`\mathrm{}`$ and $`c`$ configuration current contacts $`\mathrm{}_c=td`$) one can expect the ratio $`\mathrm{}/d`$ to be independent of the current and only a function of the anisotropy factors $`\gamma `$ and $`\mathrm{\Gamma }=J_{ab}^{\mathrm{th}}/J_c^{\mathrm{th}}`$. If we interpret the top and terrace thresholds of the experiment on sample $`C`$ as measurements of $`I_\mathrm{}_1^{\mathrm{th}}`$ and $`I_{d_1}^{\mathrm{th}}`$ and we make the hypothesis that the aspect ratio $`\mathrm{}/d`$ is independent of current, then the measurements tell us that

$$\mathrm{}/d=(\mathrm{}_1/d_1)(I_{d_1}^{\mathrm{th}}/I_\mathrm{}_1^{\mathrm{th}})2500$$

and

$$J_{ab}^{\mathrm{th}}d/\mathrm{}+J_c^{\mathrm{th}}10\mathrm{A}\mathrm{cm}^2$$
(4)
or, multiplying by $`\mathrm{}/d`$,

$$J_{ab}^{\mathrm{th}}+J_c^{\mathrm{th}}\mathrm{}/d2.5\times 10^4\mathrm{A}\mathrm{cm}^2,$$
(5)
where the last two quantities are the same combination of $`ab`$ and $`c`$ transport properties which we cannot separate without a model calculation for the aspect ratio of the resistive front. This will be the subject of a future paperKriza et al. () (unpublished) on numerical solutions of the present model.

The upper limits on $`J_c^{\mathrm{th}}`$ and $`J_{ab}^{\mathrm{th}}`$ given by relations (4) and (5) are consistent with the value $`J_c^{\mathrm{th}}2`$ A cm-2 obtained from combining $`J_J`$ measured in the mesa experiment with $`\mathrm{cos}\phi _{n,n+1}4\times 10^3`$ extrapolated from the Josephson plasma resonance experiments Shibauchi et al. (1999) and the value $`J_{ab}^{\mathrm{th}}2\times 10^4`$ A cm-2 estimated from the magnetic hysteresis.Khaykovich et al. (2000) The condition for slip between planes is well satisfied for these values indicating that the description is self-consistent.

## VI Summary

The experimental results on $`VI`$ characteristics for $`ab`$ surface contacted samples of the highly anisotropic BSCCO superconductor in the low-temperature vortex lattice phase above about $`2000`$ Oe show that the $`VI`$ responses for $`ab`$ or $`c`$ current configurations are virtually identical and that the dissipative region invades the samples in the form of a resistive/nonresistive front moving out from the current contacts with increasing current. The $`VI`$ characteristic depends only on current magnitude, voltage contact position and nearest current injection point. It is argued that current penetration in the nonresistive regime is limited by the high anisotropy to depths less than the London screening length, and that the resistive region must involve simultaneous $`ab`$ and $`c`$ breakdown and shear of successive planes of vortex segments. The shape of the resistive region is argued to be composed of rectangular segments and seems to be a simple rectangle with interfaces to the nonresistive portion consisting of pure $`ab`$ or pure $`c`$ breakdown. A different contact arrangement using an ion etched terrace to sample the potential in the depth of the sample brings extra information that allows access to the aspect ratio of the rectangle. But because the $`ab`$ and $`c`$ configuration experiments measure the same combination of $`ab`$ and $`c`$ properties, it is necessary to propose and solve a specific model to separate them. A detailed discussion of numerical solutions of the model has been reserved for future publication.

###### Acknowledgements.

We take pleasure in acknowledging discussion with I. Tüttő and L. F. Kiss. We are grateful to L. Forró and to B. Keszei for supplying the crystals used in this study. B.S. and Á.P. acknowledge the Atomic Energy Commission for a Bursary to visit the Saclay laboratory. Research in Hungary has been supported by Grant Nos. OTKA T037976 and TS040878. We are grateful also for the support of the “BALATON” Collaboration program of the French and Hungarian Foreign Affairs Ministries.


text:
# Energy loss in high energy heavy ion collisions from the Hydro+Jet model

## Abstract

We investigate the effect of energy loss of jets in high energy heavy ion collisions by using a full three-dimensional space-time evolution of a fluid combined with (mini-)jets that are explicitly evolved in space-time. In order to fit the $`\pi ^0`$ data for the Au+Au collisions at $`\sqrt{s_{NN}}=130`$ GeV, the space-time averaged energy loss $`dE/dx(\tau \text{3 fm/}c)=0.36`$ GeV/fm is extracted within the model. It is found that most energy loss occurs at the very early time less than 2 fm/$`c`$ in the QGP phase and that energy loss in the mixed phase is negligible within our parameterization for jet energy loss. This is a consequence of strong expansion of the system.

Measurements of high $`p_\mathrm{T}`$ hadrons at the Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC) may provide insight into the excited partonic matter, often called a quark gluon plasma (QGP), produced in heavy ion collisions QMproc . Because jets have to traverse the excited matter, their spectra should be changed compared to the elementary hadron-hadron data. This energy loss of jets has been proposed as a possible signature of the QGP phase Gyulassy:1990ye . Over the past year, a lot of work have been devoted to study the propagation of jets through QCD matter Wiedemann ; Zakharov ; Baier ; Levai .

Recently, hadronic transverse momentum distribution in Au + Au collisions at $`\sqrt{s_{NN}}=130`$ GeV has been measured at RHIC and found that these spectra show the depletion at high transverse momentum phenix\_pi0 . Comparison of the data with pQCD parton model based calculations Wang:2001gv ; Fai ; Levai2  and a transport theoretical study Nara:2001zu  shows the indication of energy loss in central Au + Au collisions at RHIC energies. Phenomenological studies Wang:2001gv ; jamal  suggest that the space-time averaged energy loss yields small values of $`dE/dx0.25`$$`0.3`$ GeV/fm. In addition, jet quenching is found to be the main contribution on the elliptic flow parameter in the large transverse momentum region v2jet1 ; v2jet2 . Jet energy loss is modeled by modifying the fragmentation functions in pQCD parton model and usually neglecting the dynamical effects.

It might be important to consider dynamical effects on the jet interactions. The effects of the expansion of the system on the energy loss and the azimuthal asymmetry have been studied v2jet1 ; GVWH ; EWANG . Hydrodynamics provides the space-time evolution of thermalized partonic/hadronic matter produced in heavy-ion collisions KHHH ; Teaney:2000cw ; Hirano:2001eu ; Hirano:2001yi . Recently, one of the authors (T.H.) developed a fully three-dimensional hydrodynamic model in Cartesian coordinate Hirano:2000eu  and also in Bjorken coordinate Hirano:2001eu  and, moreover, took into account the picture of the early chemical freeze-out in this model Hirano:2002ds . Hydrodynamic model calculations at the RHIC energies suggest that the lifetime of the pure QGP phase is to be about 4 fm/$`c`$ and that the mixed phase exists up to 9–10 fm/$`c`$ Hirano:2002ds . The density of the system drops rapidly due to the longitudinal expansion.

In this Letter, we explore the dynamical effects on the energy loss of jets by taking into account the full 3D space-time evolution of a fluid. The hydrodynamic model is combined with the jets (Hydro + Jet model) which are calculated from the pQCD parton model and explicitly propagated in space-time with fluid elements.

We use full 3D hydrodynamics Hirano:2001eu  in which initial parameters are fixed for Au+Au collisions at $`\sqrt{s_{NN}}=130`$ GeV Hirano:2002ds . The transverse profile of the initial energy density is assumed to scale with the number of binary collisions KHHH . For 5% central collisions, we choose the impact parameter as $`b=2.4`$ fm and the maximum initial energy density at the initial time $`\tau _0=0.6`$ fm/$`c`$ as 33.7 GeV/fm3. These parameters lead us to reproduce transverse momentum spectra of charged hadrons Adler:2001yq  up to 1.0–1.5 GeV/$`c`$. It is found that the $`p_\mathrm{T}`$ slope is insensitive to the thermal freeze-out temperature $`T^{\mathrm{th}}`$ when we take into account the early chemical freeze-out Hirano:2002ds . We fix $`T^{\mathrm{th}}=140`$ MeV throughout this Letter.

We include hard partons using pQCD parton model,

$$\frac{d\sigma _{\mathrm{jet}}}{dp_\mathrm{T}^2dY_1dY_2}=K\underset{a,b}{}x_1x_2f_a(x_1,Q^2)f_b(x_2,Q^2)\frac{d\sigma _{ab}}{d\widehat{t}},$$
(1)
where $`Y_1`$ and $`Y_2`$ are the rapidities of the scattered partons and $`x_1`$ and $`x_2`$ are the fractions of momentum of the initial partons. The parton distribution functions $`f_a(x,Q^2)`$ are taken to be CTEQ5 leading order cteq5 . We use $`Q^2=p_\mathrm{T}^2`$ for the evaluation of parton distribution. The minimum momentum transfer $`p_{\text{T,min}}=2.0`$ GeV/$`c`$ is assumed. The summation runs over all parton species and relevant leading order QCD processes

$`q+q^{}`$ $``$ $`q+q^{},q+\overline{q}q^{}+\overline{q}^{},`$ (2)
$`q+\overline{q}`$ $``$ $`g+gq+gq+g,`$ (3)
$`g+g`$ $``$ $`q+\overline{q},g+gg+g.`$ (4)
are included in addition with the initial and final state radiation to simulate the emission of multiple soft gluons. Gaussian primordial transverse momentum $`k_\mathrm{T}`$ distribution with the width of $`k_\mathrm{T}^2=1`$ GeV2 is assigned to the shower initiator in the QCD hard $`22`$ processes. We use PYTHIA 6.2 pythia  to simulate each hard scattering in the actual calculation. A factor $`K`$ is used for the higher order corrections. We chose $`K=2`$ to fit the UA1 data of $`p\overline{p}`$ at $`\sqrt{s}=200`$ GeV ua1 . In order to convert hard partons into hadrons, we use an independent fragmentation model using PYTHIA after hydro simulations. We have checked that this hadronization model provides good agreement with the transverse spectra of charged hadrons above $`p_\mathrm{T}=1`$ GeV/$`c`$ in $`p\overline{p}`$ data ua1 .

The number of hard partons is assumed to scale with the number of hard scattering which is estimated by using Woods-Saxon nuclear density. This assumption is consistent with the peripheral Au+Au collision at $`\sqrt{s_{NN}}=130`$ GeV Fai . The space coordinate of a parton in longitudinal direction is taken to be $`\eta _\mathrm{s}=Y`$, where $`\eta _\mathrm{s}=(1/2)\mathrm{log}[(t+z)/(tz)]`$, and $`Y`$ is a momentum rapidity. The transverse coordinate of a parton is specified by the number of binary collision distribution for two Woods-Saxon distributions. Most hard partons are produced with the formation time of $`1/p_\mathrm{T}`$ assuming the uncertainty relation. Since the initial time of the hydrodynamical simulation is $`\tau _0=0.6`$ fm/$`c`$, partons are assumed to travel freely up to $`\tau _0`$. We neglect the nuclear shadowing effect in the present work for simplicity, because it is small at high transverse momentum Wang:2001gv ; jamal .

We assume that the form of energy loss is simply

$$\frac{dE}{dx}=ϵ\sigma \rho (\tau ,𝒓),$$
(5)
where $`ϵ`$ is an energy loss per scattering and $`\sigma `$ is a parton-parton cross section. From a hydrodynamic simulation, we obtain the space-time evolution of temperature $`T(\tau ,𝒓)`$. A parton density $`\rho (\tau ,𝒓)`$ in the QGP phase is calculated from $`T(\tau ,𝒓)`$. For the mixed phase at $`T=T_\mathrm{c}`$(=170 MeV), we put $`\rho (\tau ,𝒓)=f_{\mathrm{QGP}}(\tau ,𝒓)\rho (T_\mathrm{c})`$ where the fraction of the QGP phase in the mixed phase $`f_{\mathrm{QGP}}(\tau ,𝒓)`$ is calculated from energy density $`\epsilon (\tau ,𝒓)`$ and the maximum and minimum energy densities in the mixed phase $`\epsilon _{\mathrm{QGP}}`$ and $`\epsilon _{\mathrm{had}}`$ as

$$f_{\mathrm{QGP}}(\tau ,𝒓)=\frac{\epsilon (\tau ,𝒓)\epsilon _{\mathrm{had}}}{\epsilon _{\mathrm{QGP}}\epsilon _{\mathrm{had}}}.$$
(6)
We assume that the energy loss for a quark jet is half of the energy loss of a gluon wang2 . In the present work, we neglect the possible space-time variation of parton-parton cross section and regard the product $`ϵ\sigma `$ as an adjustable free parameter in the model. Feedback of the energy to fluid elements is ignored, because we use the relatively small values of energy loss parameter used in the work.

In Fig. 1, the transverse momentum distributions of neutral pions for central 10% Au + Au collisions at $`\sqrt{s_{NN}}=130`$ GeV from PHENIX experiment phenix\_pi0  are compared to the results of the hydro + jet model with different energy loss parameters $`ϵ\sigma =`$ 0, 0.06 and 0.3 GeV fm2 together with the hydro result by dotted line. In the calculation, we choose the impact parameter $`b=3.35`$ fm for 0–10% central events. The corresponding number of binary collisions is $`N_{\mathrm{binary}}=906`$. In order to connect the contribution from jets and a fluid smoothly, low $`p_\mathrm{T}`$ pions from jets are cut by a switch function $`\{1+\mathrm{tanh}[3(p_\mathrm{T}1.5\mathrm{GeV})]\}/2`$ v2jet1 . We reproduce $`p_\mathrm{T}`$ spectrum for $`\pi ^0`$ by choosing $`ϵ\sigma =0.06`$ GeV fm2, while the result which is not taken into account a energy loss overestimates the neutral pion spectrum. The hydro + jet result for $`dE/dx=0.2`$ GeV/fm which is obtained by neglecting the density dependence in Eq. (5) is also shown in Fig. 1. We note that the hydro + jet result for $`ϵ\sigma =0.3`$ GeV fm2 is reproduced by the constant energy loss of $`dE/dx=1.0`$ GeV/fm. As shown in Fig. 2, $`\pi ^0`$ spectra for peripheral (60-80%) Au + Au collisions phenix\_pi0  can be reproduced only by changing the impact parameter to $`b=12.1`$ fm ($`N_{\mathrm{binary}}=20`$) and leaving the other parameters. The hydro + jet result in peripheral collisions does not change when energy loss is not included, because the parton density and the volume of the QGP phase are small. On the other hand, the hydro + jet model overpredicts the data without energy loss in the central collisions. The deviation from data in the transverse momentum range below 1.5 GeV/$`c`$ in the peripheral collisions might be reasonable, because we do not expect that local thermalization is achieved in such a peripheral collision and that hydrodynamics is reliable.

The jet quenching is quantified by the ratio of the particle yield in A+A collisions to the one in p+p collisions scaled up by the number of binary collisions

$$R_{AA}=\frac{d^2N^{A+A}/dp_Td\eta }{N_{\mathrm{binary}}d^2N^{p+p}/dp_Td\eta }.$$
(7)
The transverse momentum dependence of $`R_{AA}`$ can provide information about the mechanism of jet quenching jeon . In Fig. 3, the result from the Hydro+Jet model with $`ϵ\sigma =0.06`$ GeV fm2 is compared to the PHENIX data phenix\_pi0 . Our result is almost flat in the range $`2<p_T<4`$ GeV/$`c`$.

Let us now turn to the study of the dynamical effects on the jet energy loss. The maximum and the average thermalized parton densities at $`\eta _\mathrm{s}=0`$ from hydrodynamic simulations are plotted as a function of proper time in the collision of impact parameter $`b=3.35`$ $`(12.1)`$ fm in upper (lower) panel of Fig. 4. As seen in the figures, thermalized parton density drops rapidly due to the strong longitudinal expansion of the system produced in heavy ion collisions. We expect that jets are likely to loss their energies in the time span of less than $`\tau 2`$ fm/$`c`$. In order to see when jets lose their energies, we plot in Fig. 5, the ratio of the numbers of jets $`N(\tau )/N(\tau _0)`$ at $`p_\mathrm{T}=5`$ GeV/$`c`$. It is found that about 90% of jet quenching occurs up to $`\tau 2`$ fm/$`c`$, although the QGP phase exists up to $`\tau 4`$ fm/$`c`$ and the mixed phase lasts until $`\tau 10`$ fm/$`c`$ within the hydro parameters used in this work. This holds true, even when we take the total density neglecting the existence of the hadrons in the mixed phase in Eq. (5) in order to get the maximum energy loss. Therefore, within our model, only QGP phase is responsible for the jet energy loss.

We showed in Fig. 1 the energy loss $`dE/dx=0.2`$ GeV/fm gives the same amount of energy loss as the energy loss parameter $`ϵ\sigma =0.06`$ GeV fm2. Let us estimate a space-time averaged energy loss parameter for central collision ($`b=3.35`$ fm). If we take a space-time average parton density $`\overline{\rho }(\tau )`$, $`\overline{\rho }(\tau \text{1.0 fm/}c)=11`$ fm-3, $`\overline{\rho }(\tau \text{2.0 fm/}c)=7.7`$ fm-3, and $`\overline{\rho }(\tau \text{3.0 fm/}c)=6.0`$ fm-3, then the space-time averaged energy loss yields $`dE/dx(\tau \text{1.0 fm/}c)0.66`$ GeV/fm, $`dE/dx(\tau \text{2.0 fm/}c)0.46`$ GeV/fm, and $`dE/dx(\tau \text{3.0 fm/}c)0.36`$ GeV/fm with the best fit parameter of $`ϵ\sigma =0.06`$ GeV fm2. Those values are close to that of the result calculated by the pQCD parton model Wang:2001gv , in which the effect of expansion is not included. The reason is that the energy loss parameter can be regarded as an average of roughly a short time of $`\tau 2`$ fm/$`c`$.

In summary, we proposed the full 3D hydrodynamical model combined with (mini-) jets, where jets are explicitly propagated in space-time with the hydro simulation. In particular, this model allows us to study the dynamical effects on the jet energy loss. We estimated the energy loss of partons by the fluid elements whose temperature is above the critical value $`T_c`$. It is found that energy loss of jets occurs in the pure QGP phase of $`\tau <2`$ fm/$`c`$ and that the contribution of energy loss in the mixed phase is negligible under the assumption of the energy loss formula Eq. (5). This indicates that suppression of hadronic high $`p_\mathrm{T}`$ spectra contain information about the early stage of partonic matter. However, it has been shown that the formula for the energy loss has non-trivial energy dependence Baier ; Levai  and it is found that energy loss is sensitive to the critical point DP01 . It is very important to take into account the coherent (Landau-Pomeranchuk-Migdal) effect in the calculation of radiation spectrum. A result obtained by using a different energy loss formula will be presented elsewhere.

A prediction at higher transverse momentum region is under progress. Elliptic flow parameter should be also studied within the model. It should be studied the charged particle spectra in order to have the unified understanding of the jet quenching mechanism in the medium.

In this study, we do not include the energy loss effects before hydrodynamical evolution. It might be interesting to study to what extent partons lose energy before thermalization using non-equilibrium models Nara:2001zu ; AR99 . Because parton density is maximum in this stage, energy loss by these partons should be important.

###### Acknowledgements.

We would like to thank T. Hatsuda, J. Jalilian-Marians, and A. Dumitru for useful comments. We would like to thank R. Venugopalan for reading the manuscript.


text:
# Azimuthal anisotropy of 𝐾_𝑆⁰ and Λ production at mid-rapidity from Au+Au collisions at √𝑠__{𝑁𝑁}=130 GeV

## 1 Introduction

Measurements of azimuthal anisotropies in the transverse momentum distributions of particles probe early stages of ultra-relativistic heavy-ion collisions Sorge (1999, 1997); Ollitrault (1992). After an initial geometric anisotropy is established in an off-center or non-central collision, rescattering in the overlapping region of the colliding nuclei amongst collision participants transfers the spatial anisotropy into an anisotropy in momentum space. The extent of the transformation depends on the initial conditions and the dynamical evolution of the collision. As a result, anisotropy measurements for nucleus-nucleus collisions at RHIC energies may increase our understanding of the processes governing the evolution of the collision system and in particular, may provide information about an early partonic stage in the evolution of the system Sorge (1999); Zhang et al. (1999); Huovinen et al. (2001); Gyulassy et al. (2001); Teaney et al. (2001); Lin and Ko (2002).

For the purpose of studying azimuthal anisotropies it is advantageous to write the triple differential distribution of particles in the form of a Fourier series

$$E\frac{d^3N}{d^3p}=\frac{1}{2\pi }\frac{d^2N}{p_tdp_tdy}(1+\underset{n=1}{\overset{\mathrm{}}{}}2v_ncos(n\varphi )),$$
(1)
where $`p_t`$ is the transverse momentum of the particle, $`y`$ is its rapadity and $`\varphi `$ denotes its azimuthal angle of emmision with respect to the true reaction plane angle111The reaction plane is the plane defined by the beam line and the line connecting the centers of the colliding nuclei. Voloshin and Zhang (1996); Poskanzer and Voloshin (1998). The harmonic coefficients $`v_n`$ are anisotropy parameters and the second coefficient $`v_2`$ is called the elliptic flow parameter. Recent experimental results from RHIC Ackermann et al. (2001); Adler et al. (2001a); Lacey (2002); Park et al. (2002) include measurements of $`v_2`$ as a function of collision centrality and $`p_t`$ for charged particles up to a $`p_t`$ of about 2.0 GeV/c, and for identified $`\pi ^\pm `$, $`K^\pm `$ and $`p(\overline{p})`$ up to a $`p_t`$ of about 0.8 GeV/c.

We report on the first measurement of the azimuthal anisotropy parameter $`v_2`$, as a function of $`p_t`$ and collision centrality, for the strange particles $`K_S^0`$ and $`\mathrm{\Lambda }`$ from minimum bias Au + Au collisions at $`\sqrt{s_{_{NN}}}=130`$ GeV. Our measurements of $`v_2`$ for identified particles using the Solenoidal Tracker At RHIC (STAR) are the first to extend beyond the $`p_t`$ range where particles are identified by their specific energy loss (dE/dx) in the gas of the Time Projection Chamber (TPC), and up to $`p_t3.0`$ GeV/c. Previously $`v_2`$ in this higher $`p_t`$ range had only been measured for unidentified charged particles Snellings (2002).

## 2 Analysis

The STAR detector Ackermann et al. (1999), due to its azimuthal symmetry and large acceptance, is ideally suited for measuring elliptic flow. For collisions in its center, the STAR TPC measures the tracks of charged particles in the pseudo-rapidity range $`|\eta |<1.8`$ with $`2\pi `$ azimuthal coverage. A scintillator barrel, the Central Trigger Barrel (CTB), surrounding the TPC that measures the charged particle multiplicity within $`|\eta |<1`$ was used for a central trigger. Two Zero-Degree Calorimeters Adler et al. (2001b) at both ends of the TPC in coincidence provided a minimum bias trigger.

The masses and kinematic properties of both $`K_S^0\pi ^++\pi ^{}`$ and $`\mathrm{\Lambda }(\overline{\mathrm{\Lambda }})p+\pi ^{}(\overline{p}+\pi ^+)`$, are reconstructed via their decay topologies in the TPC Wieman et al. (1997); Betts et al. (1997); Klein et al. (1996). Figure 1 shows the invariant mass distributions for a $`\pi ^+\pi ^{}`$ mass hypothesis and a $`p\pi ^{}(\overline{p}\pi ^+)`$ hypothesis. The background is dominated by combinatoric counts and the observed masses, $`496\pm 4`$ MeV/$`c^2`$ and $`1116\pm 7`$ MeV/$`c^2`$, are consistent with values listed in the PDG Groom et al. (2000) for $`K_S^0`$ and $`\mathrm{\Lambda }`$ respectively. The $`K_S^0`$ and $`\mathrm{\Lambda }`$ particles used in the $`v_2`$ analysis are from the kinematic region of $`0.20p_t3`$ GeV/c and $`|y|1.0`$. To compensate for limited statistics, $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ are summed together. To reduce the combinatoric background, for $`K_S^0`$, pion-like tracks are required to have a distance of closest approach dca $`>1.0`$ cm, while for $`\mathrm{\Lambda }`$, the pion-like tracks have a dca $`>1.5`$ cm and the proton-like tracks have a dca $`>0.8`$ cm. Tracks are determined to be either proton-like or pion-like based on their energy loss (dE/dx) in the TPC gas. The yield from the enhancement in the invariant mass peaks in each $`\varphi `$, $`p_t`$ bin is used to evaluate $`v_2=cos(2\varphi )`$ as a function of $`p_t`$. This method enables us to measure identified particle flow beyond the $`p_t`$ range where dE/dx particle identification fails.

The event plane, an experimental estimator of the true reaction plane Poskanzer and Voloshin (1998), is calculated from the azimuthal distribution of tracks using cuts similar to those used in reference Ackermann et al. (2001). To avoid auto-correlations, only tracks excluded from the neutral vertex reconstruction are used in the event plane calculation. The observed $`v_2`$ is corrected to account for the imperfect event plane resolution estimated using the method of subevents described previously Poskanzer and Voloshin (1998). The maximum resolution correction factor for the $`K_S^0`$ and $`\mathrm{\Lambda }`$ analysis is found to be $`0.681\pm 0.004`$ and $`0.582\pm 0.007`$ respectively and is reached in the centrality corresponding to 25 – 35% of the measured cross section, where the relative multiplicity distribution is used to estimate the event centrality as in reference Ackermann et al. (2001).

For this analysis, three sources contribute to systematic errors in the measured anisotropy parameters: (1) particle identification; (2) background subtraction; (3) non-reaction plane related correlations contributing to $`v_2`$ such as resonance decays or Coulomb and Bose-Einstein correlations Dinh et al. (2000); Borghini et al. (2000). The first two sources are estimated by examining the variation in $`v_2`$ after changing several track, event and neutral vertex cuts and are found to contribute an error of less than $`\pm `$ 0.005 to $`v_2`$. A previous study used the correlation of event plane angles from subevents to estimate the magnitude of non-reaction plane related correlations Adler et al. (2001a). That analysis showed that these effects which always act to increase the measured value of $`v_2`$ above its true value typically contribute a systematic error to $`v_2`$ of -0.005, but that the magnitude is larger in the more peripheral events where the error increases to about -0.035 for the centrality corresponding to 58 – 85% of the measured cross section.

## 3 Results

The centrality dependence of $`v_2`$ as a function of transverse momentum calculated from 201 thousand minimum bias and 180 thousand central events is shown in figure 2. The two particles show a similar $`p_t`$ dependence in the respective centralities with more flow in the more peripheral collisions. This is similar to observations made previously at the same energy Adler et al. (2001a) where the agreement with hydrodynamic calculations in the lower $`p_t`$ region was interpreted as evidence for early local thermal equilibrium in all but the most peripheral events (45 – 85% of the measured cross section).

In figure 3 (left) we plot $`v_2(p_t)`$ for $`K_S^0`$ and $`\mathrm{\Lambda }`$ from minimum bias collisions with results from hydrodynamic model calculations Huovinen et al. (2001) and $`v_2(p_t)`$ for negatively charged particles Snellings (2002). The $`K_S^0`$ results are in agreement with the $`v_2`$ for $`K^\pm `$ in the $`p_t`$ range they share ($`300p_t700`$ MeV/c) Adler et al. (2001a). We observe that $`v_2`$ for both strange particles increases as a function of $`p_t`$ similar to the hydrodynamic model prediction, up to about 1.5 GeV/c. In the higher $`p_t`$ region however ($`p_t2`$ GeV/c), the values of $`v_2`$ seem to be saturated. It has been suggested Gyulassy et al. (2001) that the shape and height of $`v_2`$ above 2 – 3 GeV/c is related to energy loss in an early, high parton-density stage of the evolution.

The $`p_t`$ integrated $`v_2`$ from minimum bias collisions for $`K_S^0`$, $`\mathrm{\Lambda }`$ and negatively charged particles are shown in figure 3 (right). The integrated values of $`v_2`$ are calculated by parameterizing the yield with the inverse slope parameter of exponential fits to the $`K_S^0`$ or $`\mathrm{\Lambda }`$ transverse mass distributions and are dominated by the region near the particles mean $`p_t`$. The relatively larger $`v_2`$ of $`\mathrm{\Lambda }`$ reflects the higher mean $`p_t`$ of the $`\mathrm{\Lambda }`$ compared to the $`K_S^0`$. Hydrodynamic model calculations Huovinen et al. (2001), shown as a gray-band and central line, are, within errors, in agreement with this result. The width of the gray-band indicates the uncertainty of the model calculation, mostly due to the choice of the freeze-out conditions. The increase of $`v_2`$ with particle mass in figure 3 points to a significant commonality in velocities between particles of different masses that is perhaps, established early in the collision. The nature of the particles during this process however, whether parton or hadron, and the degree of thermalization remains unclear.

## 4 Summary

We have reported the first measurement of $`v_2`$ for $`K_S^0`$ and $`\mathrm{\Lambda }`$ from Au + Au collisions at $`\sqrt{s_{_{NN}}}=130`$ GeV and the first measurement, at this energy, of $`v_2`$ above $`p_t0.8`$ GeV/c for any identified particle. For both particles more flow is seen in mid-central than in central collisions. The integrated values of $`v_2`$ show a mass dependence consistent with the development of a common velocity, a feature of hydrodynamic models, where local thermalization is assumed. In $`v_2(p_t)`$ however, we see a strong $`p_t`$ dependence only up to $`p_t`$ 2.0 GeV/c where $`v_2`$ seems to saturate, suggesting that the hydrodynamic picture is incomplete for particle production above $`p_t=`$ 1.5 – 2.0 GeV/c.

We thank P. Huovinen for hydrodynamic model calculations and the RHIC Operations Group at Brookhaven National Laboratory for their support and for providing collisions for the experiment. This work was vsupported by the Division of Nuclear Physics and the Division of High Energy Physics of the Office Science of the U.S. Department of Energy, the United States National Science Foundation, the Bundesministerium für Bildung und Forschung of Germany, the Institut National de la Physique Nucleaire et de la Physique des Particules of France, the United Kingdom Engineering and Physical Sciences Research Council, and the Russian Ministry of Science and Technology.


text:
# Relaxation Phenomena in Supercomputer Job Arrivals

## Abstract

We show that the distribution of supercomputer job submission interarrival times can be understood as a relaxation process. The process of deciding when to submit a job involves a complicated set of interactions between the users themselves, the queuing algorithm, the supercomputer, and a hierarchy of other decision makers. This is analogous to the hierarchically constrained dynamics found in glassy relaxation modelled by a stretched exponential. Empirical supercomputer log data shows that the tails of the distributions are well fit by a stretched exponential.

Today’s supercomputers have thousands of processors and perform sophisticated simulations on a wide variety of problems in material science, structural, and thermal dynamics. Supercomputers are an integral and enabling component in the complex system of Big Science. Among the most powerful supercomputers are those from the Advanced Simulation and Computing Initiative (ASCI)ASC . These machines were built for specific purposes to primarily serve a small group of users who end up dominating the cycles on the machine.

Supercomputers represent the largest single computing resources in the world and they must perform over a staggering range of conditions spanning small interactive jobs to very large jobs, both in terms of the number of processors involved (in the thousands) and for long time periods(on the order of a day or more for a single run). Similar to other complex systems, the workflow of jobs through a supercomputer system is a dynamic and complicated cycle of phases involving submission, dispatch, running, analysis, and resubmission. Often the “output” of a phase depends critically on one or more of the other phases. For example, the submission of a particular job at a particular time by a particular user depends on the time the user has to spend setting up the next run and the previous runs the user has to analyze. These in turn depend upon when it finished running on the machine, which depend upon when it was dispatched, which depend upon the prioritization constraints imposed by the facility managers via the queuing system. On top of these conditions is the laboratory hierarchy who approve projects and above that the governmental funding agencies and finally the elected officials who fund the facilities.

In 1854 Kohlrausch adapted Weber’s famous elasticity equation to explain the residual charge in a Leyden jar as a function of time and discovered the stretched exponential distributionKohlrausch (1854), namely, that the decay time probability of the relaxation process is given by

$$\varphi (t)=e^{(t/\tau )^\alpha }.$$
(1)
Since then his equation has found application in numerous relaxation processes of complex systems in nature including colloids, polymers, glasses, and more recently, radio emission from galaxies, earthquakes, oilfield reserve sizes as well as man-made phenomena such as certain market price variations and numbers of citationsLaherrere and Sornette (1998). In this paper we will show how job arrivals at a supercomputer can be mapped to a hierarchical relaxation process and therefore to Kohlrausch’s result.

Heavy-tailed distributions, defined here as those that drop off more slowly than an exponential, including the stretched exponential and power laws, have been reported in a number of manmade phenomena, specifically computer systems. Some examples of heavy tail distributions in computer systems include: computer networks both in terms of their connectivityWillinger et al. (2002) and their traffic patternsWillinger et al. (1996), file systemsGribble et al. (1998), video trafficBeran et al. (1995), software cachesVoldman et al. (1983), and the job size distributions on a single processorHarchol-Barter and Downey (1997). Ultimately, these computer systems are driven by some form of human activity interacting with algorithms hardcoded in the hardware or programmed into the software.

Heavy-tail distributions have important implications for both physical and manmade systems. In particular, heavy tails indicate a significant probability of very large events. In the case of earthquakes it means a meaningful chance for very large and damaging events. In the case of supercomputers it means the possibility that the machine may become overloaded for significant periods of time even if the average turnaround time is moderate. Significantly, the confluence of many large jobs impinging on a supercomputer as a consequence of heavy-tailed distributions both in job size and interarrival time can have serious consequences on the timeliness of the important work done at these facilities. Thus it is important to these facilities that the implications of these heavy-tails be characterized so that they may be taken into account in the design of queuing algorithms and in funding decisions for new hardware.

Despite the work in networks and single processor systems, little is known about the scaling behavior in the largest supercomputers. Given log data from thousands of jobs over a period of several months we can examine these issues quantitatively.

The complicated set of conditions for determining which job gets submitted and when is exactly the sort of conjunctive, i.e., multiplicative, process that is described by the stretched exponential distributionFrisch and Sornette (1997). More specifically, we can think of all the hierarchy of agents interacting in getting a job submitted to be in a discrete set of $`N`$ pseudospins arrayed in different levels for each agent class. In the case of Big Science the hierarchy is something like user $``$ project $``$ group $``$ facility $``$ laboratory $``$ government agency $``$ executive and legislative entities.

The relaxation function, $`\varphi `$, the probability of the system being in a state at time $`t`$ is given by

$$\varphi (N,t)=1/N\underset{n=0}{\overset{N}{}}S_i(0)S_i(t)$$
(2)
where $`S_i(t)`$ is the state of the $`i^{th}`$ pseudospin at time $`t`$ and $`N`$ is the number of levels. In terms of an ensemble of relaxation times we have

$$\varphi (N,t)=\underset{n=0}{\overset{N}{}}w_n\mathrm{exp}(t/\tau _n)$$
(3)
where $`w_n`$ is the relative number of pseudospins for level $`n`$. Following the arguments in Palmer et al. (1984) and Klafter and Shlesinger (1986), only $`\mu _nN_n`$ actually contribute to the decision at the $`n^{th}`$-level of the hierarchy. Under this scenario the $`\mu _n`$ spins in the $`n^{th}`$-level are free to change only when spins in level $`n1`$ have relaxed into one of their $`2^{\mu _{n1}}`$ possible states. If we ignore intralevel correlations then

$$\tau _{n+1}=2^{\mu _n}\tau _n.$$
(4)
Defining $`\stackrel{~}{\mu }_k=\mu _k\mathrm{ln}2`$ then

$$\tau _{n+1}=\tau _0\mathrm{exp}\left(\underset{k=0}{\overset{n}{}}\stackrel{~}{\mu }_k\right).$$
(5)
For a glassy relaxation, $`\mu _n`$ should decrease rapidly enough to make Eq.(5) converge. One such condition is given by $`\mu _n=\mu _0/n`$. For the job submission hierarchy we are considering, $`\mu _0`$ is about $`1001000sec`$, the average time between job submissions, the next level might be weekly meetings, a factor of ten thousand. For the highest levels in the hierarchy, the decisions are made on a much compressed scale relative to the next highest level. For example, the penultimate level meets quarterly and the highest level on a yearly scale, a difference of only a factor of $`4`$.

We also need to model the branching ratio between levels, or “span of control” in bureaucracy parlance. We model this as

$$w_n=w_0\lambda ^n.$$
(6)
Converting the sum to an integral we have

$$\varphi (t)=w_0_0^{\mathrm{}}\lambda ^n\mathrm{exp}(tn^{\stackrel{~}{\mu }_0}/\tau _0)𝑑n$$
(7)
This equation cannot be solved in closed form so by the method of steepest descent expanding around the point $`nt^\alpha `$ we finally obtain the desired result, Eq.(1), where $`\tau `$ defines a characteristic scale to the distribution (contrast this with a power law’s scale free behavior) and finally

$$\alpha =1/(1+\stackrel{~}{\mu }_0)$$
(8)
is a measure of the heaviness of the tail. The smaller the value of $`\alpha `$, the heavier the tail.

Another way of looking at a relaxation process is as a random walk in a fractal spaceJund et al. (2001). When the relaxation process is described by a stretched exponential this is seen as the signature of a fractal morphology of the configuration space at the current temperature of the system. In this view the complex morphology of the job submission landscape as the set of necessary steps needed for submission fall into place is what drives the system into its heavy-tailed relaxation. Table 1 shows an analogy between job events and a spin relaxation process.

To demonstrate that supercomputer job submissions can be understood as a stretched exponential relaxation process, we analyzed job logs from the ASCI supercomputers ASCI-BlueMountain (Los Alamos National Laboratory), and ASCI-BluePacific (Lawrence Livermore National Laboratory)ASC . Each lab has devised its own method for queuing jobs based in part on the historical political realities at each labClearwater and Kleban (2002a). The important thing to keep in mind is that the queuing algorithm through its prioritization and “backfilling,” (running jobs that are not first in the queue but can run now without slowing down the first job in the queue) acts to alter the order that jobs were submitted and thus when they will be dispatched, run, and finally analyzed, all affecting the next job to be submitted and thus the interarrival submission times.

For all the analysis shown below we tried to fit other distributions such as the exponential, lognormal, and power law functions, but none provided as good a fit and over such a long range as the stretched exponential. Qualitatively, the exponential fell off more rapidly than the data and the power law not fast enough. The lognormal fit well for smaller values, but did poorly at larger values, as one would expect from its functional form. Intuitively, we might expect the stretched exponential to be applicable and fill in this intermediate range with a moderately heavy tail and a characteristic scale.

Blue Mountain at Los Alamos has 5418 processors in its large partition.There were 8171 jobs in this sample taken over a period of 83 days. The distribution of interarrival times is shown in Fig. 1. The best fit (short dashes) to a stretched exponential is shown with a characteristic time of $`\tau =524`$s and $`\alpha =.57`$. Fits to lognormal and exponential are also shown. As can clearly be seen, only the stretched exponential is able to model the data well over its entire range.

The results from Blue Pacific at Livermore consisted of 57,430 jobs taken over a period of 63 days. Unlike Blue Mountain, Blue Pacific did not have any partitions and used about 1000 CPUs, although the full machine has more. The part of Blue Pacific we used was no longer fulfilling its primary mission to the ASCI program and is involved in more academic research.

Fig. 2 shows the cumulative distribution function for interarrival times. We have truncated our fit at 10,000 seconds because beyond that time the interarrival times are likely due to system issues and not user issues. For example, these events may correspond to outages in the machine or logging errors(about 10% of the log had bogus entries and were not used) that could anomalously effect the very long portion of the tail. For interarrival times up to 10,000s the parameters for the stretched exponential are a characteristic time of $`\tau =1655`$s and $`\alpha =0.61`$.

The parameters we found for the stretched exponential fit are shown in Table 2. It is interesting to note that for the interarrival time distribution, both LSF and DPCS have similar exponents for large jobs, $`\alpha _{LSF}=0.57`$ and $`\alpha _{DPCS}=0.61`$.

Our results have shown the applicability for the first time of the stretched exponential to describing distributions from supercomputer systems. Remarkably, the stretched exponential provided a good fit over the entire range of values for some of the cases we studied, spanning up to 8 orders of magnitude.

One interesting implication of the constrained hierarchical model we are using is the relationship between levels in the hierarchy, Eq.(4), which implies that relaxation (or response times in our case) take much longer as one gets farther from those doing the actual work. Interpreting the job submission process as a relaxation phenomenon where “barriers” to the decision to submit a certain sized job at a certain time must be overcome, the exponent $`\alpha `$ may be understood as related to the number of levels in a hierarchy that underlies the overall job submission processFrisch and Sornette (1997); Laherrere and Sornette (1998); Palmer et al. (1984).

We can test the convergence of 3 by using parameters derived from the supercomputers. Using the empirically determined value of $`\alpha _{LSF}=0.57`$, so that $`\stackrel{~}{\mu }_0=0.75`$ from Eq.(8). Values of $`\mu _n<1`$ correspond to weak constraints between levelsPalmer et al. (1984), not surprising in a scientific environment. Since we are plotting cumulative distribution functions, the $`w_n`$’s themselves don’t matter, but the span of control is critical. We choose $`\lambda =5`$ which is a typical span of control, $`\lambda `$ Eq. (7), in a high tech research lab. For $`\tau _0`$ we use the average time between job submissions, $`787`$s. We then plot $`\varphi (N,t)`$ for $`N=1`$ to $`4`$ in Fig. 3. The sum converges quickly and approximates that of the exact stretched exponential. From an organizational standpoint this tells us that no more than 4 or so levels in the hierarchy are having any effect on the time scale at which work gets done.

Both the supercomputers utilized in this research run under a “Fair Share”Henry (1984); Kay and Lauder (1988) algorithm (user priorities are decreased if they go over their “share”) so it will be interesting to see, when data becomes available, if another queuing algorithm, such as NQS (essentially first-in first-out) at Sandia has a similar characteristic exponent for job sizes and job interarrival times.

The characteristic scale implied by the stretched exponential distribution may prompt another look at some computer phenomena previously thought to exhibit scale-free behavior. It also tells us that the deviations from a power law are a fundamental part of the phenomenaLaherrere and Sornette (1998). After all, as big as these supercomputers are, they are still finite and their operators have put in additional constraints as well to satisfy administrative requirements, i.e., political realities. Together these constraints act to define a characteristic size of the distribution as well as the heaviness of the tail. For example, the size of jobs measured in terms of number of processors and run time was also found to be well modelled by a stretched exponential.Clearwater and Kleban (2002b)

In conclusion, we have shown that the interarrival time of jobs are not exponential nor do they posses pure power-law tails, but are somewhere in between and can be well fit by stretched exponentials over a large and important part of their range. These are indicative of finite scaling behaviors and have implications for the ultimate performance of these facilities because they relate to the frequency, and therefore the turnaround of big jobs that are the bread and butter of the ASCI supercomputers.

This paper would not have been possible with the expertise of many people who provided us with detailed knowledge of algorithms, data formats, as well as providing the job log data. The authors gratefully acknowledge the assistance of Stephany Boucher, Charles Hales, Michael Hannah, Steve Humphreys, Moe Jette, Wilbur Johnson, Tom Klingner, Jerry Melendez, Amy Pezzoni, Randall Rheinheimer, Phil Salazar, Bob Wood, and Andy Yoo. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy under Contract DE-AC04-94AL85000.