text:
# Particle production in central Pb+Pb collisions at 158 𝐴 GeV/𝑐

## I INTRODUCTION

The goal of heavy-ion collision research at ultra-relativistic energies is to understand the fundamental properties of nuclear matter under extreme conditions of high particle and energy density. When these densities and temperatures in heavy-ion collisions are sufficiently high, quarks should no longer be confined in hadrons but instead move freely through the entire volume. The colliding system expands and cools, driven by the large internal pressure generated by the strong interaction of the constituents. The quarks eventually become confined again, forming hadrons that can be detected by experiments. The hadron yields and momentum distributions provide information about the temperature, density, and dynamics of the later stages of the collision. From this inferences about the conditions at earlier stages can be made. QM (1)

The multiplicity of the produced hadrons ranges from several hundred to a few thousand in relativistic heavy-ion collisions at $`\sqrt{s_{NN}}`$20 GeV, and their rapidity density cannot be explained by simple superposition of $`p`$+$`p`$ collisions, as final state rescattering may play an important role. The transverse momentum distributions of the hadrons display a Boltzmann-like shape, inspiring a thermodynamical interpretation of the momentum distributions.

In this report, we present invariant cross sections of $`\pi ^\pm `$, $`K^\pm `$, $`p`$, and $`\overline{p}`$, as a function of transverse mass near mid-rapidity in central 158 A GeV/$`c`$ Pb+Pb collisions. NA44 has previously reported the particle yields for pions, kaons, and protons in $`p`$+A to S+A collisions Bearden et al. (1998); Bøggild et al. (1999), kaon production in Pb+Pb Bearden et al. (1999), and antiproton and antideuteron production in Pb+Pb Bearden et al. (2000). Because of the excellent particle identification and systematic measurements from small to large systems, NA44 is well-suited to study the thermal and chemical properties of the hadronic system at freeze-out. The dependence of the thermal parameters on the colliding system size will also be discussed.

## II NA44 EXPERIMENT

The NA44 focusing spectrometer is designed to measure one and two particle momentum distributions of charged hadrons near mid-rapidity Bøggild et al. (1993). A schematic view of the spectrometer is shown in Fig. 1. The target is a lead disk of 2 mm thickness (3.4% interaction probability for a lead beam). A Cherenkov gas beam counter (CX) Maeda et al. (1994) and a beam veto counter (CX-veto) are placed upstream of the target. CX selects single beam particles (rejecting events with multiple simultaneous beam particles) up to a rate of $`2\times 10^6`$ ions/second while the CX-veto signals the presence of any beam halo. A scintillation multiplicity counter (T0) is 10 mm downstream of the target, and produces a pulse height proportional to the number of charged particles impinging upon it. T0 covers 20% of the azimuthal angle and has a pseudo-rapidity range of 0.6$``$$`\eta `$$``$3.3. This contains mid-rapidity ($`\text{y}_{mid}=`$2.9) for the Pb+Pb collision. T0 is used to select collision centrality at the trigger level.

Two dipole magnets (D1 and D2) analyze the particle momenta. Three super conducting quadrupole magnets (Q1, Q2, and Q3) focus the particles. The PID/tracking detectors measure hit positions after the magnetic analysis, and consist of strip chambers and scintillation hodoscopes, along with two threshold gas Cherenkov counters, C1 and C2. Three hodoscopes, H2, H3, and H4 consist of 60, 50, and 60 scintillation counters, respectively. Only H2 and H3 are used in this analysis. The hodoscopes provide position information for tracking, with vertical hit information determined from photomultiplier tubes on each end of every counter. The hodoscopes also measure time-of-flight for hadron identification Kobayashi and Sugitate (1990). One pad chamber (PC) and two strip chambers (SC1 and SC2) provide 0.3 mm position resolution for track reconstruction. C1 and C2 contain Freon-12 at 2.7 atm and N2 at 1.3 atm respectively, and identify particles at the trigger level. The trigger requires a valid beam, a high multiplicity detected in T0 (central collision), at least one track in each hodoscope, and the absence of an electron (all runs) or pion (for $`K`$ and $`p`$ runs).

The timing signal from CX is used for the time-of-flight (TOF) start with a resolution of 35 ps. The average of the top and bottom PMT of the hodoscope is used for the stop signal. The overall TOF resolution is approximately 100 ps and 80 ps for H2 and H3, respectively.

Two magnetic field settings of the spectrometer are used, yielding nominal momentum of 4 GeV/$`c`$ ($`3.2p5.2`$ GeV/$`c`$) and 8 GeV/$`c`$ ($`6.0p9.8`$ GeV/$`c`$). Two angular settings, 44 mrad and 129 mrad are used to give a large range in transverse momentum ($`p_T`$). Combining these settings yields transverse momenta from $`p_T=0`$ to 0.9 GeV/$`c`$ (0 to 1.6 GeV/$`c`$) for 4 GeV/$`c`$ (8 GeV/$`c`$) setting. Figure 2 shows the NA44 acceptance in $`p_T`$ and rapidity for pions, kaons, and protons. The NA44 acceptance includes mid-rapidity ($`\text{y}_{mid}`$=2.9) for Pb+Pb collisions at 158 GeV/$`c`$ per nucleon.

The invariant cross section is determined by the following analysis procedures: track and momentum reconstruction, particle identification of tracks, a selection of high multiplicity events, an acceptance correction via a Monte-Carlo simulation, normalization by summing of beam scalers, and multiplicative correction factors for cuts.

The hit positions on PC, H2, H3, SC1, and SC2 were used to reconstruct tracks by fitting them with a straight line using a $`\chi ^2`$-minimization method. As is clear from Fig. 1, the tracking detectors lie outside the magnetic field. The momentum resolution $`\delta p/p`$ is 0.2% for all particle species.

Particle identification uses a combination of pulse heights in C1 and C2 with time-of-flight from CX to the hodoscopes. The Cherenkov threshold momenta of C1(C2) for $`e`$, $`\pi `$, $`K`$, and $`p`$ are 0.014, 2.5, 8.9, and 16.9 (0.020, 5.3, 18.3, and 34.8) GeV/$`c`$, respectively. In the 4 GeV/$`c`$ spectrometer setting, pions are separated from $`K`$ and $`p`$ by C1 and C2, and then $`K`$/$`p`$ separation is done by TOF. In the 8 GeV/$`c`$ setting, $`e`$ and $`\pi `$ are vetoed by C2, and then $`K`$/$`p`$ separation is done by C1. The $`e`$/$`\pi `$ separation is not able to be done by C1 and C2 in the 8GeV/$`c`$ setting, however the contamination from electrons is negligible for $`p_T`$$`>`$0.35 GeV/$`c`$. The statistics for each particle are shown in Table 1.

Because it has a longer flight path in addition to its better timing resolution, H3 rather than H2 provides the primary TOF information. Figure 3 shows an example of particle identification in the 4 GeV/$`c`$ setting. The left-hand side shows a scatter plot of pulse heights in C1 and C2. Events containing an electron fire C2 – these events are rejected. The events selected by the box labeled (a) and (b) on the left-hand side of figure are shown on plots of TOF as a function of momentum on the right-hand side of the figure. Particles between the curves are selected. The contamination for each particle species is less than 3% in both momentum settings.

We trigger on central collisions using the interaction counter, T0, which measures the total amount of light produced by charged particles crossing its scintillators. Figure 4 shows the pulse height distribution of the T0 for valid beam events (1) and central triggers (2). The shaded region is used for the analysis. The width of the valid beam distribution comes from the production of delta electrons. Since T0 covers only 20% of the azimuthal angle, it is not possible to count all low multiplicity Pb+Pb interactions and we cannot measure the Pb+Pb interaction cross section. Therefore we scale the measured $`p`$+Pb cross section Groom et al. (2000) by a geometry factor: this fraction is denoted by $`\lambda `$. The value of $`\lambda `$ is obtained by scaling the interaction probability for $`p`$+Pb collisions (the nuclear interaction length) to the interaction probability for Pb+Pb collisions by a geometry factor of $`(208^{\frac{1}{3}}+208^{\frac{1}{3}})^2`$/$`(1^{\frac{1}{3}}+208^{\frac{1}{3}})^2`$, i.e., $`\lambda `$=2.93. The nuclear interaction length of Pb is 194 g/cm2 and density of Pb is 11.35 g/cm3. Therefore, the interaction probability of $`p`$+Pb of 2 mm Pb target is $`0.2/(194/11.35)=0.0117`$. The total Pb interaction probability is 3.4%$`(=\lambda \times 0.0117)`$ for 2 mm Pb target.

We then find that the shaded region in Fig. 4 corresponds to $`3.7\pm 0.1`$% (statistical only) of the Pb+Pb interaction cross section. This region corresponds to a region of impact parameter $`b`$$`<`$3 (RMS$``$1) fm, estimated by applying the T0 acceptance cut to the multiplicity distribution from RQMD (v2.3) Sorge (1995). We note that the T0 covers only 20% of the azimuthal angle, therefore the event centrality selected by T0 has a weak sensitivity to multiplicity, that is, to impact parameter.

## III ANALYSIS

In order to construct the invariant cross section, the effect of the acceptance is estimated by a Monte-Carlo (MC) simulation based on TURTLE TUR (1974), which passes generated tracks through the detectors. The acceptance, reconstruction efficiencies, momentum resolution, and loss of pions and kaons due to in-flight-decays are estimated by comparing generated tracks with reconstructed tracks. The momentum spectrum and rapidity distribution for each particle species in Pb+Pb collisions is taken from RQMD (v2.3). The Monte-Carlo events were analyzed with the tracking programs exactly as applied to the physics events. The corrections are given by the ratio of the generated yield to the reconstructed yield as a function of $`m_T`$ inside the rapidity acceptance. The absolute cross section as a function of $`m_T`$ in a rapidity range $`\mathrm{\Delta }y`$ is given by simplify this equation.

$`{\displaystyle \frac{1}{\sigma }}{\displaystyle \frac{Ed^3\sigma }{dp^3}}|_{\mathrm{\Delta }y}`$ $`=`$ $`{\displaystyle \frac{1}{N_{\text{event}}}}{\displaystyle \frac{1}{2\pi }}{\displaystyle \frac{dN}{m_Tdm_T\mathrm{\Delta }y}}`$
$`=`$ $`{\displaystyle \frac{1}{N_{\text{event}}}}{\displaystyle \frac{dn}{dm_T\mathrm{\Delta }y}}\times {\displaystyle \frac{\frac{1}{2\pi }\frac{1}{m_T}\frac{d\stackrel{˘}{n}_{mc}}{dm_T\mathrm{\Delta }y}}{\frac{d\stackrel{˘}{n}_r}{dm_T\mathrm{\Delta }y}}}`$
where $`n`$ is the number of measured tracks, $`N_{\text{event}}`$ is the number of collisions, $`\stackrel{˘}{n}_{mc}`$ is the number of particles as input of the MC simulation and $`\stackrel{˘}{n}_r`$ is the number of tracks reconstructed in the simulation.

Secondary hadronic interactions of protons and antiprotons in the spectrometer material does not distort the shape of the momentum distributions but reduces the yieldBearden et al. (1998). The correction factors are estimated to be 1.11 for protons and 1.17 for antiprotons, respectively.

The corrected momentum distribution is then normalized using the number of beam particles, the interaction probability, the fraction of interactions satisfying the trigger, the measured live time of the data-acquisition system, and losses due to particle identification cuts.

Separation of protons and kaons from pions is done via Cherenkov veto trigger, but in the process some protons and kaons are rejected in events that also include a pion. The veto factor (i.e., the correction for this loss) for each particle and momentum setting is estimated by comparing the number of the particles with no Cherenkov veto to the number of particles with Cherenkov veto, using real data without the Cherenkov veto trigger. Since the veto factor depends on relative particle ratios in spectrometer, it must be independently estimated for each momentum and angle setting.

The transverse mass distribution is obtained from the two angle settings (44 mrad and 129 mrad). The data at 129 mrad have smaller (typically one half to one third) systematic errors in the absolute normalization than those at 44 mrad. Therefore, in this analysis, the data at 44 mrad were scaled to match those at 129 mrad in the $`m_T`$ region where they overlap. The distributions from both settings were fit by exponentials with the same inverse slope parameter. The dN/dy is the sum of the normalized cross section in each measured $`m_T`$ bin within the acceptance plus an extrapolation that uses the fitted exponential distribution beyond the measurement region. This extrapolation contributes less than 1% of the yield for pions but up to 7% for protons with the 4 GeV/$`c`$ setting.

## IV SYSTEMATIC ERRORS

The inverse slope parameters are extracted by fitting the experimental distribution with a single $`m_T`$ exponential function (see Eq. (1) in the section VI). Systematic errors on the inverse slope parameters are estimated by dividing the data into two $`p_T`$ ranges for the fits. The systematic errors are shown with the inverse slope parameters in Table 2.

The absolute normalization procedure results in significant uncertainties in dN/dy. Contributions include uncertainties on the following; the data-acquisition live time ($`5\%`$), the stability of the T0 pulse height used to select the event fraction ($`2\%`$), the Cherenkov veto factor ($`<`$$`5\%`$), the correction for particle loss due to quality cuts ($`<`$$`5\%`$), and the extrapolation beyond the measurement region ($`<`$$`17\%`$). The last three vary with particle species and spectrometer setting. We assume that the errors are uncorrelated and so add them in quadrature to get the total systematic error. The errors for each particle species and setting are shown with results in Table 2.

## V FEED-DOWN FROM WEAK DECAYS

NA44 has no vertex detectors in the target region, so decay protons from some fraction of the short-lived particles such as strange baryons will be detected as protons in the spectrometer. This contamination of protons and antiprotons is not negligible.

The same Monte-Carlo used for acceptance corrections was also used to estimate the weak decay yields from $`\mathrm{\Lambda }`$, $`\mathrm{\Sigma }^+`$ and $`\mathrm{\Sigma }^0`$. Since the $`\mathrm{\Sigma }^0`$ lifetime is short, its yield was included in the $`\mathrm{\Lambda }`$ yield. The dN/dy for various particles near mid-rapidity in central collisions was taken from RQMD(v2.3).

The correction factor as a function of $`m_T`$ is

$$C_p(m_T)=\frac{N_{p_{\text{orig}}}(m_T)}{N_{p_{\text{orig}}}(m_T)+N_{p_\mathrm{\Lambda }}(m_T)+N_{p_\mathrm{\Sigma }}(m_T)}$$

where $`N_{p_{\text{orig}}}`$ is the number of the original protons and $`N_{p_\mathrm{\Lambda }}`$ and $`N_{p_\mathrm{\Sigma }}`$ represent the protons from $`\mathrm{\Lambda }`$ decays, and $`\mathrm{\Sigma }^+`$ decays that track through the spectrometer. Figure 5 shows the correction factors as a function of $`m_T`$mass. Since these depend on the ratio of $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }^+`$ to proton, the production ratios were varied. The circles are the correction factor estimated from $`\mathrm{\Lambda }/p`$ and $`\mathrm{\Sigma }/p`$ ratios of RQMD, while the lines assume yield ratios larger or smaller by a factor 1.5.

The correction procedure contributes to the systematic uncertainties on the inverse slopes of proton and antiproton by about 7% for the 4 GeV/$`c`$ setting and about 3% for the 8 GeV/$`c`$ setting. The uncertainties due to feed-down correction on dN/dy are approximately 15% for protons and 19% for antiprotons in both 4 GeV/$`c`$ and 8 GeV/$`c`$ settings.

## VI RESULTS

The invariant cross sections measured by the NA44 spectrometer in Pb+Pb collisions are shown in Fig. 6. They are plotted as a function of $`m_T`$ in the specified rapidity regions. The $`m_T`$ distribution is generally well described at these energies by

$$\frac{1}{\sigma _{trig}}\frac{Ed^3\sigma }{dp^3}=Ae^{(m_Tm)/T},$$
(1)
where $`A`$ is a constant and $`T`$ is the inverse slope parameter. Table 2 lists the inverse slope parameters, the fit region, and the values of dN/dy.

The slope parameters are consistent with previous NA44 results (Ref. Bearden et al. (1997)) within errors. The proton dN/dy is consistent with the values in Ref. Bearden et al. (1996). The statistics in our previous publication are $`1/3`$ to $`3/4`$ of the statistics used in this paper. We note that the NA44 collaboration published a paper on pion and kaon spectra focusing on “strangeness enhancement” Bearden et al. (1999). This analysis used a different definition of the event fraction, based on the silicon multiplicity counter. The results are consistent within errors.

The NA49 collaboration has reported dN/dy distributions of $`\pi ^+,\pi ^{},K^+,K^{},p`$, and $`\overline{p}`$ for the 5$`\%`$ most central Pb+Pb reactions Appelshäuser et al. (1999); Siklér (1999); Cooper (1999); Veres (1999). These values are slightly higher than our data (top 3.7$`\%`$ event fraction). However, since they trigger on the energy of forward going spectators their selection of low impact parameter is better than ours. The top 5$`\%`$ NA49 data corresponds to mean impact parameter $`b`$=2.2 fm Cooper (1999). On the other hand, the NA44 data corresponds to $`b`$$``$3 fm. Comparing similar impact parameter regions, the values of dN/dy are consistent with the NA49 results within the errors.

We have obtained the invariant cross sections for pions, kaons, and protons in Pb+Pb collisions using data analysis procedures compatible with these used in Refs. Bearden et al. (1998); Bøggild et al. (1999). Consequently, we can study the systematics of particle yields at CERN-SPS energy for many collision systems (450 GeV/$`c`$ $`p`$+A, 200 A GeV/$`c`$ S+A, and 158 A GeV/$`c`$ Pb+Pb).

Figure 7 shows the inverse slope parameters for pions, kaons, and protons versus the $`dN_{ch}/dy`$ near central rapidity. All of the inverse slopes increase with $`dN_{ch}/dy`$ but this effect is more rapid for protons than for kaons and pions. The slopes of particles and antiparticles are equal except for protons have slightly larger inverse slopes than antiprotons for $`p`$+A and S+A collisions. While the uncertainties on the slope parameters for pions make a definitive statement difficult, the pion slope parameters in the region $`m_T`$mass$`<`$0.6 GeV/$`c^2`$ in Pb+Pb are similar to, or less than, the slope in $`p`$+A and S+A collision systems. Perhaps, the lower values in the low $`p_T`$ region for Pb+Pb reflect a larger enhancement at low $`p_T`$.

The kaon, proton and antiproton inverse slope parameters increase with collision system size slowly, but more rapidly than the pion inverse slope parameters. The proton inverse slopes are higher at mid-rapidity (y=2.3–2.9) than at more backward rapidities (y=1.9–2.3) in $`p`$+A and S+Pb collisions.

Figure 8 shows the collision system dependence of charged pion, charged kaon, proton, and antiproton yields, and ratios of negative to positive charged particle. The pion rapidity density in Pb+Pb collisions increases by a factor of $`71\pm 10`$ compared to $`p`$+Be collisions, and by a factor $`5.1\pm 0.6`$ compared to S+S collisions. dN/dy for $`\pi ^{}`$ and $`\pi ^+`$ is compatible, within the errors for all collision systems.

The mid-rapidity dN/dy of $`K^+`$ ($`K^{}`$) in Pb+Pb collisions increases by a factor of $`133\pm 20`$ ($`105\pm 10`$) compared to $`p`$+Be collisions and $`6.4\pm 0.9`$ ($`6.2\pm 0.8`$) compared to S+S collisions. The system size dependence is stronger than for pions, i.e., the $`K/\pi `$ ratio increases from $`p`$+Be to Pb+Pb.

Proton dN/dy increases with charged multiplicity at mid-rapidity more rapidly than antiprotons. In heavy-ion collisions, the $`\overline{p}/p`$ ratio decreases with system size, and falls more rapidly at mid-rapidity than at backward rapidity. This is exactly what is expected from increased stopping power, and therefore larger net baryon density, with larger colliding system size.

## VII DISCUSSION

In heavy-ion collisions, we use the notions of chemical freeze-out and thermal freeze-out to describe the moment when the inelastic and elastic interactions stop. In this section, we will describe the parameters extracted from the measured spectra and yields.

### VII.1 Transverse momentum distributions

The transverse momentum distributions of hadrons have been measured by several experimental groups at BNL-AGS and CERN-SPS including NA44 Bearden et al. (1997). The $`m_T`$ distributions are studied from the viewpoint of collective flow, with local thermal equilibrium at mid-rapidity, based on Refs. Schnedermann et al. (1993); Cörgo and Lörstad (1996). The NA44 collaboration has reported a thermal freeze-out temperature $`T_{th}140`$ MeV Bearden et al. (1997) in Pb+Pb collisions at $`\sqrt{s_{NN}}`$ 18 GeV. Recently two-pion interferometry results have been combined with single particle spectra to measure the transverse flow and thermal freeze-out temperature Rayford Nix (1998); Appelshäuser et al. (1999). These studies find $`T_{th}120`$ MeV Appelshäuser et al. (1999) and 100 MeV Rayford Nix (1998). Evaluating those reports equally, one concludes that the thermal freeze-out temperature is in the range 100–140 MeV.

### VII.2 Particle yields

If we assume a thermodynamic description is valid for these systems, then at chemical freeze-out, the density of different particles can be characterized by macroscopic parameters: the chemical freeze-out temperature and chemical potentials. Here, we adopt a chemical freeze-out model based on Ref. Sollfrank et al. (1999), using the grand canonical ensemble. We assume no difference between $`u`$ and $`d`$ quark chemical potentials. Therefore, the hadron gas is described by a chemical freeze-out temperature ($`T_{ch}`$), light ($`u`$ and $`d`$) quark potential ($`\mu _q`$), strange quark potential ($`\mu _s`$), and strangeness saturation factor ($`\gamma _s`$). The density of a particle $`i`$ in the hadron gas is given by:

$$\rho _i=\gamma _{s}^{}{}_{}{}^{s+\overline{s}_i}\frac{g_i}{2\pi ^2}T_{ch}^{}{}_{}{}^{3}\left(\frac{m_i}{T_{ch}}\right)^2K_2(m_i/T_{ch})\lambda _{q}^{}{}_{}{}^{Q_i}\lambda _{s}^{}{}_{}{}^{s_i}$$
(2)
where $`m_i`$ is the mass of the hadron $`i`$, $`g_i`$ is the number of spin-isospin degree-of-freedom, $`K_2`$ is the second-order modified Bessel function and,

$$\lambda _q=\mathrm{exp}(\mu _q/T_{ch}),\lambda _s=\mathrm{exp}(\mu _s/T_{ch}).$$

Note that we use the Boltzmann approximation for all hadrons except pions, where the Bose distribution is applied. The potential $`\mu _q`$ is for u/d/$`\overline{\text{u}}`$/$`\overline{\text{d}}`$ quarks, and $`\mu _s`$ is for s/$`\overline{\text{s}}`$ quarks. $`Q_i`$ and $`s_i`$ correspond to the net number of valence u/d quarks ($`Q_i=u\overline{u}+d\overline{d}_i`$), and s quark ($`s_i=s\overline{s}_i`$) of particle species $`i`$, respectively. The factor $`\gamma _s`$ ($`0\gamma _s1`$) is introduced to allow of possibility of incomplete chemical equilibration for strange particles Rafelski (1991). $`\gamma _s`$=1 signifies full strangeness chemical equilibrium. The power factor of $`\gamma _s`$, $`s+\overline{s}_i`$, is total number of s and $`\overline{\text{s}}`$ quarks in the particle. Ideally the system would live long enough that strange quarks would become fully equilibrated with the light quarks and $`\gamma _s`$ would be unity. However, calculations in Refs. Becattini et al. (2001); Becattini (1998); Letessier and Rafelski (1999) support $`\gamma _s`$$`<`$1. Therefore, we adopt $`\gamma _s`$ as a free parameter. The particle densities are computed by Eq. (2) for the hadron gas including known resonances up to mass of 1.7 GeV/$`c^2`$, and the effect of broad resonance mass widths is taken into account Sollfrank et al. (1999). The resonance decay to lower mass hadrons after chemical freeze-out is also taken into account. Therefore, the final particle densities for the pions, kaons, and protons are given by the sum of the direct production and the decays of the resonances. Here, we apply the above procedure for the particle yields observed near mid-rapidity by NA44 for $`p`$+A, S+A, and Pb+Pb collisions Bearden et al. (1998); Bøggild et al. (1999). The model is fit to the data using MINUIT James (1984).

Figure 9 and Table 3 shows the collision system dependence of the parameters. The chemical freeze-out temperatures ($`T_{ch}`$) are approximately 140 MeV in $`p`$+Be through Pb+Pb collisions: Within errors, there is no dependence on system size. Comparing to other results, the $`T_{ch}`$ is consistent with Ref. Letessier and Rafelski (1999) ($`T_{ch}`$=133–144 MeV) and smaller than Ref. Becattini et al. (2001) ($`T_{ch}`$=158$`\pm `$3 MeV) and Ref. Braun-Munzinger et al. (1999) ($`T_{ch}`$=168$`\pm `$2.4 MeV). As reported in Refs. Sollfrank et al. (1999); Kaneta and Xu (2001) there is a rapidity dependence of the chemical freeze-out parameters. Those calculations show that a lower value of the temperature is obtained when only mid-rapidity data are used. The limited number of particle species used in this analysis may also influence the value of $`T_{ch}`$.

The values of $`\mu _q`$ increase with the size of the collision system. Since the chemical freeze-out temperature is independent of collision systems ($`p`$+A and A+A), the system in $`p`$+A collisions has smaller $`\mu _q/T_{ch}`$, that is, smaller baryon density than A+A collisions. $`\mu _s`$ is 10$``$20 MeV for $`p`$+A through Pb+Pb. $`\mu _s`$ should be zero in a equilibrated quark gluon plasma and on the phase boundary of a hadron gas Rafelski (1991); Lee et al. (1988); Asprouli and Panagioutou (1995). Consequently, the small value we obtain may suggest that chemical freeze-out takes place near the phase transition.

Finally we find that $`\gamma _s`$ increases logarithmically with $`dN_ch/dy`$ from $`0.49\pm 0.06`$ in $`p`$+Be collisions to $`0.74\pm 0.05`$ in Pb+Pb. Thus, we seem to be approaching full chemical equilibrium for strange particles in heavy-ion collisions. The use of the grand canonical ensemble for analysis of strange particle production in $`p`$+A collisions does not take into account phase space effects arising from production of $`s`$ and $`\overline{s}`$ pairs. Consequently, $`\gamma _s`$ is expected to be lower than in heavy-ion collisions where more strange quark pairs are produced Redlich et al. (2001); Redlich (2002).

## VIII CONCLUSION

The NA44 experiment has measured single particle distributions for charged pions, kaons, and protons as functions of transverse mass ($`m_T`$) near mid-rapidity in 158 A GeV/$`c`$ Pb+Pb collisions. We studied particle ratios near mid-rapidity from the viewpoint of chemical freeze-out for $`p`$+A (A = Be, S, and Pb), S+A (A = S, Pb), and Pb+Pb collisions. We find that chemical freeze-out occurs at a temperature of $``$140 MeV for all collisions. The values of the light quark chemical potential ($`\mu _q`$) are approximately 30 MeV in $`p`$+A and 54 to 63 MeV in S+S to Pb+Pb collisions. The trend of $`\mu _q`$ shows that the baryon stopping power increases with collision system size. The values of the strange quark chemical potential are about 10 to 20 MeV. The strangeness saturation factor is 0.61 to 0.75 in heavy-ion collisions and 0.46 to 0.58 in $`p`$+A, indicating that $`p`$+A collisions are further away from strangeness equilibration than A+A collisions.

###### Acknowledgements.

The NA44 collaboration wishes to thank the staff of the CERN PS-SPS accelerator complex for their excellent work. We are also grateful for support given by the Danish Natural Science Research Council; the Japanese Society for the Promotion of Science; the Ministry of Education, Science and Culture, Japan; the Swedish Science Research Council; the Austrian Fond für Förderung der Wissenschaftlichen Forschung; the National Science Foundation, and the US Department of Energy.


format: text

text:
# The Optical Afterglow of the Gamma-Ray Burst GRB 0112111footnote 11footnote 1Based on observations collected at the OGLE 1.3m telescope, the Magellan 6.5m Walter Baade telescope, the FLWO 1.2m telescope, and the Vatican Advanced Technology Telescope.

## 1 Introduction

The gamma-ray burst (GRB) GRB 011211 (catalog ) was detected in the constellation Crater by the BeppoSAX satellite at 19:09:21 UT on 2001 Dec. 11. The burst was a shallow, long event with two peaks and a total duration of approximately 270 s making it the longest event that has been localized by BeppoSAX. The temporal profiles of the event were similar in gamma- and $`X`$-rays. BeppoSAX measured a gamma-ray fluence of $`5\times 10^6`$ erg cm-2 between 40 keV and 700 keV (Frontera et al., 2002). Approximately ten hours after the burst occurred Grav et al. (2001) identified an optical source within the BeppoSAX error circle that was not present in the Digital Sky Survey 2. Bloom & Berger (2001) and Jensen et al. (2001a) reported that this source was fading and Soszyński et al. (2001) estimated that it had a power law decay with a slope of $`\alpha =0.93\pm 0.06`$. Fruchter et al. (2001) found, and Gladders et al. (2001) confirmed, a redshift of $`z=2.14`$ based on several absorption lines in the spectrum of the optical afterglow (OA). Burud et al. (2001) identified a host galaxy with $`R_{\mathrm{host}}=25.0\pm 0.3`$ and found that the OA is offset $`0\stackrel{}{\mathrm{.}}5`$ southeast from the center of this object.

In this paper we present photometry and spectroscopy of the OA of GRB 011211 (catalog ) taken between $`0.5`$ and 2.7 days after the burst occurred. Our data suggests that there were rapid fluctuations in the optical flux approximately 12 hours after the burst occurred. We adopt a cosmology with a Hubble parameter of $`H_0=65`$ km s-1 Mpc-1, a matter density of $`\mathrm{\Omega }_m=0.3`$, and a cosmological constant of $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. For this cosmology a redshift of $`z=2.140`$ corresponds to a luminosity distance of 18.18 Gpc and a distance modulus of 46.30. One arcsecond corresponds to 28.06 comoving kpc, or 8.94 proper kpc. The lookback time is 11.28 Gyr.

## 2 The Photometric Data

The OA for GRB 011211 (catalog ) is located at R.A. $`=`$ 11:15:17.98, Dec. $`=`$ $``$21:56:56.2 (J2000) (Grav et al., 2001), which corresponds to Galactic coordinates of $`(b^{\mathrm{II}},\mathrm{}^{\mathrm{II}})=(+35\stackrel{}{\mathrm{.}}695,333\stackrel{}{\mathrm{.}}491)`$. The reddening maps of Schlegel et al. (1998) give a Galactic reddening of $`E_{BV}=0.043\pm 0.020`$ mag in this direction. The corresponding Galactic extinctions are $`A_V=0.142`$, $`A_R=0.114`$, and $`I=0.083`$. We estimate the extinction in the host galaxy in § 7.

Our photometry data was collected at three telescopes: the Optical Gravitational Lensing Experiment (OGLE) 1.3m telescope at the Las Campanas Observatory, the Magellan 6.5m Walter Baade telescope at Las Campanas, and the Fred Lawrence Whipple Observatory (FLWO) 1.2m telescope. The OGLE 1.3m telescope was equipped with the 8k $`\times `$ 8k OGLE-III CCD, which a had a scale of $`0\stackrel{}{\mathrm{.}}26`$ per pixel, a gain setting of 1.0 e-/ADU, and a read-out noise of 5 e-. The Magellan images were taken with the LDSS-2 imaging spectrograph in its imaging mode. It had a scale of $`0\stackrel{}{\mathrm{.}}378`$ per pixel, a gain of 1.0 e-/ADU, and a read-out noise of 7 e-. The FLWO images were taken with the “4Shooter” CCD mosaic (Szentgyorgyi et al., 2002). The pixel scale was $`0\stackrel{}{\mathrm{.}}67`$ (binned $`\times 2`$), the gain was 4.0 e-/ADU, and the read-out noise was 10 e-. A log of the observations is given in Table 1. Fig. 1 shows the field containing the OA.

The field stars around GRB 011211 (catalog ) were calibrated on 2002 Jan. 12 (UT) using images taken with the Vatican Advanced Technology Telescope. Landolt standard fields (Landolt, 1992) were observed in $`UBVRI`$ filters throughout the night under photometric conditions. We derived airmass and color corrections and computed magnitudes for eight stars near GRB 011211 (catalog ). These magnitudes and colors are given in Table 2. The magnitudes for the USNO-A2 star U0675\_11427359 (catalog ) (our star S2) are consistent with those of Henden (2002).

We used DAOPhot II (Stetson, 1987) and AllStar (Stetson & Harris, 1988) to perform point-spread function fitting photometry on the calibrated field stars, and on the OA, on each image. A zero-point offset was calculated for each filter on each night and applied to the observed magnitudes. Our calibrated photometry for the OA is presented in Table 1. The two FLWO images in each band were averaged to obtain the final images. The OA was not detected in any of our 2001 Dec. 14 FLWO images. Therefore, we estimated upper limits to the brightness of the OA on 2001 Dec. 14.5 from the limiting magnitudes of these images, which are listed in Table 1.

The OA is consistent with a point source in all of our images. We aligned and co-added all of our $`R`$-band images to obtain a deep image (shown in Fig. 1) with a limiting magnitude of $`R_{lim}23.5\pm 0.2`$. There is no evidence for a host galaxy in this image.

## 3 The Spectral Data

A spectrum of the OA was obtained with the Magellan 6.5m Walter Baade telescope using the LDSS-2 imaging spectrograph on 2001 Dec. 13.3 UT, approximately 1.5 days after the burst. The slit width was $`1\mathrm{}`$ and the total exposure time was $`4\times 620`$ s. The resolution is $`12`$ Å. The spectra were reduced in the standard manner and the wavelength calibration was done using a He–Ne arc. Several weak metal lines were identified to determine the redshift of the system. These are listed in Table 3. The mean redshift is $`z=2.140\pm 0.001`$ (standard error \[se\]). Table 3 also lists the observed and rest-frame equivalent widths for each line.

The measured absorption line widths for GRB 011121 (catalog ) are quite similar to those observed in GRB 000301C (catalog ) (Jensen et al., 2001b) but weaker than the exceptional case of GRB 000926 (catalog ) (Castro et al., 2001; Fynbo et al., 2002). Both of these GRBs have redshifts similar to that of GRB 011211 (catalog ). We also observe an unidentified, broad feature at approximately 4600 Å. This may be an additional absorption system along the line of sight, but the signal-to-noise ratio of our spectrum is not high enough to confirm this hypothesis. A plot of the spectrum, with our line identifications, is given in Fig. 2.

## 4 The Optical Light Curve

The $`VRI`$ light curves of the OA of GRB 011211 (catalog ) are shown in Fig. 3. The flux from the host galaxy (Burud et al., 2001) has not been subtracted. The pixel sizes in our images are comparable to the reported offset between the OA and the host ($`0\stackrel{}{\mathrm{.}}5`$), so our quoted magnitudes include the flux from both the OA and the host galaxy. However, the host is faint enough that it will contribute only $`1`$% of the flux when the OA has $`R20`$ and $`10`$% of the flux when the OA has $`R22`$. This corresponds to systematic errors of only 0.01–0.10 mag. Therefore, we believe that the flux from the host galaxy does not significantly affect our results.

We fit all of the data in Table 1 up to 1.6 days after the burst with a power law of the form

$$M_j=K_j2.5\mathrm{log}_{10}(C_jt^\alpha )$$
(1)
where $`M_j`$ is the magnitude for filter $`j`$, $`K_j`$ is the photometric zero point for filter $`j`$ taken from Fukugita et al. (1995), $`t`$ is the time in days since the burst occurred, $`\alpha `$ is the slope, and $`C_j`$ is a normalization constant representing the flux from the OA in filter $`j`$ one day after the burst. The best-fitting power law to our data has $`\alpha =0.83\pm 0.04`$ with $`\chi ^2/\mathrm{DOF}=19.222/14`$ and a root-mean-square (RMS) residual of 0.09 mag.

Our power-law fit suggests that the OA should have $`R=21.9\pm 0.1`$ at $`t=2.72`$ days. However, we find that the OA has $`R>22.4`$ at this time, $`5\sigma `$ fainter than the predicted magnitude. This non-detection suggests that a break occurred in the light curve between 1.52 and 2.72 days after the burst. If the break occurred at more than 1.52 days then the late-time decay has a slope of $`\alpha _21.4`$, which corresponds to the decay becomeing steeper by $`\mathrm{\Delta }\alpha =\alpha _2\alpha 0.6`$. The FLWO $`V`$-band data are consistent with this, but the FLWO $`I`$-band data are consistent with there being no break in the light curve. Burud et al. (2001) found $`R=24.8\pm 0.3`$ for the OA ten days after the burst. If we fit a broken power law to our $`R`$-band data and their point we find a late-time slope of $`\alpha _2=1.70\pm 0.26`$ and a break time of $`t_b=1.80\pm 0.62`$ days, which are consistent with the estimates of the late-time slope and break time made using only our data.

Fig. 3 shows that a single power law with a slope of $`0.83\pm 0.04`$ is a good fit to the $`VRI`$ data up to 1.52 days after the burst. Using Eq. 1 we find colors of $`VR=0.27\pm 0.13`$, $`VI=0.70\pm 0.13`$, and $`RI=0.43\pm 0.13`$ before the break. We stress that these colors are corrected for Galactic extinction, but not for extinction in the host galaxy of GRB 011211 (catalog ). In § refSECTION:extinct we show that the extinction in the host along the line of sight to the burst is probably small compared to the Galactic extinction in that direction. Assuming a spectrum of the form $`f_\nu \nu ^\beta `$ we find a weighted mean dereddened spectral slope from these colors of $`\beta =0.61\pm 0.15`$ (se). This is similar to the dereddened spectral slope found for many other GRBs.

## 5 The Energy in the Burst

The BeppoSAX fluence is $`5\times 10^6`$ erg cm-2 in the 40–700 keV band (Frontera et al., 2002). Applying a cosmological $`k`$ correction (Bloom et al., 2001) this fluence corresponds to an isotropic equivalent energy of $`E_{\mathrm{iso}}=7\times 10^{52}`$ erg between 20 keV and 2000 keV. Our $`k`$ correction assumes that the gamma-ray spectrum has the same form as the Band et al. (1993) spectrum with $`\alpha =1`$, $`\beta =2`$, and $`E_0=150`$ keV. These are the cannonical values for a GRB’s gamma-ray spectrum but Band et al. (1993) found that they are not universal values and that individual bursts can have very different spectral shapes. We estimated the uncertainty in our $`k`$ correction to be $`20`$% based on changing the spectral shape parameters by factors of two.

The opening angle, $`\theta _0`$, of a GRB jet is related to the time of the break in the light curve (Rhoads, 1999; Sari et al., 1999). Frail et al. (2001) cast this relation as

$$\theta _0=0.057\left(\frac{t_b}{1\mathrm{day}}\right)^{3/8}\left(\frac{1+z}{2}\right)^{3/8}\left(\frac{E_{\mathrm{iso}}}{10^{53}\mathrm{erg}}\right)^{1/8}\left(\frac{\eta _\gamma }{0.2}\right)^{1/8}\left(\frac{n}{0.1\mathrm{cm}^3}\right)^{1/8}\mathrm{rad}$$
(2)
where $`\eta _\gamma `$ is the efficiency of converting energy in the ejecta into gamma rays and $`n`$ is the circumburst particle density. Frail et al. (2000) and Panaitescu & Kumar (2001) find $`0.001n3`$ cm-3 in the vicinity of five GRBs with a median number density of $`n=0.1`$ cm-3. Therefore, we adopt this as the circumburst number density. We will assume, as did Frail et al. (2001), that $`\eta _\gamma =0.2`$. Eq. 2 gives $`\theta _0=3\stackrel{}{\mathrm{.}}4`$$`4\stackrel{}{\mathrm{.}}2`$ for break times of $`t_b=1.5`$$`2.7`$ days. From this we estimate that GRB 011211 (catalog )’s total beamed energy in gamma rays, after correcting for the beam geometry, was $`E_\gamma 1.2`$$`1.9\times 10^{50}`$ erg. This is only $`2\sigma `$ smaller than the “standard” total beamed energy in gamma rays of $`5\times 10^{50}`$ erg (Frail et al., 2001; Piran et al., 2001; Panaitescu & Kumar, 2002). This agreement suggests that our assumptions of $`\eta _\gamma 0.2`$ and $`n0.1`$ cm-3 are reasonable.

In order for $`E_\gamma `$ to equal the “standard” value the opening angle of the beam needs to be $`7\mathrm{°}`$. This corresponds to a break time of $`10`$ days, which is inconsistent with the observed brightness of the OA on 14 Dec. 2001. An opening angle of $`7\mathrm{°}`$ can be made consistent with our estimate of the break time by increasing the circumburst particle density to $`n5`$–30 cm-3. This is at the high end of the range of particle densities found by Frail et al. (2000) and Panaitescu & Kumar (2001) for several GRBs, but not inconsistent with there results. Therefore, we believe that the environment of GRB 011211 (catalog ) was similar to the environments of other GRBs.

## 6 The Ambient Medium Near the Burst

For a collimated outflow into an ambient medium with a number density distribution of the form $`n(r)r^\delta `$ (Panaitescu et al., 1998; Mészáros et al., 1998) the power-law index, $`\delta `$, is related to the observed magnitude of the break in the light curve by $`\delta =(4\mathrm{\Delta }\alpha 3)/(\mathrm{\Delta }\alpha 1)`$. We find $`\mathrm{\Delta }\alpha 0.6`$, which corresponds to $`\delta 1.5`$. This is not consistent with expansion into an ambient medium that is dominated by a pre-existing stellar wind ($`\delta =2`$), but might be consistent with expansion into a homogeneous medium ($`\delta =0`$). If the ambient medium is, on average, uniform then the magnitude of the break in the light curve should be $`\mathrm{\Delta }\alpha =3/4`$, which implies a late-time slope of $`\alpha _2=1.6`$. This is consistent with the late-time slope that we find ($`1.4`$) from our non-detection of the OA on 2001 Dec. 14. We note that it is also in good agreement with the late-time slope that we find ($`1.70\pm 0.26`$) using the data of Burud et al. (2001). Further evidence that the ambient medium is not dominated by a pre-existing stellar wind comes from the relationships between the time and spectral evolutions of the flux in a homogeneous medium (Sari et al., 1999) and in a wind (Chevalier & Li, 1999). Using the observed pre-break slope of the light curve a homogeneous medium predicts $`\beta =0.55\pm 0.03`$ for $`\delta =0`$ in the fast cooling regime, which is within $`0.5\sigma `$ of the observed value. In the slow cooling regime the predicted slope is $`\beta =0.22\pm 0.03`$, which is a worse agreement. A pre-existing wind predicts $`\beta =0.89\pm 0.03`$ for the fast cooling case and $`\beta =0.22\pm 0.03`$, for the slow cooling case. Both are worse fits than the fast-cooling homogeneous medium model (although all four cases are within $`2.5\sigma `$ of the observed value). Therefore, we believe that a homogeneous medium with fast-cooling electrons provides better agreement to the observations than the wind model does. The electron power law distribution index for this case is $`p=2.1\pm 0.1`$, which is consistent with what is seen in other bursts.

Fig. 4 suggests that there may be rapid variations in the $`R`$-band light $`0.5`$ days after the burst. To test this we subtracted the best-fitting power-law from the $`R`$-band magnitudes and computed the residuals. The $`R`$-band residuals have a weighted mean of $`0.02\pm 0.03`$ with an RMS of 0.10 mag. This says that we can reject the hypothesis that the residuals are consistent with random scatter about a power-law decay at the 95% confidence level. Therefore, we believe that the small-scale variations in the $`R`$-band light $`0.5`$ days after the burst may be real. We find no evidence for variability in the $`I`$-band data.

Garnavich et al. (2000) have postulated that the rapid variability seen in GRB 000301C (catalog ) was due to gravitational microlensing. However the rapid $`R`$-band variations seen in GRB 011211 (catalog ) are much smaller than those seen in the OA of GRB 000301C (catalog ). An alternate explanation is that there is small-scale structure in the local medium around the burst. Wang & Loeb (2000) find that linear density variations on spatial scales of $`1`$$`10^3`$ AU in the ambient medium that a GRB is expanding into can cause rapid fluctuations in the optical flux similar to those seen in Fig. 4. Their methodology predicts fluctuations of approximately 5%–10% over time periods of $`0.25`$–1.50 hours for GRB 011211 (catalog ). The observed scatter in the $`R`$-band magnitudes of several field stars is 0.06 mag and the RMS residual in the OA magnitudes is 0.10 mag. Therefore the RMS fractional variability in the $`R`$ band at about 0.5 days after the burst is $`8`$% on time-scales of $`1`$–2 hours. This is consistent with the predicted variability.

Using Table 3 of Wang & Loeb (2000), and the values for the isotropic equivalent energy of the burst and the density of the circumburst medium that we find in Sect. 5, we can estimate the size of the density fluctuations near the burst and the radius of the shock front. For $`E_{\mathrm{iso}}=7\times 10^{52}`$ erg and $`n=0.1`$–10 cm-3 the observed $`R`$-band fluctuations correspond to typical density variations on spatial scales of $`40`$–125 AU. At 0.5 days after the burst the radius of the shock front is $`\mathrm{12\hspace{0.17em}000}`$–40 000 AU1. (The apparent faster-than-light motion of the shock front is an illusion caused by the highly relativistic velocity of the shock front.) Therefore, the density fluctuations are small compared to the region that has been swept out by the GRB and our finding that the circumburst environment can be represented by a homogeneous medium is, on average, valid. The density fluctuations are between $`\delta n/n0.3`$ (for $`n0.1`$ cm-3) and $`3`$ (for $`n10`$ cm-3). Smaller regions will have larger density fluctuations than larger regions will. These fluctuations are similar to those seen in the interstellar medium in our Galaxy (Diamond et al., 1989; Faison & Goss, 2001).

## 7 Extinction in the Host

We can estimate the amount of extinction in the host galaxy along the line of sight to GRB 011211 (catalog ) by comparing the observed spectral slope after correcting for Galactic reddening ($`\beta =0.61\pm 0.15`$) to the predicted intrinsic slope found in § 6 ($`\beta _0=0.55\pm 0.03`$). This difference corresponds to $`E_{BV}=0.02\pm 0.05`$ in the host galaxy, which yields an extinction is $`A_V=0.06\pm 0.15`$ in the observer’s frame. Using the extinction law of Cardelli et al. (1989) we find $`A_V0.03`$ in the rest-frame of the burst. This result depends somewhat on the details of the extinction law used, but it suggests that there is no significant extinction along the line of sight to GRB 011211 (catalog ) in its host galaxy.

## 8 Conclusions

We present early-time $`VRI`$ photometry and spectroscopy of the OA of GRB 011211 (catalog ) starting approximately 0.5 days after the burst. The spectrum contains several narrow metal lines that are consistent with a redshift of $`2.140\pm 0.001`$. There is an unidentified broad feature at $`4600`$ Å that may be an absorption system along the line of sight to the GRB.

The OA is red with a spectral slope between 5505 Å and 8060 Å of $`\beta =0.61\pm 0.15`$ after correcting for Galactic extinction, but ignoring extinction in the host. This corresponds to a color of $`VI=0.70\pm 0.13`$. The magnitude of the break in the optical decay is consistent with the burst expanding into an approximately homogeneous medium. If we assume that GRBs have a the standard energy suggested by Frail et al. (2001), Piran et al. (2001), and Panaitescu & Kumar (2002) then the ambient medium near the burst has a particle density of 0.1–10 cm-3. Comparing the observed color of the OA with predictions of the standard beaming model suggests that the rest-frame $`V`$-band extinction in the host galaxy is $`A_V0.03`$ mag.

GRB 011211 (catalog ) follows the same broad pattern established by other long-duration GRBs. The OA decays as a power law with a slope of $`\alpha =0.83\pm 0.04`$ for the first $`2`$ days, and there is evidence for a break occurring between 1.52 and 2.72 days after the burst. The slope of the light curve after the break is $`\alpha _21.4`$. We find evidence for variations of $`8`$% in the flux on time scales of 1–2 hours occurring approximately half a day after the burst. Interpreting this in the framework of Wang & Loeb (2000) suggests that there are density fluctuations with scales of approximately 40–125 AU within $`0.05`$–0.20 pc of the GRB’s progenitor. The discovery of rapid variations in the optical light from a GRB highlights the importance of continuous, high-precision observations of GRBs at early times.

We wish to thank the BeppoSAX team, Scott Barthelmy, and the GRB Coordinates Network (GCN) for rapidly providing precise GRB positions to the astronomical community. We also wish to thank Arne Henden for providing precision photometry of stars in GRB fields. The authors would like to thank the anonymous referee for their helpful comments on this paper. STH and PMG acknowledge support from the NASA LTSA grant NAG5–9364. DB has been supported by NSF grand AST–9979812. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA.


text:
# Giant hole and circular superflow in a fast rotating Bose-Einstein condensate

## Abstract

A fast rotating Bose-Einstein condensate confined in a quadratic-plus-quartic potential is found to dynamically generate a “giant vortex” that absorbs all phase singularities into a central low density hole, thereby sustaining a quasi-one-dimensional circular superflow at a supersonic speed.

Gaseous Bose-Einstein condensates (BECs) are an extremely versatile testing ground for superfluidity particularly when an externally driven rotation exists. Through a number of experimental Madison ; Abo ; Haljan ; Hodby  and theoretical studies Fetter1  on rotating BECs, the argument now extends to the fast rotating regime Ho ; Fetter2 ; Feder ; Fisher ; Lundh . For a harmonic trapping potential $`(1/2)m\omega _{}^2r^2`$, the centrifugal force prevents a BEC from rotating at frequency $`\mathrm{\Omega }`$ beyond the radial trapping frequency $`\omega _{}`$. However, a recent experimental development Kuga  should enable one to create a confinement potential tighter than harmonic, thus opening a possible method to explore the nature of fast rotating BECs.

In this article, we extend our previous studies Tsubota  to the dynamics of vortex lattice formation of a BEC in a quadratic-plus-quartic potential with $`\mathrm{\Omega }`$ greater than $`\omega _{}`$. Here, the stable vortex configuration, realized at a value of $`\mathrm{\Omega }`$ much greater than $`\omega _{}`$, is not triangular but rather a circular array around a central low density hole \[see Fig. 1(e) below\]. This hole is created where singly quantized vortices are close together but do not completely overlap. Surprisingly, as $`\mathrm{\Omega }`$ is further increased, all the vortices are absorbed into the central low density hole, around which persistent current circulates. This superflow is supersonic, which offers an avenue of research in superfluidity, because it has been difficult to obtain a circular superflow in superfluid helium Donnely . For an atomic-gas BEC, a circular superflow can occur in an ideal situation and the superfluid speed is easily controlled by changing $`\mathrm{\Omega }`$. This opens up possibilities of studying phase slippage Anderson  and the stability of a controlled superflow.

We consider a two-dimensional system subject to rotation $`𝛀=\mathrm{\Omega }\widehat{𝐳}`$ by assuming translation symmetry along the $`z`$ axis. In a frame rotating with frequency $`\mathrm{\Omega }`$ around the $`z`$ axis, the dynamics of a condensate “wave function” $`\mathrm{\Psi }`$ is described by the Gross-Pitaevskii equation

$$i\frac{\mathrm{\Psi }}{t}=\left(^2+V_{\mathrm{trap}}\mu +C|\mathrm{\Psi }|^2\mathrm{\Omega }\widehat{L}_z\right)\mathrm{\Psi },$$
(1)
where $`\mu `$ is the chemical potential and the wave function is normalized to unity: $`𝑑x𝑑y|\mathrm{\Psi }|^2=1`$. The units of energy, length, and time are given by the corresponding scales of the harmonic potential, that is, $`\mathrm{}\omega _{}`$, $`b_{}=\sqrt{\mathrm{}/2m\omega _{}}`$, and $`\omega _{}^1`$, respectively, where $`m`$ is the atomic mass and $`\omega _{}`$ the frequency of the harmonic potential. The strength of the mean-field interaction between atoms is characterized by $`C`$ which is proportional to the $`s`$ wave scattering length $`a_s`$ as $`C=8\pi N_{2\mathrm{D}}a_s`$, where $`N_{2\mathrm{D}}`$ is the particle number per unit length along the $`z`$ axis. The trapping potential is $`V_{\mathrm{trap}}(r)=r^2/4+kr^4/16`$, where $`r^2x^2+y^2`$ and $`k`$ is the strength of the quartic term. The quartic part of $`V_{\mathrm{trap}}`$ allows us to increase $`\mathrm{\Omega }`$ above $`\omega _{}`$. For a rotating condensate an effective potential in the radial direction is

$$V_{\mathrm{trap}}^{\mathrm{eff}}(r)=\frac{1}{4}\left[(1\mathrm{\Omega }^2)r^2+\frac{1}{4}kr^4\right].$$
(2)
This potential has a “Mexican hat” structure for $`\mathrm{\Omega }>1`$, which changes the dynamics of vortex lattice formation as shown below.

Now, we investigate the dynamics of vortex lattice formation numer . The stationary solution of Eq. (1) for a nonrotating trap is used as the initial state. The rotating drive is suddenly turned on by introducing a small anisotropy of the harmonic trap as $`(ϵ_xx^2+ϵ_yy^2)/4`$ Madison ; Tsubota . Because the transition from a nonvortex state to a vortex state requires energy dissipation, a phenomenological damping constant $`\gamma `$ is put in Eq. (1) by replacing $`i/t`$ with $`(i\gamma )/t`$ Tsubota ; Choi . Finally, the small trap anisotropy $`ϵ_{x,y}`$ is turned off adiabatically to obtain the axisymmetric stationary state.

Numerical simulations were done for several values of $`\mathrm{\Omega }`$ at fixed $`C=250`$ and $`k=1`$. Dynamics that are unique to the quadratic-plus-quartic potential occur when $`\mathrm{\Omega }>1`$. Figure 1 shows the time development of the angular momentum per atom $`\mathrm{}_z=𝑑x𝑑y\mathrm{\Psi }^{}\widehat{L}_z\mathrm{\Psi }`$, and the density and phase profiles for $`\mathrm{\Omega }=2.5`$. The ripples are excited on the surface \[Fig. 1(b)\], and some ripples penetrate into the condensate and develop into vortex cores \[Fig. 1(c)\]. Some penetrating vortices move toward the center, merging together to make a hole at the center \[Fig. 1(d)\]. During the period between that shown in Fig. 1(c) and that in Fig. 1(d), additional vortices entered the condensate successively, which gradually increased the angular momentum. The density and phase profiles of the final steady state are shown in Figs. 1(e) and 1(f), respectively. A hole appears at the center, around which some vortices form a circular array. Furthermore, near the center of the condensate, some phase defects come very close to each other without overlapping \[Fig. 1(f)\].

As $`\mathrm{\Omega }`$ is further increased, all vortices generated from the condensate surface are absorbed into the central low density hole as shown in Fig. 2(a). Note that this central hole is composed of 24 singly quantized vortices packed together \[Fig. 2(b)\]. This packing is possible for high rotation frequencies ($`\mathrm{\Omega }>1`$) because the centrifugal force decreases the condensate density near the center; thus packing vortices costs little energy. The minimum of Eq. (2) determines the radius of the ring condensate in Fig. 2 as $`R=\sqrt{2(\mathrm{\Omega }^21)/k}`$. Here we call a set of vortices such as that in Fig. 2 a “giant vortex” to indicate that a number of phase defects are contained in a single hole \[see Fig. 2(b)\]. Such a density depression with a concentration of vorticity occurs in rotating superfluid helium Marston .

Figures 1(f) and 2(b) show that the phase singularities do not completely overlap. The condensate inside the hole has a very low density and hence a long coherence length, so that the vortex cores have a nearly uniform density. Thus, only the locations of singularities matter and one can assume that the interaction energy between the singularity points is proportional to $`\rho \kappa ^2\mathrm{ln}d`$, where $`\rho `$ is the condensate density, $`\kappa `$ the quantum circulation, and $`d`$ the distance between defects. The logarithmic divergence at $`d0`$ prevents the phase singularities from merging together.

The appearance of such a density hole is relevant to recent studies on fast rotating BECs that used the Thomas-Fermi and Wigner-Seitz approximation Fetter2 ; Fisher . Fischer and Baym evaluated the vortex core size self-consistently by minimizing the Gross-Pitaevskii energy functional and using a variational wave function Fisher ; they showed that vortex cores never overlap, in agreement with our finding. They also showed that the giant vortex state is energetically favorable at large $`\mathrm{\Omega }`$. Although the Thomas-Fermi approximation accurately describes stationary solutions such as that in Fig. 1(e), it fails when the width of the ring condensate becomes thin such as that in Fig. 2. Even if $`\mathrm{\Omega }`$ is not large enough, a pinning potential might be used to stabilize the giant vortex Simula .

For sufficiently large $`\mathrm{\Omega }`$, the ring structure of Fig. 2(a) offers some insights into circular superflow, which has been studied for a torus Donnely . The superfluid velocity is given by the gradient of the phase of the wave function as $`𝐯_s=(\mathrm{}/m)\theta `$; the circulation along a closed loop is quantized as $`𝐯_s𝑑𝐥=\kappa n`$ with $`\kappa h/m`$ being the quantum circulation. The quantum number $`n`$ is equal to the number of branch cuts of the phase field along the loop or, equivalently, the number of phase defects in the low density hole. If there is sufficient dissipation, some vortices in the low density hole should spiral out of the condensate when $`\mathrm{\Omega }`$ is abruptly decreased from that of Fig. 2(a). The motion of these vortices across the circular superflow causes $`v_s`$ to decay. This is known as phase slippage Anderson , which is a characteristic feature of superfluidity.

Measurement of the angular momentum of the giant vortex allows the phase slippage to be observed well under control. For example, phase slippage was observed in helium superflow through microapertures Varoquaux , but it was difficult to control single phase slips and giant slips. Figure 3 shows the time development of the angular momentum and phase field resulting from a decrease in $`\mathrm{\Omega }`$ from $`\mathrm{\Omega }=3.5`$ ($`n=32`$) to $`\mathrm{\Omega }=3.2`$. Here we assumed a small anisotropy in the trapping frequencies to break the rotational symmetry. Some phase defects are released from the low density hole, jostling at the inner border of the ring condensate. We see that four vortices cross the superflow by phase slippage and then escape outside (only one escape is shown in Fig. 3), thus decreasing the angular momentum.

When the ring radius $`R`$ is sufficiently large, the circular flow becomes quasi-one-dimensional for the following reason. The radial wave function of the stationary giant vortex state is obtained by substituting $`\mathrm{\Psi }f(r)e^{in\theta }`$ into Eq. (1):

$$\left[\frac{^2}{r^2}\frac{1}{r}\frac{}{r}+\frac{n^2}{r^2}+V_{\mathrm{trap}}+C|f|^2n\mathrm{\Omega }\right]f=\mu f.$$
(3)
The excitation of the radial component is easily obtained by writing $`f+\delta f`$; the effective potential for $`\delta f`$ is then given by

$$V_{\mathrm{eff}}(r)=V_{\mathrm{trap}}+\frac{n^2}{r^2}+2C|f(r)|^2\mu n\mathrm{\Omega }.$$
(4)
Because $`n=\mathrm{\Omega }R^2/2`$ in the limit of a rigid body rotation Fetter2 ; Feder , the radial zero-point energy around the minimum of the potential $`V_{\mathrm{trap}}+n^2/r^2`$ becomes larger than the mean-field interaction strength $`C|f(r)|^2`$ for large $`R`$. Such a high excitation gap freezes the motion of the radial wave function while keeping the motion free along the $`\theta `$ direction. To investigate the stability of this one-dimensional superflow, we substitute $`\mathrm{\Psi }(r,\theta ,t)=f(r)\psi (\theta ,t)`$ into Eq. (1) to obtain

$$i\frac{\psi }{t}=\left(\frac{1}{R^2}\frac{^2}{\theta ^2}+E_r+CP|\psi |^2i\mathrm{\Omega }\frac{}{\theta }\right)\psi ,$$
(5)
where

$`E_r={\displaystyle \frac{1}{A}}{\displaystyle r𝑑rf^{}\left[\frac{^2}{r^2}\frac{1}{r}\frac{}{r}+V_{\mathrm{trap}}\mu \right]f},`$
$`R^2{\displaystyle \frac{1}{A}}{\displaystyle r𝑑r|f|^2r^2},P={\displaystyle \frac{1}{A}}{\displaystyle r𝑑r|f|^4}`$
with $`A=r𝑑r|f|^2`$. The wave function is linearized as $`\psi (\theta ,t)=e^{i\mu t}\{e^{in\theta }+\eta [e^{i\omega t}u(\theta )+e^{i\omega ^{}t}v^{}(\theta )]\}`$. The resulting eigenvalue problem with the ansatz $`u(\theta )=u_qe^{i(n+q)\theta }`$ and $`v(\theta )=v_qe^{i(n+q)\theta }`$ yields the dispersion relation

$$\omega _{q,\pm }=\frac{q}{R}\left\{\left(\frac{2n}{R}R\mathrm{\Omega }\right)\pm \sqrt{\left(\frac{q}{R}\right)^2+2CP}\right\},$$
(6)
where $`q/R`$ is the wave number of the excitation. The second term on the right-hand side of Eq. (6) gives the sound velocity $`c=\sqrt{2CP}`$. The superflow is stable when $`v_s<c`$ according to the Landau criterion. The values of the parameters and the form of $`f(r)`$ of Fig. 2(a) yield $`c=3.75`$ and the circular flow velocity $`v_s=2n/R12`$ (in units of $`b_{}\omega _{}`$). Because $`v_s>c`$, the circular flow is supersonic and energetically unstable.

Impurities and walls can cause dissipation for a supersonic superflow. For example, Ref. Jackson  used the Gross-Pitaevskii equation for bulk condensates and ordinary trapped condensates to show that the dissipation of the superflow is caused by vortex-pair creation at an impurity object. This was shown experimentally Raman . However, this does not apply to the present case of quasi-one-dimensional superflow because of the tight radial confinement omake . To show the influence of impurities, we introduce an impurity potential of Gaussian form as $`V_{\mathrm{imp}}=V_0\mathrm{exp}[\{(xx_0)^2+(yy_0)^2\}/\sigma ^2]`$ (at rest in the laboratory frame), where $`\{x_0,y_0\}`$ is the impurity position. Figure 4 shows an example of the destabilization process of the giant vortex state by $`V_{\mathrm{imp}}`$, for which $`V_0`$ is chosen to be larger than $`V_{\mathrm{trap}}(R)`$ and $`\sigma `$ smaller than the width of the ring condensate. Because $`v_s`$ exceeds the Landau critical velocity $`c`$, the impurity excites the longitudinal density waves along the ring condensate, thus producing a drag force. Due to the presence of dissipation and the rotating drive, the system reaches a steady state as shown in Fig. 4(b). Therefore, the quadratic-plus-quartic system offers a scheme for studying the stability of a quasi-one-dimensional superflow.

In conclusion, a fast rotating BEC in a quadratic-plus-quartic potential can generate a giant vortex and circular superflow around it by absorbing singly quantized vortices into a central low density hole. The circular superflow is quasi-one-dimensional and goes at a supersonic speed. We thus hope that this system will encourage experiments to study low dimensional superfluidity in atomic-gas BECs.

M.T. acknowledges support by a Grant-in-Aid for Scientific Research (Grant No. 12640357) by the Japan Society for the Promotion of Science. M.U. acknowledges support by a Grant-in-Aid for Scientific Research (Grant No. 11216204) by the Ministry of Education, Science, Sports, and Culture of Japan, and by the Toray Science Foundation.


text:
# Subgap anomaly and above-energy-gap structure in chains of diffusive SNS junctions

\[

## Abstract

We present the results of low-temperature transport measurements on chains of superconductor–normal-constriction–superconductor (SNS) junctions fabricated on the basis of superconducting PtSi film. A comparative study of the properties of the chains, consisting of 3 and 20 SNS junctions in series, and single SNS junctions reveals essential distinctions in the behavior of the current-voltage characteristics of the systems: (i) the gradual decrease of the effective suppression voltage for the excess conductivity observed at zero bias as the quantity of the SNS junctions increases, (ii) a rich fine structure on the dependences $`dV/dIV`$ at dc bias voltages higher than the superconducting gap and corresponding to some multiples of $`2\mathrm{\Delta }/e`$. A model to explain this above-energy-gap structure based on energy relaxation of electron via Cooper-pair-breaking in superconducting island connecting normal metal electrods is proposed.

\]

In the past few years, the mesoscopic systems, consisting of a normal metal (N) or heavily doped semiconductor being in contact with a superconductor (S), have attracted an increased interest mainly because of the richness of the involved quantum effects . The key mechanism governing the carrier transport through the NS contact is the Andreev reflection . In this process, an electron-like excitation with an energy $`ϵ`$ smaller than the superconducting gap $`\mathrm{\Delta }`$ moving from the normal metal to the NS interface is retro-reflected as a hole-like excitation, while a Cooper pair is transmitted into the superconductor. This phenomenon is the basis of the proximity effect, which generally implies the influence of a superconductor on the properties of a normal metal when being in electrical contact. The consequences of Andreev reflection on the current-voltage characteristics (CVC) of a NS junction are studied in detail in the so-called BTK theory , which describes the subgap current transport in terms of ballistic propagation of quasiclassical electrons through the normal-metal region, accompanied by Andreev and normal reflections from NS interface. The probability of Andreev reflection and normal reflection are energy-dependent quantities and the relation between them depends on the barrier strength at NS interface, that is characterized by a parameter $`Z`$ ranging from $`0`$ for a perfect metallic contact to $`\mathrm{}`$ for a low transparency tunnel barrier. For a perfect contact ($`Z=0`$), probability of Andreev reflection is equal to unity for particles with an energy $`ϵ`$ smaller than the superconducting gap $`\mathrm{\Delta }`$ and the subgap conductance are found to be twice the normal state conductance thus demonstrating the double charge transfer. In the other limit, when $`Z\mathrm{}`$, the conductivity is very small.

When a normal metal is placed between two superconducting electrods another mechanism is involved in charge transfer. It is multiple Andreev reflection process (MAR). The concept of MAR was first introduced by Klapwijk, Blonder and Tinkham  in order to explain the subharmonic energy gap structure observed as dips in the differential resistance $`dV/dI`$ of SNS junctions at the voltages $`V_n=2\mathrm{\Delta }/ne`$ ($`n=1,2,3,\mathrm{}`$). In all systems, diffusive as well as ballistic, the MAR mechanism relies on the quasiparticles to be able to add up energy gained from multiple passages. This energy gain results in strong quasiparticle nonequilibrium. Calculations of the CVC, within the approach developed in Ref. , and taking into account MAR process was performed in  (OTBK theory). It was shown that MAR process results not only in SGS at voltages $`V_n=2\mathrm{\Delta }/ne`$ but affects the general form of CVC at all voltages as well. Although BTK theory is not intended to describe diffusion transport in normal-metal region, a suitable value $`Z0.55`$ of the barrier strenght gives result very close to those obtained by Green’s-function method for microconstriction in the dirty limit . Experimentally, in some cases BTK theory has been successful in providing a method to extract the junction reflection coefficients based on experimental CVC .

At present, a large variety of SNS junctions (most of them are diffusive) is being fabricated and studied . The investigations of these junctions are primarily focused on the nonlinear behavior of the current-voltage characteristics, which exhibit the anomalous resistance dip at zero bias \[zero bias anomaly (ZBA)\], the subharmonic energy gap structure (SGS), etc. These experiments have revealed a number of peculiar properties of diffusive SNS junctions, the explanation of those being beyond the ballistic BTK and OTBK theories. Nowadays the properties of diffusive SNS junctions are considered to be determined not only by the parameters of the NS interface, but nonlocal coherent effects as well, namely (i) the superposition of multiple coherent scattering at the NS interfaces in the presence of disorder (so-called reflectionless tunneling ) and (ii) electron-electron interaction in the normal part. The latter is one of the important points of a recently developed “circuit theory” when applied to diffusive superconducting hybrid systems. Within this approach, which is based on the use of nonequilibrium Green’s functions, the electron-electron interaction induces a weak pair potential in the normal metal. It results not only in a change of the resistance, but also in a non-trivial distribution of the electrostatic and chemical potentials in the structure, which implies non-local resistivity. Following the spirit of Nazarov’s circuit theory Bezuglyi et al.  have developed a consistent theory of carrier transport in long diffusive SNS junctions with arbitrary transparency of NS interfaces. Although much work both theoretical and experimental has been done on single SNS junctions, it is a challenge to fabricate and carry out comparative measurements on multiply connected SNS systems. In this paper, we present the results of low-temperature transport measurements on chains of SNS junctions and on single SNS junctions and perform the comparative analysis of their properties.

The design of our samples is based on the fabrication technique of SNS junctions that has been recently proposed and realized by us for preparation of single SNS junctions and two-dimensional arrays of SNS junctions . The main idea employed in those experiments was to use superconductive and normal region made of the same material. The point is the suppression of superconductivity in the submicron constrictions made in an ultrathin polycrystalline PtSi superconducting film by means of electron beam lithography and subsequent plasma etching. Here we use the same fabrication technique to prepare chains of SNS junctions.

The original PtSi film (thickness – 6 nm) was formed on a Si substrate. The film was characterized using Hall bridges $`50\mu `$m wide and $`100\mu `$m long. The film had a critical temperature $`T_\mathrm{c}=0.56`$ K. The resistance per square was 104 $`\mathrm{\Omega }`$. The carrier density obtained from Hall measurements was $`710^{22}`$ cm-3, corresponding to a mean-free path $`l=1.2`$ nm and a diffusion constant $`D=6`$ cm2/s, estimated using the simple free-electron model.

Single SNS junction is a constriction between two holes made in the film and placed in the Hall bridge. The hole diameters are $`1.7\mu `$m and the distance between centres of holes is $`2.1\mu `$m, resulting in a width of narrowest part of the constriction of $`0.4\mu `$m. A scanning electron micrograph and a schematic view of a sample are presented in Fig.1a,b. The chains of SNS junctions are a series of constrictions connected by the islands of the film (Fig.1c,d). They are designed to provide the possibility of the comparative study. The dimensions of constrictions are identical with those of the single SNS junction, and the characteristic size of the islands of the film is $`1.3\mu `$m. As the constrictions are not superconducting we have a chain of SNS junctions (Fig.1e).

The measurements were carried out with the use of a phase sensitive detection technique at a frequency of 10 Hz that allowed us to measure the differential resistance ($`dV/dI`$) as a function of the dc current ($`I`$). The ac current was equal to 1–10 nA. Figure 2 shows typical dependences of $`dV/dI`$-$`V`$ for the samples with single constriction (1D-1) and chains, consisting of 3 (sample 1D-3) and 20 (sample 1D-20) SNS junctions in series. On the abscissa the value of voltage obtained by numerical integration of the experimental dependences $`dV/dI`$-$`I`$ and afterward divided by the number of SNS junctions is plotted. Measured differential resistance is also divided by the number of SNS junctions (by 3 for structure 1D-3 and by 20 for structure 1D-20). Thus for the chains data presented in Fig.2 are dependences the differential resistance as a function of dc bias voltage per SNS junction in the average. As is seen at high voltage for chains the average resistances of SNS unit of chains obtained in such a way are close to each other and to the value $`dV/dI`$ of the structure with single SNS junction. A common feature for all structures studied is a minimum of the differential resistance at zero bias voltage. It points out the high transparency of NS interfaces in these samples. The data for single SNS junction (1D-1) exhibit a behavior very similar to that reported in our previous work . In Ref.  the differential resistance of long diffusive SNS junction with perfect NS interfaces at zero voltage was estimated by the following expression:

$$dV/dI(0)=R_N(12.64\xi _N/L),$$
(1)
where $`\xi _N=\sqrt{\mathrm{}D/2\pi kT}`$ is the decay lenght for the pair amplitude, $`L`$ is the lenght of the normal-metal region, and $`R_N`$ is its resistance at $`T>T_\mathrm{c}`$. The expression (1) has obtained on condition that the inelastic lenght $`l_ϵ`$ exceeds the junction lenght $`L`$. The energy relaxation is described by the time $`\tau _ϵ^1=\pi ϵ^2/(8\mathrm{}E_F)`$. This is the time it takes for a “hot” quasiparticle with energy $`ϵ`$ much larger than temperature $`T`$ to thermalize with all the other electrons. For the excitation energy $`ϵ=2\mathrm{\Delta }`$ ($`\mathrm{\Delta }/e270`$ $`\mu `$V) we obtain $`l_ϵ=\sqrt{D\tau _ϵ}4.6`$ $`\mu `$m. The magnitude of $`l_ϵ`$ really exceeds the characteristic lenght of the constriction. As $`R_N`$ we take the difference between the resistance of the whole structure with constriction and the resistance of the original film without constriction measured between the same probes at $`T>T_\mathrm{c}`$. It gives $`R_N530`$ $`\mathrm{\Omega }`$. Assuming $`L1`$ $`\mu `$m, we find from (1) at $`T=100`$ mK the differential resistance at zero bias $`dV/dI(0)0.8R_N=420`$ $`\mathrm{\Omega }`$. As is seen in Fig.2, this value is close fit to measured $`dV/dI(0)`$ for the sample with single constriction 1D-1. A clear-cut distinction between dependences $`dV/dI`$-$`V`$ of the single SNS junction and the chains consists in non-monotone behavior of $`dV/dI`$-$`V`$ characteristics of the latter. The differential resistance of the chains has a minimum at zero bias voltage and shows a maximum at a finite bias voltage of about 250 $`\mu `$V for the sample 1D-3 and 30 $`\mu `$V for 1D-20. The second feature is the gradual decrease of the effective suppression voltage for the excess conductivity observed at zero bias as the quantity of the SNS junctions increases. The same behavior – non-monotonic $`dV/dI`$-$`V`$ characteristics and the decrease of the effective suppression voltage for the excess conductivity in comparison of that in single SNS junctions – has been observed in two-dimensional arrays of SNS junctions .

The most unexpected result obtained for chains of SNS junction is the presence of a fine and fully symmetrical structure as dips in the differential resistance at high biases (Fig.3a,b). The positions of these dips approximately correspond to 4, 6, and 12 multiples of the superconducting gap. As is seen in Fig.3c, the temperature dependence of the positions of these dips actually reflects the dependence $`\mathrm{\Delta }(T)`$. To our knowledge the results presented in this paper is the first investigation of chains of SNS junctions and the first observation such above energy gap structure (AGS). To offer an explanation of the AGS let us consider a current driven in the normal region at high voltages. As the inelastic lenght more than the normal metal region, at finite voltage quasiparticle passing through N1 between two superconductors S1 and S2 (Fig.4b,d) results in a nonequilibrium energy distribution of quasiparticles. For a start we address this problem to the ballistic OTBK theory . In this approach the quasiparticles are divided into two subpopulation which depend on their direction of motion $`f_{}(E)`$ and $`f_{}(E)`$, with all energies measured with respect to the local chemical potential. The current through the junction is

$$I_{\mathrm{}}^+\mathrm{}[f_{}(E)f_{}(E)]𝑑E.$$
(2)
In Fig.4a,c we have plotted $`f_{}(E)`$, $`f_{}(E)`$, and $`[f_{}(E)f_{}(E)]`$, assuming a parameter $`Z=0.55`$. The voltage across the junction equals $`4\mathrm{\Delta }/e`$ and $`6\mathrm{\Delta }/e`$ in Fig.4a and Fig.4c, respectively. The function $`[f_{}(E)f_{}(E)]`$ is depicted by the blackness (high values are black) in Fig.4b,d, with the energy scale being plotted vertically. These calculations clearly show that the nonequilibrium distributions are sharply peacked at the four gap edges. What happens to quasiparticle injected to the superconducting electrode S2 and having the energy more than $`2\mathrm{\Delta }`$ measured with respect to the gap edge of S2? They can decay under spontaneous phonon emission into states of lower energy . Phonons emitted in the process of energy relaxation of injected quasiparticles can be reabsorbed via Cooper-pair-breaking. Depending on the volume of finite states, the probability of this process is non-monotone function of the primary energy of injected quasiparticles. It takes the peak values at the energies multipled to $`2\mathrm{\Delta }`$, for density of quasiparticle states to have the singularity at the gap edge. As a result one quasiparticle injected to S2 with energy high the gap edge by $`m2\mathrm{\Delta }`$ can cause up to $`1+2m`$ quasiparticles injected to N2 (Fig.4b,d) and consequently the increase of the differential resistance. The realization of this scenario requires the lenght of energy relaxation to be comparable to dimensions of superconductor. It is really fulfilled in the SNS systems studied. Proposed mechanism supposes the sistem consisting of two SNS junctions in series to be sufficient in order to observe the AGS.

In summary, we performed the first comparative study of the properties of the chains and single SNS junctions. Our experiments reveal essential distinctions in the behavior of the current-voltage characteristics of these systems. A detailed quantitative analysis needs to take into account the contributions of nonlocal coherent transport and an effect of quasiparticle injection in both the superconducting and normal regions of chains of the SNS junctions.

We thank R. Donaton and M. R. Baklanov (IMEC) for providing us the PtSi films, and A. E. Plotnikov for the performing of electron lithography. We acknowledge useful discussions with V. F. Gantmakher, V. V. Ryazanov, Ya. G. Ponomarev, and M. Feigel’man. This work has been supported by the program “Low-dimensional and mesoscopic condenced systems” of the Russian Ministry of Science, Industry and Technology and by RFBR Grant No. 00-02-17965.


text:
# 1 Introduction

## 1 Introduction

The quark-gluon plasma (QGP) has been expected to appear at high temperature and/or high density, due to asymptotic freedom of QCD. Indeed, not a few of the experimental results upto CERN SPS incident energies suggest the existence of the QGP phase , and the recent experiment at BNL RHIC may give us cleaner and richer signals .

Theoretically, in the lattice QCD calculations , no phase transition is observed for the two light or two light and one medium light flavor cases at zero baryon density. In both cases only a continuous cross over is observed for realistic quark mass values.

On the other hand, in more phenomenological approach, state equations of QGP and hadrons are assumed, and the transition point is determined by Gibbs condition. In this approach, the cross over cannot be reproduced, and that may give rise to a reason to doubt its validity. However, the authors believe that the gross structure of the phase diagram obtained in this approach should still remain valid far away from the cross over region. It should also be noted that at a finite baryon chemical potential, not much is known from lattice QCD because of a well known technical difficulty which makes the Monte Carlo technique inapplicable.

The purpose of the present paper is to investigate what type of the state equation of hadrons is preferable in the phenomenological approach, where the discussion is devoted exclusively to gross phase structure of hadrons and QGP, neglecting the multi-quark states. It means that the detail structure around phase transition or cross over region is neglected.

In the phenomenological approach, there are following problems concerning the gross structure of phase diagrams. When a simple ideal gas of nucleons is used for the equation of hadron state and MIT bag model is used for the equation of QGP state, hadron phase appears at high density and at zero temperature (see, e.g. ). This is because

$$p_h^0/p_q^027\text{as}\mu _B\mathrm{},(T=0)$$
(1.1)
where $`p_h^0`$ and $`p_q^0`$ denote pressures of ideal nucleon phase and (ideal) quark phase, respectively, and $`\mu _B`$ is baryochemical potential. To circumvent this difficulty, volume exclusion effects of hadrons of van der Waals type were taken into account (see, e.g. , see also ). However, the treatments in  are thermodynamically inconsistent , for example, $`n(p/\mu )_T`$, $`s(p/T)_\mu `$. Although the formulation in  is thermodynamically consistent, it has still a difficulty that pressure diverges as $`n1/v_0`$, where $`n`$ is number density and $`v_0`$ is fixed hadron volume, and acausality emerges (there is no relativistic rigid body) as the authors themselves pointed out. Thus we were led to ‘soft core’ model, that is the compressible bag model .

In addition to the difficulty that hadron phase appears at high density and at zero temperature, an equation of simple ideal gas of hadrons again suffers from the same difficulty at another region, that is, at high temperature and zero baryon number density . The argument in  is simple. At $`\mu _B=0`$ and in high temperature region, the pressure $`p_q`$ of QGP phase and the pressure $`p_h`$ of hadron phase are given by

$$p_q=g_q(1/90\pi ^2)T^4B,p_h=g_h(1/90\pi ^2)T^4,$$
(1.2)
where $`g_q`$ and $`g_h`$ are the degeneracy factors in each phase and $`B`$ is a bag constant. When we consider a mixed gas of lower lying nonstrange hadrons, that is, $`\pi `$, $`\rho `$, $`N`$ and $`\mathrm{\Delta }`$, the degeneracy factor $`g_h`$ becomes $`3+9+(7/4)(4+16)=47`$ and it exceeds $`g_q=28+(7/4)43=37`$. Thus hadron phase appears again. In the compressible bag model, however, the effective degeneracy factor of hadron phase is greatly reduced so that QGP phase is realized at high temperature even if infinitely many nonstrange mesons are taken into account .

The authors of  suggested that the compressible bag model gives the phase diagram free from above difficulties, but it has not shown explicitly: In  the equations of state only at $`T=0`$ are used for discussing neutron stars, and in  used are those at $`\mu _B=0`$ for multiple production in heavy ion collisions. In  the phase diagram with finite $`T`$ and $`\mu _B`$ is calculated, but rich hadron spectrum is neglected.

The unique purpose of the present paper is, then, to show that the compressible bag model indeed gives the naively expected phase diagram in all regions. As far as the authors know, it is the only model that satisfies the followings:

It is built in thermodynamically consistent formalism.

It is valid in very high temperature regions even if rich hadron spectrum is taken into account.

It is explicitly calculable in whole regions, and it gives well-behaved phase diagrams, i.e. there are no extra hadron phase.

Since the hadron level structure affects the phase diagram, we consider the following two models:

Hadron phase consists of $`N\overline{N}\pi `$ system while QGP phase consists of nonstrange quarks, their anti-quarks and gluons.

Hadron phase consists of a system of $`102`$ species of hadrons of which masses are less than $`10.6\mathrm{GeV}`$, except for hadrons with top quark . QGP phase consists of u-, d-, s-, c-, b-quarks and their anti-quarks and gluons.

Model I is a model with a few hadrons, while model II is an example which includes many hadrons.

While in  the continuous level-density function of  is used, in this paper the phase diagrams are calculated with real spectrum of hadrons in order to avoid the model dependence. The effect of the cutoff in the hadron spectrum in model II is discussed in the final section.

In section 2, for comparison the phase structure in free point-like model is shown with the above two types of models, where one can explicitly see the difficulties explained in this section. In section 3, the equation of state in the compressible bag model  is briefly reviewed, and with the resulting phase diagram it is shown that the difficulty is indeed removed. Final section is devoted to concluding remarks.

## 2 A difficulty in the phase structure for models of free point-like hadrons

In this section, we assume that all hadrons, quarks, and gluons are free point-like objects and discuss their phase structure. For this purpose, we present expressions of pressures of free point-like particles in order to establish notations in the first subsection. On the basis of these expressions, the phase diagram is discussed and a difficulty is pointed out in the second subsection.

### 2.1 Pressures for hadron gas and QGP

Total pressure $`p`$ of a mixed gas of $`N`$-species of free point-like particles is given by

$$p(T,\mu _B)=\underset{n=1}{\overset{N}{}}p_n(T,\mu _n),$$
(2.1)
where $`T`$ ($`\mu _B`$) is the temperature (baryochemical potential) of the gas and $`p_n`$ ($`\mu _n`$) is the pressure (baryochemical potential) of the $`n`$-th species of particles. The pressure $`p_n`$ is given by

$`p_n`$ $`=`$ $`\eta _ng_nT{\displaystyle \frac{d^3𝒌}{(2\pi )^3}\mathrm{log}\left[1+\eta _n\mathrm{exp}\left(\frac{E_n\mu _n}{T}\right)\right]},`$ (2.2)
$`E_n`$ $`=`$ $`\sqrt{𝒌^2+m_n^2},`$ (2.3)
$`\mu _n`$ $`=`$ $`a_n\mu _B,`$ (2.4)
where $`\eta _n`$ is a statistical factor that takes the value $`1`$ for bosons and $`+1`$ for fermions, respectively. The quantities $`g_n`$, $`a_n`$, and $`m_n`$ denote the degeneracy factor, the baryon number, and the mass of the $`n`$-th species of particles, respectively.

In low temperature region, $`p_n`$ is expanded as follows;

$`p_n`$ $`=`$ $`{\displaystyle \frac{g_n}{24\pi ^2}}\left[\mu _nk_{nf}^3{\displaystyle \frac{3}{2}}m_n^2\mu _nk_{nf}+{\displaystyle \frac{3}{2}}m_n^4\mathrm{log}\left({\displaystyle \frac{\mu _n+k_{nf}}{m_n}}\right)\right]`$ (2.6)
$`+{\displaystyle \frac{g_n}{12}}\mu _nk_{nf}T^2+{\displaystyle \frac{7\pi ^2g_n}{720}}\left[{\displaystyle \frac{3\mu _n}{k_{nf}}}\left({\displaystyle \frac{\mu _n}{k_{nf}}}\right)^3\right]T^4,`$
$`k_{nf}`$ $`=`$ $`\sqrt{\mu _n^2m_n^2},`$ (2.7)
when the $`n`$-th species of particles are fermions. In the above $`k_{nf}`$ is a fermi momentum. Eq. (2.6) is applicable for $`\mu _n>m_n`$. Otherwise the pressure $`p_n`$ vanishes identically.

Now let us present expressions for pressures of hadrons and QGP. As for total pressure $`p_h`$ of a mixed gas of $`N`$-species of free point-like hadrons, it is given by

$$p_h(T,\mu _B)=\underset{i=1}{\overset{N}{}}p_{h_i}(T,\mu _i),$$
(2.8)
where $`p_{h_i}`$ follows (2.2)$``$(2.7). In the followings we consider two cases for hadronic systems; that is, the model I and the model II which are defined in the introduction.

On the other hand, total pressure of QGP $`p_q`$ is given by

$$p_q(T,\mu _B)=\underset{j=1}{\overset{N^{}}{}}p_{q_j}(T,\mu _j)B,$$
(2.9)
in the bag model. Here $`p_{q_j}`$ follows (2.2)$``$(2.7) and $`B`$ is a bag constant. The index $`j`$ denotes the species of particles, and we consider two cases, the model I and the model II.

### 2.2 Phase structure

Given the expressions of pressures of hadron phase and QGP phase, we can consider their phase structure. Critical points are determined by the following Gibbs condition;

$$p_h(T,\mu _B)=p_q(T,\mu _B).$$
(2.10)
In Fig. 1, $`p_h`$ and $`p_q`$ are plotted as functions of $`\mu _B^4`$, where $`T`$ is fixed to zero. In Fig. 2, they are plotted as functions of $`T^4`$, where $`\mu _B`$ is fixed to zero. In both figures, calculations are performed for the model I. As shown in Fig. 1, there is no critical point for $`B^{1/4}>0.1536\mathrm{GeV}`$. For $`B^{1/4}=0.1536\mathrm{GeV}`$, there is one critical point. For $`0.1534\mathrm{GeV}B^{1/4}<0.1536\mathrm{GeV}`$, there are two critical points. For $`B^{1/4}<0.1534\mathrm{GeV}`$, there is one critical point since $`p_h`$ vanishes for $`\mu _B<m_N`$ and unphysical is an apparent critical point in the region $`\mu _B<m_N`$, where $`m_N`$ is nucleon mass. Therefore, at $`T=0`$, abnormal hadron phase always appears at high densities.

The reason for this difficulty is as follows. In higher density region of $`\mu _Bm_N`$, masses are negligible so that pressures $`p_h`$ and $`p_q`$ are approximated as

$`p_h`$ $``$ $`g_h\mu _B^4,`$ (2.11)
$`p_q`$ $``$ $`g_q\mu _q^4,`$ (2.12)
where $`g_h`$ ($`g_q`$) is a statistical degree of freedom of hadrons (QGP). Since $`\mu _q=(1/3)\mu _B`$, (2.12) reads

$$p_q\frac{g_q}{81}\mu _B^4.$$
(2.13)
Then $`p_h`$ becomes large faster than $`p_q`$ as $`\mu _B`$ increases, since an effective statistical degree of freedom of QGP is greatly reduced. In other words, hadrons get a share of energy three times larger than QGP, when $`\mu _B`$ increases, so that $`p_h`$ becomes larger than $`p_q`$. As a consequence hadron phase appears in high density region.

Since the difficulty mentioned above stems from the qualitative nature of the model of hadrons that their statistical degree of freedom is too large, we should consider to modify the model of hadrons. So far we consider the model I at present, the situation will become worse in the model II, since the statistical degree of freedom becomes further large.

In Fig. 3, $`p_h`$’s for the model I and the model II are plotted as functions of $`\mu _B`$, where $`T`$ is fixed to zero and in Fig. 4, $`p_h`$’s for those are plotted as functions of $`T`$, where $`\mu _B`$ is fixed to zero. As seen in Fig. 3 (Fig. 4), $`p_h`$ for the model II becomes considerably larger than $`p_h`$ for the model I in high density (temperature) region, as expected.

In Fig. 5, $`p_h`$ and $`p_q`$ are plotted as functions of $`\mu _B^4`$ with $`T=0`$. In Fig. 6, they are plotted as functions of $`T^4`$ with $`\mu _B=0`$. In both figures, calculations are performed for the model II. Qualitative features in Fig. 1 and Fig. 5 are the same. However, qualitative features in Fig. 2 and Fig. 6 are different. As shown in Fig. 6, hadron phase appears in high temperature region. The resulting phase diagram for the model I  is shown in Fig. 7 and that for the model II is shown in Fig. 8.

The situation for the model II has become worse than for the model I. In order to modify this situation, we should replace the model of hadrons with alternative ones. As one of them, we formulate the compressible bag model in the next section.

## 3 The compressible bag model

First, we briefly summarize the results of . Let us suppose a mixed gas of $`n`$-species of hadrons enclosed in volume $`V`$ at finite temperature $`T`$. In the compressible bag model, the free energy function $`\widehat{F}`$ of the gas is given by

$`\widehat{F}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}F_f(N_i,V^{},T,M_i(v_i)),`$ (3.1)
$`V^{}`$ $`=`$ $`Vb{\displaystyle \underset{j=1}{\overset{n}{}}}N_jv_j,`$ (3.2)
where $`N_i`$ is a number of $`i`$-th species of hadrons and $`v_i`$ is its volume and $`M_i`$ its mass. The constant $`b`$ is a volume exclusion efficiency parameter. The function $`F_f`$ is a free energy function of free point-like hadron gas. As for mass function $`M_i(v_i)`$, we assume that of MIT bag model:

$`M_i(v_i)=A_iv_i^{1/3}+Bv_i,`$ (3.3)
where $`B`$ is bag constant.

Under the approximation that the average of inverse Lorentz factor of $`i`$-th species of hadrons equals to unity ($`\gamma _i^11`$) in the rest frame of the system, basic requirement of the compressible bag model that $`\widehat{F}/v_i=0`$ and the requirement that the chemical potential of hadron should be, $`\mu _i=\widehat{F}/N_i=a_i\mu _B`$, where $`a_i`$ is the baryon number of the $`i`$-th species of hadrons, determine the pressure $`p`$ of the system as a function of $`T`$ and $`\mu _B`$ by the following equation:

$`p`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}\eta _ig_iT{\displaystyle \frac{d^3𝒌}{(2\pi )^3}\mathrm{log}\left\{1+\eta _i\mathrm{exp}[(E_i\mu _i^{}(p,m_i))/T]\right\}},`$ (3.5)
$`E_i=\sqrt{𝒌^2+M_i(p,m_i)^2},`$
$`M_i(p,m_i)`$ $`=`$ $`m_i\left(1+{\displaystyle \frac{3bp}{4B}}\right)\left(1+{\displaystyle \frac{bp}{B}}\right)^{3/4},`$ (3.6)
$`\mu _i^{}(p,m_i)`$ $`=`$ $`\mu _ibv_ip=a_i\mu _B{\displaystyle \frac{bm_ip}{4B}}\left(1+{\displaystyle \frac{bp}{B}}\right)^{3/4},`$ (3.7)
where $`g_i`$ is a degeneracy factor of the $`i`$-th species of hadrons and $`m_i=4(A_i/3)^{3/4}B^{1/4}`$ is its mass in the vacuum.

Here two comments are in order. First $`p`$ determined as above depends on $`b`$ and $`B`$ only in the combination of $`b/B`$ since $`M_i`$ and $`\mu _i^{}`$ depend on $`b`$ and $`B`$ in that combination as seen in (3.6) and (3.7). Second the pressure in the compressible bag model does not become so large in high density and/or high temperature region in contrast to the one in the model of free point-like hadrons. This is because masses of hadrons become large in high density and/or high temperature region in the compressible bag model so that kinetic energies of hadrons increase rather slowly and the pressure of hadrons do not become so large. Thus the compressible bag model have a chance to evade the difficulty discussed in the previous section. In the following, we assure that the difficulty is indeed removed, by numerical calculations.

In order to do numerical calculations, we have to fix bag constant $`B`$ and volume exclusion parameter $`b`$. The parameter $`b`$ is determined by the relation

$`bv_N(p=0)={\displaystyle \frac{bm_N}{4B}}={\displaystyle \frac{4\pi }{3}}(0.82\mathrm{fm})^3,`$ (3.8)
when the bag constant is given. The value $`0.82\mathrm{fm}`$ is the proton charge radius. It should be noted that equations of state for hadron phase are determined uniquely at this stage since the equation of states for hadron phase only depend on $`b/B`$. We tentatively put

$$B^{1/4}=0.25\mathrm{GeV}$$
(3.9)
and then get $`b=5.0`$. The resulting $`p_h`$ and $`p_q`$ are shown as functions of $`\mu _B`$ ($`T`$) in Fig. 9 (Fig. 10). We can recognize that the compressible bag model does not suffer from the difficulty that exist for free point-like hadrons. This point is more clearly seen by the phase diagrams shown in Fig. 11 with $`B^{1/4}=0.25\mathrm{GeV}`$. Critical temperature $`T_c`$ and critical baryochemical potential $`\mu _c`$ are given by

$`T_c=0.18\mathrm{GeV},\mu _c=2.2\mathrm{GeV}\text{in model I},`$ (3.10)
$`T_c=0.17\mathrm{GeV},\mu _c=1.8\mathrm{GeV}\text{in model II}.`$ (3.11)
It is noted that the state equation of hadrons are determined by the ratio $`b/B`$ and that of QGP by $`B`$, the critical values, then, depends on $`B`$ and $`b`$. If the constraint (3.8) is taken, the values depend only on $`B`$. As a reference, some cases are shown below;

$`T_c=0.14\mathrm{GeV},\mu _c=1.7\mathrm{GeV}\text{in model I},`$ (3.12)
$`T_c=0.14\mathrm{GeV},\mu _c=1.4\mathrm{GeV}\text{in model II},`$ (3.13)
for $`B^{1/4}=0.20\mathrm{GeV}`$ and

$`T_c=0.22\mathrm{GeV},\mu _c=2.7\mathrm{GeV}\text{in model I},`$ (3.14)
$`T_c=0.21\mathrm{GeV},\mu _c=2.1\mathrm{GeV}\text{in model II},`$ (3.15)
for $`B^{1/4}=0.30\mathrm{GeV}`$.

From Figs. 9$``$11 one can see that the compressible bag model gives an expected and reasonable phase diagram in whole regions. This result is stable in the sense that one can choose the bag constant from rather wide range.

## 4 Concluding remarks

In this paper, it has been shown by explicit calculations that the compressible bag model gives well-behaved phase diagram in whole regions, even if many hadron states are taken into account. In the present calculation it is true that, even in the model II, the infinite series of hadrons are cut off at finite mass, but the compressible bag model may give an expected phase diagram if infinite series of hadrons are included. This is because, in the model, the masses of hadrons increase and the effect of higher mass state is much suppressed in high-temperature or high-density regions. Indeed in  a continuous level-density function are used for hadrons and it is shown that the abnormal hadron phase does not appear in the compressible bag model, although limited at $`\mu _B=0`$.

Free point-like models produce unnatural hadron state in high-temperature or high-density regions as shown in Sections 1 and 2. A way to avoid this difficulty is the compressible bag model, and it is worth examining in various phenomenological analyses.