text:
# Abundances in Stars from the Red Giant Branch Tip to the Near the Main Sequence Turn Off in M71: II. Iron Abundance Based on observations obtained at the W.M. Keck Observatory, which is operated jointly by the California Institute of Technology and the University of California

## 1 INTRODUCTION

Abundance determinations of stars in Galactic globular clusters can provide valuable information about important astrophysical processes such as stellar evolution, stellar structure, Galactic chemical evolution and the formation of the Milky Way. Surface stellar abundances of C, N, O, and often Na, Mg, and Al are found to be variable among red giants within a globular cluster. The physical process responsible of these star-to-star element variations is still uncertain (see Cohen et al., 2001, Paper I).

Of particular importance to the present study are the results of King et al. (1998), who found that \[Fe/H\] among M92 subgiants is a factor of two smaller than \[Fe/H\] from red giants in the same cluster (Cohen, 1979; Sneden et al., 1991). If this result is not due to systematic differences arising from the analysis procedures of the different groups that handle the red giant branch (RGB) and the subgiant samples, then the results of King et al. (1998) would suggest some modification of photospheric Fe abundances which would be quite difficult to explain. Clearly the present day clusters stars are incapable of modifying their own Fe abundances through nucleosynthesis. Yet evolution-driven changes in \[Fe/H\] could be possible if the outer envelopes of the subgiants were somehow infused with Fe-poor external material. Thus as evolution progesses up the RGB the deeping convective envelope would dilute the photospheric \[Fe/H\] with more Fe-rich “unmodified” material from the interior causing an increase in \[Fe/H\] with luminosity. We note the unlikelyhood of this scenario as Fe, unlike the lighter elements, is not involved in normal mass transfer processes that might ocurr in binary systems or in normal stellar winds, and is only produced in supernovae. Moreover, it is difficult to understand how cluster stars could be contaminated by material with less Fe, presumably at a time after the present subgiants had formed.

In order to study the origin of the star-to-star abundance variations and to address the issues raised by King et al. (1998), we have started a program to determine chemical abundances of the nearer galactic globular cluster stars. In this paper, we present our results for the iron abundance of M71, the nearest globular cluster reachable from the northern hemisphere. Similar programs are underway at ESO by Castilho et al. (2000) and Gratton et al. (2001) taking advantage of the fact that the nearest globular cluster accessible from a southern site (NGC 6397) is 1.8 times close than M71 and has a lower reddening. Our M71 sample includes stars over a large range in luminosity: 19 giant stars, 3 horizontal branch stars, and 3 stars near the main sequence turnoff, in order to study in a consistent manner red giants, horizontal branch stars, and stars at the main sequence turnoff. Details on the star sample, observations, data reduction and determination of stellar parameters are described in Paper I. Previous high dispersion abundance analysis for M71 involve studies of red giants only, resulting in \[Fe/H\] of $``$0.70 (Cohen, 1983), $``$0.80 (Gratton et al., 1986), between $``$0.6 and $``$1.0 (Leep et al., 1987), and $``$0.79 (Sneden et al., 1994).

## 2 ANALYSIS

We begin our analysis with Fe, as many Fe lines are identified in our HIRES spectra over a wide range of excitation potentials and line strengths, as well as over two ionization states (Fe I and Fe II). The iron abundance analysis is done using a current version of the LTE spectral synthesis program MOOG (Sneden, 1973). A line list specifying the wavelengths, excitation potentials, gf–values, damping constants, and equivalent widths for the observed Fe I and Fe II lines is required. The provenance of the gf–values and the measurement of equivalent widths are discussed below. The damping constants for all Fe I and Fe II lines were set to twice that of the Unsöld approximation for van der Waals broadening following Holweger et al. (1991). The use of the Blackwell approximation for the damping constants gives the same result within the errors, when comparing \[Fe/H\] obtained with the “good line set” of Fe I (see definition of line sets on Sec. 2.2).

In addition, a model atmosphere for the effective temperature and surface gravity appropriate for each star and a value for the microturbulent velocity are also required. We use the grid of model atmospheres from Kurucz (1993b) with a metallicity of \[Fe/H\] = $``$0.5 dex, based on earlier high dispersion abundance analysis of M71 red giants (Cohen, 1983; Gratton et al., 1986; Leep et al., 1987; Sneden et al., 1994). The final result for \[Fe/H\] is not sensitive to small changes in the metallicity of the model atmosphere. In particular, the error introduced by using a model with \[Fe/H\] = $``$0.5 instead of $``$0.7 dex is very small (see Tables 2 & 3). The effective temperatures and surface gravities are derived from the photometry of the stars as described in Paper I. The error in the photometric $`T_{eff}`$ is 75 K for giants and 150 K for the dwarfs and the error in the photometric log($`g`$) is 0.2 dex (Paper I). The microturbulent velocity is derived spectroscopically (see below). The stellar parameters are listed in Table 1.

### 2.1 Transition Probabilities

Transition probabilities for the Fe I lines were obtained from several laboratory experiments, including studies of Fe I absorption lines produced by iron vapor in a carbon tube furnace (Blackwell et al., 1979, 1982a, 1982b, 1986) (Oxford Group), measurement of radiative lifetimes of Fe I transitions by laser induced fluorescence (O’Brian et al., 1991; Bard et al., 1991; Bard & Kock, 1994), Fe I emission line spectroscopy from a low current arc (May et al., 1974), and emission lines of Fe I from a shock tube (Wolnik et al., 1971). We also considered solar gf–values from Thévenin (1989, 1990) when needed.

We compare the gf–values obtained by the different experiments in an attempt to place them onto a common scale with respect to the results from O’Brian et al. (1991), who provided the longest list of gf–values. We considered for the comparison only the set of lines present in our data, which have a wavelength coverage roughly from 5380 Å to 7900 Å. We found that the values of O’Brian et al. (1991) and of the Oxford Group were on the same scale; the mean difference in log($`gf`$) between the two experiments is $`0.02\pm 0.01`$ for 21 lines in common. Similar results are found when considering 21 lines in common between O’Brian et al. (1991) and Bard et al. (1991) (see also Bard & Kock, 1994), where the mean difference in log($`gf`$) is $`0.02\pm 0.04`$. Considering 34 lines in common between O’Brian et al. (1991) and May et al. (1974), the mean difference in log($`gf`$) is $`0.03\pm 0.03`$. When comparing the 14 lines in common between O’Brian et al. (1991) and Wolnik et al. (1971) we found a mean difference in log($`gf`$) of $`0.07\pm 0.02`$. We also compared the results from O’Brian et al. (1991) with solar gf–values obtained by Thévenin (1989, 1990). We found that the mean difference in log($`gf`$) is $`+0.05\pm 0.02`$ when comparing 68 lines in common between these latest works. Similar offsets are found by Lambert et al. (1996).

The gf–values for our Fe I lines were taken when possible from laboratory data in the following order of priority: from O’Brian et al. (1991), from the Oxford Group, from May et al. (1974), and from Wolnik et al. (1971) corrected by $`0.07`$ dex. In the cases where no laboratory data were available, we used solar gf–values from Thévenin (1989, 1990) corrected by +0.05 dex.

Transition probabilities for the Fe II lines were taken from the solar analysis of Blackwell et al. (1980), Biémont et al. (1991), and from the semiempirical calculations of Kurucz (1993b). While restricting ourselves to lines present in our spectra, we compared the gf–values from Blackwell et al. (1980) and Kurucz (1993b) to the study of Biémont et al. (1991). There are 9 lines in common between Biémont et al. (1991) and Blackwell et al. (1980). The mean difference between the two in log($`gf`$) is $`0.14\pm 0.02`$. A similar offset was found by Lambert et al. (1996). We found no significant difference between the results of Biémont et al. (1991) and Kurucz (1993b), since the 9 lines in common result in a mean difference of $`0.03\pm 0.02`$.

The gf–values for our Fe II lines were taken in the following order of priority: from Biémont et al. (1991), from Blackwell et al. (1980) corrected by 0.14 dex, and from Kurucz (1993b).

### 2.2 Measurement of Equivalent Widths

Our sample contains many stars observed in mutiple orders, with many detectable absorption features in these high S/N spectra. For example, in the coolest M71 star in our sample, M71-1-45, 1407 absorption lines have been identified. A FORTRAN code to automatically search for absorption features and measure their equivalent width ($`W_\lambda `$), EWDET, was developed for this project. The code is available upon request to SVR. EWDET determines the continuum location of the HIRES spectra by fitting a curve to the spectra performing several iterations of point rejection above and below sigma levels given by the user. Then, EWDET identifies lines above the noise level defined by a factor of two of the continuum rms dispersion. Each of the identified lines is fit by a Gaussian profile and then the $`W_\lambda `$ are computed by the integration of the fitted Gaussian. The error in $`W_\lambda `$ is computed by adding quadratically the error at each point of integration, $`\sigma _i`$, times the step of the integration. The error at each point of integration is given by:

$$\sigma _i^2=g(\lambda _i)^2\times \left[\frac{\sigma _P^2}{P^2}+\frac{(\lambda _i\lambda _{cen})^2}{\sigma ^4}\sigma _{\lambda _{cen}}^2+\frac{(\lambda _i\lambda _{cen})^4}{\sigma ^6}\sigma _\sigma ^2+\sigma _{cont}^2\right]$$

where $`g(\lambda _i)`$ is the Gaussian profile, $`P`$ is the peak of the Gaussian, $`\sigma _P`$ is the error in the peak of the Gaussian, $`\lambda _{cen}`$ is the central wavelength of the Gaussian, $`\sigma _{\lambda _{cen}}`$ is the error in the central wavelength, $`\sigma `$ is the dispersion of the Gaussian, $`\sigma _\sigma `$ is the error in the dispersion of the Gaussian, and $`\sigma _{cont}`$ is the error in the continuum. The errors of the Gaussian parameters are from the covariance matrix of the Gaussian fit. The expression for the error in $`W_\lambda `$ is deduced by propagating the errors of the Gaussian parameters and assuming that the continuum level is equal to one (see Appendix). The fit by a Gaussian profile is reasonable even for the strongest lines we use, as shown in Figure 1, where the observed line is plotted with a solid curve and its corresponding Gaussian profile fit is shown with a dashed curve.

The line list identified and measured by EWDET is then correlated to the line list with the atomic parameters to specifically identify the Fe I and Fe II lines. The detailed lists of $`W_\lambda `$ and gf–values will be given in the next paper.

The spectral resolution, $`\lambda /\mathrm{\Delta }\lambda `$, of an echelle is fixed, unlike a low incidence, angle of low order grating spectrograph where $`\mathrm{\Delta }\lambda `$ is constant, independent of $`\lambda `$. Hence a line of constant central depth, $`D`$, will have an equivalent width proportional to $`\lambda D`$. We construct the relationship between $`\lambda D`$ and $`W_\lambda `$ for the Fe I lines to look for possible blends and for saturation effects. The $`\lambda D`$ versus $`W_\lambda `$ relation for three stars are plotted in Figure 2; M71-I is one of the most luminous and coolest stars of our sample, M71-G53476\_4543 is a medium luminosity and temperature star, and M71-G53392\_4624 is one of the faintest and hottest stars of the sample. We fit a second order polynomial to the $`\lambda D`$$`W_\lambda `$ relationship for each star, performing several iterations of point rejection above and below the two sigma level. The second order term is needed to account for line saturation at large $`W_\lambda `$ and the rejection of points is needed to eliminate blended lines. The second order fit is plotted as a solid line in Figure 2, and the points considered for such a fit are shown in black.

For Fe I, we use two sets of lines. The first set, subsequently called “the weak line set”, contains those Fe I lines which are within two sigma levels of the $`\lambda D`$$`W_\lambda `$ fit, have $`W_\lambda `$ $`<`$ 60 mÅ, and have errors less than a third of the $`W_\lambda `$. This set of lines produces a sample of the best weak Fe I lines with the most accurate $`W_\lambda `$ and the abundances derived using them will have a minimal dependence on the choice of microturbulent velocity. The second set, subsequently called “the good line set”, consists of all the Fe I lines with errors less than a third of the $`W_\lambda `$ and with $`W_\lambda `$ computed from the fit for $`\lambda D`$$`W_\lambda `$ determined for each star. In future papers, the $`W_\lambda `$ of lines of other elements will be determined in the same manner as the Fe I lines of “the good line set”. This way a consistent comparison can be done among the resulting abundances without a restriction on the strength of the lines used. The “weak line set” is different for each star. Actually the weak lines for the stars near the main sequence are no longer weak for the cooler stars in our sample. We compare the results from 20 lines common to 15 stars over almost the whole range in $`T_{eff}`$ with the results from the “weak line set”; there is no difference within the errors, nor a trend in $`T_{eff}`$.

For Fe II, the $`W_\lambda `$ of the lines are also determined using the fit to the $`\lambda D`$$`W_\lambda `$ relation of the Fe I lines of “the good line set”. The Fe II lines follow the relationship determined from Fe I lines well, as shown in Figure 3. Additional Fe II lines, not picked up automatically, were measured by hand for the stars near the turn off and for M71-G53425\_4612 and M71-G53457\_4709. The set of $`W_\lambda `$ for the hand selected Fe II lines is computed from the $`\lambda D`$$`W_\lambda `$ relation of the Fe I lines, after determining their observed depth from the spectra.

The number of Fe I lines, for both sets of lines, and the number of Fe II lines utilized in this analysis are listed in Table 1.

### 2.3 Spectroscopic Effective Temperature

The effective temperature ($`T_{eff}`$) of a star can be determined spectroscopically by requiring the abundance to be independent of the lower excitation potential. This technique can be applied to 20 of our stars where we have Fe I lines with enough range in lower excitation potential to do so. For the spectroscopic $`T_{eff}`$ determination we are using “the weak line set” of Fe I lines, because its resulting abundance and spectroscopic $`T_{eff}`$ will be only weakly dependent on the choice of microturbulent velocity. We find that the spectroscopic $`T_{eff}`$ is in good agreement with the photometric $`T_{eff}`$ derived in Paper I, as shown in Figure 4. The solid line in Figure 4 shows the ideal case when the spectroscopic and the photometric $`T_{eff}`$ are equal. The scatter around the solid line is about 150 K, which is comparable to the error of the photometric $`T_{eff}`$ of 75 K for giants and of 150 K for dwarfs (Paper I), also shown in Figure 4. The scatter around the solid line is symmetric, not above or below, indicating the lack of systematic effects with the photometric temperatures.

### 2.4 Microturbulent Velocity

The microturbulent velocity ($`\xi `$) of a star can be determined spectroscopically by requiring the abundance to be independent of the strength of the lines measured as the equivalent width. We apply this technique for both sets of Fe I lines. The resulting $`\xi `$ and the \[Fe/H\] computed with it for “the weak line set” of Fe I lines are listed in Table 1. Only 18 of our stars have enough weak Fe I lines to derive $`\xi `$ spectroscopically. We also compute $`\xi `$ using all the good Fe I lines for 20 of our stars. The difference of the resulting \[Fe/H\] computed with the derived $`\xi `$ for the two set of lines is plotted with respect to $`T_{eff}`$ in Figure 5. The mean difference is $`0.004\pm `$ 0.011, hence, as expected, the \[Fe/H\] results from both sets of lines show a very good agreement.

The relationship between $`\xi `$ determined for the set of all good Fe I lines and the photometric $`T_{eff}`$ is shown in Figure 6. The solid line corresponds to a linear least squares fit of the data, excluding the red horizontal branch (RHB) stars, marked with circles. The best fit line is given by:

$$\xi =3.304.66\times 10^4\times T_{eff}$$

The scatter around the solid line is about 0.2 km s-1, which is a reasonable estimation of the error in $`\xi `$. For the rest of the analysis, we will use the set of all the good Fe I and Fe II lines, with $`\xi `$ computed from the $`\xi `$-$`T_{eff}`$ fit. For the RHB stars we use a value of 1.61 km s-1, which corresponds to the mean value determined for the three RHB stars. The microturbulent velocity used for our stars is listed in Table 1.

## 3 RESULTS

Given the stellar parameters from Table 1, we determined the iron abundance using the equivalent widths of all the good Fe I and Fe II lines identified in the HIRES spectra. We employ the grid of stellar atmospheres from Kurucz (1993b) to compute the iron abundance using the four stellar atmosphere models with the closest $`T_{eff}`$ and log($`g`$) to each star’s parameters. The \[Fe/H\] listed in Table 1 is an interpolation of the results from the closest stellar model atmospheres to the appropriate $`T_{eff}`$ and log($`g`$) for each star.

### 3.1 \[Fe/H\] from Fe I lines

The results of \[Fe/H\] from Fe I lines are listed in column 9 of Table 1 and plotted against the photometric $`T_{eff}`$ in Figure 7a. $`T_{eff}`$ is used for the x-axis as a convenient parameter for characterizing the position of the stars in the color-magnitude diagram as it also ranks the stars in luminosity (except for the RHB stars). The errors listed in Table 1 correspond to the larger of the statistical uncertainty, given by the standard deviation of the iron abundance from different lines divided by the square root of the number of lines, or a minimum value of 0.03 dex. These errors are lower limits to the actual uncertainties in the abundances, since they do not include uncertainties due to the stellar parameters nor any systematic effects that might be present. We estimate the sensitivity of \[Fe/H\] derived from Fe I lines with respect to the stellar parameters in three cases 4250/1.0/1.0, 5000/2.5/1.0 and 5500/4.0/1.0, where the three numbers correspond to $`T_{eff}`$/log($`g`$)/$`\xi `$. The results are listed in Table 2, where the range adopted for each parameter is representative of its uncertainty. Our determination of \[Fe/H\] from Fe I lines is most sensitive to errors in $`T_{eff}`$, which is less than $``$ 0.1 dex for $`\mathrm{\Delta }`$$`T_{eff}`$ of $`\pm `$ 100 K, and have a minimal sensitivity on the choice of metallicity of the model atmosphere grid for plausible changes in \[Fe/H\] ($`\pm `$0.2 dex).

The solid line, shown in Figure 7a, is a linear fit weighted by the errors of \[Fe/H\] versus $`T_{eff}`$. The slope of the fit is $`(0.8\pm 3.6)\times 10^5`$ dex/K, which is consistent with \[Fe/H\] being constant, independent of $`T_{eff}`$ (ie, of luminosity or equivalently position in the color-magnitude diagram). We divide our sample in four groups of stars: giant stars at or above the RHB, stars on the RHB, giant stars below the RHB, and main sequence stars near the turnoff. The mean \[Fe/H\] for each group is listed on Table 4. We found no significant difference in the mean \[Fe/H\] obtained from Fe I lines among the defined groups of stars.

The mean \[Fe/H\] weighted by the errors of all 25 stars is $`0.71\pm 0.08`$, in very good agreement with earlier determinations (Cohen, 1983; Gratton et al., 1986; Leep et al., 1987; Sneden et al., 1994).

### 3.2 \[Fe/H\] from Fe II lines

The determinations of \[Fe/H\] from Fe II lines are listed in column 11 of Table 1 and plotted against the photometric $`T_{eff}`$ in Figure 7b. The errors listed in Table 1 corresponds to the statistical uncertainty or a value of 0.05 dex, whichever is larger. We estimate the sensitivity of \[Fe/H\] derived from Fe II ines with respect to the stellar parameters in the same manner as the sensitivity of \[Fe/H\] from Fe I lines. The results are listed in Table 3, where the range adopted for each parameter is representative of its uncertainty. We see a stronger sensitivity on the stellar parameters from the Fe II lines than from the Fe I lines. The \[Fe/H\] determination from Fe II lines is most sensitive to the systematic error (note that the internal uncertainty in log($`g`$) is $``$ 0.1 dex) in log($`g`$), as well as to $`T_{eff}`$ among the coolest M71 giants. The sensitivity on the choice of metallicity of the model atmosphere grid is small for reasonable changes in metallicity.

The solid line, shown in Figure 7b, is a linear fit weighted by the errors of \[Fe/H\] versus $`T_{eff}`$. The slope of the fit is $`(+3.1\pm 5.2)\times 10^5`$ dex/K, which is consistent with \[Fe/H\] being constant, independent of $`T_{eff}`$. We found that there is no significant difference in the mean \[Fe/H\] obtained from Fe II among stars from different luminosity groups, listed in Table 4. The mean \[Fe/H\] weighted by the errors is $`0.84\pm 0.12`$, in very good agreement with our result from Fe I lines and earlier determinations (Cohen, 1983; Gratton et al., 1986; Leep et al., 1987; Sneden et al., 1994).

### 3.3 NLTE effects

The iron abundance could be affected by departures from LTE. The main NLTE effect in late-type stars is caused by overionization of electron donor metals by ultraviolet radiation (Auman & Woodrow, 1975). Recently, Gratton et al. (1999) and Thévenin & Idiart (1999) studied NLTE effects in Fe abundances in metal-poor late-type stars. Gratton et al. (1999) found that NLTE corrections for Fe lines are very small in dwarfs of any $`T_{eff}`$, and only small corrections ($`<`$ 0.1 dex) are expected for stars on the red giant branch. Thévenin & Idiart (1999) found that NLTE corrections become less important as \[Fe/H\] increases, being less than 0.1 dex for stars with \[Fe/H\]$`>`$0.75 dex, and that ionized lines are not significantly affected by NLTE.

One way to explore possible NLTE effects present in our data is by comparing the results from Fe I and Fe II lines. The difference between \[Fe/H\] from Fe II and Fe I lines is plotted in Figure 8 against $`T_{eff}`$. The solid line is a linear fit weighted by the errors. The slope of the fit is $`(+2.0\pm 8.2)\times 10^5`$ dex/K, which is nearly flat. The mean difference is $`0.13\pm 0.18`$. We conclude that NLTE effects are negligible in our iron abundance determination, as expected from results of earlier studies (Gratton et al., 1999; Thévenin & Idiart, 1999).

## 4 DISCUSSION

Our \[Fe/H\] abundance results involve for the first time a wide luminosity sample of stars, which includes at the same time stars from luminous giants to stars near the turnoff. We find that the \[Fe/H\] abundance, from both Fe I and Fe II lines, is independent of $`T_{eff}`$, and equivalently luminosity.

Our result is in agreement with the work of Gratton et al. (2001). They present abundances from high dispersion spectra from the VLT of stars in NGC 6397 and NGC 6752. They found that the \[Fe/H\] obtained for stars at the base of the subgiant branch agrees within a few percent with the \[Fe/H\] obtained for stars at the main sequence turnoff, and further compare this value with analysis of the RGB stars in this cluster by other groups. Note that the luminosity range of the sample presented in our work is several orders of magnitude wider than the luminosity range of Gratton et al. (2001)’s sample.

Our results, and those of Gratton et al. (2001), appear to be in disagreement with inhomogeneities in \[Fe/H\] found earlier by King et al. (1998). They obtained \[Fe/H\]=$``$2.52 dex for a sample of subgiant stars in M92, which is a factor of two smaller than \[Fe/H\] measurements using red giants in the same cluster (Cohen, 1979; Sneden et al., 1991). King et al. (1998) compare their result for the M92 subgiants with analysis of RBG stars by other groups, who may have determined the stellar parameters and performed the abundance determinations in a different way. This possible difference in the analysis of the giant and subgiant sample may account for the difference in \[Fe/H\] found by King et al. (1998) or perhaps the determination of the stellar parameters by King et al. (1998) is flawed. Our result, on the other hand, is robust, because we have determined both the stellar parameters and the Fe abundance in a homogeneous and consistent manner for all our stars.

Gratton et al. (2001) also found that NGC 6397 is homogeneous in both O and Fe, while an O-Na anticorrelation is present among unevolved stars in NGC 6752, which is very difficult to explain by the deep mixing scenario. Lines from many additional species, including O, Na, Mg, Ti, Sc among others, are observed in our HIRES spectra. We plan to present additional information in the matter of light, iron-peak, and heavy elements in the near future.

## 5 CONCLUSIONS

We present results of a high dispersion analysis of Fe I and Fe II lines to obtain \[Fe/H\] for 25 members of the Galactic globular cluster M71. Our sample of stars includes 19 giant stars, 3 horizontal branch stars, and 3 stars near the turnoff. Our conclusions are summarized as follows:

* The \[Fe/H\] obtained from Fe I lines agrees very well with the \[Fe/H\] obtained from Fe II lines.
* The mean \[Fe/H\] obtained from Fe I and Fe II lines of all 25 stars is in good agreement with earlier determinations.
* The \[Fe/H\] obtained from both Fe I and Fe II lines is independent of $`T_{eff}`$, and equivalently luminosity.
* No difference is found among the mean \[Fe/H\] from giant stars located at or above the RHB, RHB stars, giant stars located below the RHB and stars near the main sequence turnoff.

In the near future, we will present the result from the analysis now underway of additional elements.

The entire Keck/HIRES and LRIS user communities owes a huge debt to Jerry Nelson, Gerry Smith, Steve Vogt, Bev Oke, and many other people who have worked to make the Keck Telescope and HIRES and LRIS a reality and to operate and maintain the Keck Observatory. We are grateful to the W. M. Keck Foundation for the vision to fund the construction of the W. M. Keck Observatory. We thank R. Gratton for providing a detailed description of his automatic equivalent width measuring program and R. Pogge for providing the Gaussian profile fitting routine. Partial support to MMB was provided by a Theodore Dunham, Jr. grant for Research in Astronomy and the National Science Foundation under grant AST-9624680 to MMB and grant AST-9819614 to JGC.

## Appendix A Errors in the Equivalent Widths

The equivalent width of a line is defined as:

$$W_\lambda =\frac{g(\lambda )}{cont(\lambda )}𝑑\lambda =\underset{i}{}\frac{g(\lambda _i)}{cont(\lambda _i)}\mathrm{\Delta }\lambda $$

where, $`g(\lambda )`$ is the flux of the spectrum, $`cont(\lambda )`$ is the continuum level, and $`\mathrm{\Delta }\lambda `$ is the step that can be made arbitrarily small to better approximate the integral. In this case, $`g(\lambda )`$ is the Gaussian profile of the line, given by:

$$g(\lambda )=P\times exp\left[\frac{(\lambda \lambda _{cen})^2}{2\sigma ^2}\right]$$

where $`P`$ is the peak of the Gaussian, $`\lambda _{cen}`$ is the central wavelength, and $`\sigma `$ is the dispersion of the Gaussian. If the equivalent width is computed through out the summation then the error in the equivalent width, $`\sigma _w`$ will be the quadratic summation of the errors in the individual points times the step:

$$\sigma _{W_\lambda }^2=\underset{i}{}(\sigma _i\times \mathrm{\Delta }\lambda )^2$$

where:

$$\sigma _i^2=\left[\frac{w}{g}\sigma _{gi}\right]^2+\left[\frac{w}{cont}\sigma _{cont}\right]^2$$

$$\sigma _i^2=\left[\frac{1}{cont(\lambda _i)}\sigma _{gi}\right]^2+\left[\frac{g(\lambda _i)}{cont(\lambda _i)^2}\sigma _{cont}\right]^2$$

If we assume that the continuum level is equal to one then:

$$\sigma _i^2=\sigma _{gi}^2+g(\lambda _i)^2\sigma _{cont}^2$$

Since the Gaussian depends on three parameters:

$$\sigma _{gi}^2=\left[\frac{g}{P}\sigma _P\right]^2+\left[\frac{g}{\lambda _{cen}}\sigma _{\lambda _{cen}}\right]^2+\left[\frac{g}{\sigma }\sigma _\sigma \right]^2$$

$$\sigma _{gi}^2=\left[\frac{g(\lambda _i)}{P}\sigma _P\right]^2+\left[\frac{(\lambda _i\lambda _{cen})}{\sigma ^2}g(\lambda _i)\sigma _{\lambda _{cen}}\right]^2+\left[\frac{(\lambda _i\lambda _{cen})^2}{\sigma ^3}g(\lambda _i)\sigma _\sigma \right]^2$$

$$\sigma _{gi}^2=g(\lambda _i)^2\times \left[\frac{\sigma _P^2}{P^2}+\frac{(\lambda _i\lambda _{cen})^2}{\sigma ^4}\sigma _{\lambda _{cen}}^2+\frac{(\lambda _i\lambda _{cen})^4}{\sigma ^6}\sigma _\sigma ^2\right]$$

Finally:

$$\sigma _i^2=g(\lambda _i)^2\times \left[\frac{\sigma _P^2}{P^2}+\frac{(\lambda _i\lambda _{cen})^2}{\sigma ^4}\sigma _{\lambda _{cen}}^2+\frac{(\lambda _i\lambda _{cen})^4}{\sigma ^6}\sigma _\sigma ^2+\sigma _{cont}^2\right]$$


format: text

text:
# Superconducting Gap and Strong In-Plane Anisotropy in Untwinned YBa2Cu3O7-δ

\[

## Abstract

With significantly improved sample quality and instrumental resolution, we clearly identify in the ($`\pi `$,0) ARPES spectra from YBa2Cu3O6.993, in the superconducting state, the long-sought ‘peak-dip-hump’ structure. This advance allows us to investigate the large $`a`$-$`b`$ anisotropy of the in-plane electronic structure including, in particular, a 50% difference in the magnitude of the superconducting gap that scales with the energy position of the hump feature. This anisotropy, likely induced by the presence of the CuO chains, raises serious questions about attempts to quantitatively explain the YBa2Cu3O7-δ data from various experiments using models based on a perfectly square lattice.

\]

High-temperature superconductivity (HTSC) is intimately related to the CuO2 plane, which is the only common structural feature in all cuprates. This fact has led most of the proposed microscopic theories to assume a CuO2 square planar structure. However, for the practical reason of sample quality, some of the most important and defining experiments have been performed on YBa2Cu3O7-δ (Y123), which does not have a square lattice, but rather an orthorhombic structure ($`b/a1.015`$), caused by the presence of a CuO chain layer . This orthorhombicity, according to LDA calculation , should result in significant anisotropy in the in-plane electronic structure (this term will be used throughout this paper to refer to the electronic states associated with the CuO2 plane), making it problematic to compare theories based on a square lattice with experimental data from Y123. Therefore, it is crucial to quantify the effect of orthorhombicity, if any, on the in-plane electronic structure in Y123. The problem is that angle-resolved photoemission spectroscopy (ARPES), being a uniquely powerful tool for this important task, has, until now, not been particularly effective for the study of Y123 . The important ‘peak-dip-hump’ structure, which is seen routinely in Bi2Sr2CaCu2O8+δ (Bi2212)  has never been observed in Y123. This absence, together with the presence of a surface state , raises questions about ARPES data from Y123, and the universality of the superconducting peak in the cuprates.

This paper reports a breakthrough in this important issue, made possible by significantly improved sample quality and instrumental resolution. By isolating a surface state peak near the Fermi energy ($`E_F`$), we can clearly resolve a ‘peak-dip-hump’ structure in the ARPES spectra around ($`\pi ,0`$) in Y123 that resembles the superconducting peak observed in Bi2212 . More significantly, we find a strong $`a`$-$`b`$ asymmetry of the in-plane electronic structure, such as the superconducting gap magnitude, which differs by about 50%. We argue that such a strong in-plane $`a`$-$`b`$ anisotropy should be taken into account when interpreting experiments performed on Y123.

ARPES experiments were carried out at beamline 5-4 at SSRL, which is equipped with a normal-incidence-monochromator and a SCIENTA SES-200 analyzer. Untwinned YBa2Cu3O6.993 single crystals (overdoped, $`T_c`$ = 89 K) with superior chemical purity ($`99.9999.995\%`$) and crystallinity (FWHM = 0.007 in the X-ray rocking curve) were grown by the self-flux method in BaZrO3 crucibles . Single crystals were cleaved in-situ at 10 K with a base pressure better than 5x10-11 torr. ARPES spectra were recorded at a photon energy of 28 eV with an energy resolution of 10 meV and an angular resolution of 0.3, corresponding to a momentum resolution of 1.5% of the Brillouin zone (BZ). Low energy electron diffraction (LEED) patterns were routinely acquired at the end of experiments to check the quality of the cleaved surfaces and to confirm the sample orientation.

Fig. 1 presents the energy distribution curves (EDCs) along high symmetry directions in the BZ. Data were taken at 10 K under two distinct sample orientations with respect to the polarization of the incoming radiation (insets), which are referred to as O4 and O2, following the existing notation . In order to show the dispersion clearly, in particular for the broad and weak features, the intensity plots of the second derivatives of the EDCs are displayed in Fig. 2a-2f. Four features can be clearly identified, two of which appear at X as sharp peaks close to $`E_F`$, along with two broad features at higher binding energy (BE). As we will justify in the following paragraphs, we assign them, from the lower to higher BE, as surface state, superconducting peak, hump, and chain state. Fig. 2g summarizes all of the band dispersions.

The most pronounced feature in Fig. 1 is a narrow, intense peak just below $`E_F`$ at both X and Y, whose high sensitivity to surface degradation at elevated temperatures indicates its surface state character. The precise origin of this surface state is, however, not well understood at present. One possibility could be the surface termination effect. An earlier scanning tunneling microscopy (STM) study  demonstrated that Y123 cleaves between the CuO chain and BaO layers, resulting in both CuO chain-terminated and BaO-terminated areas on the cleaved surface. The scenario of a chain-related surface state is favored by a previous ARPES study , and seems to be consistent with more recent STM data . However, this assignment is difficult to reconcile with the two-dimensional dispersion of this feature along $`\mathrm{\Gamma }`$-X and $`\mathrm{\Gamma }`$-Y, as seen in the current data. For this paper, as we will focus on the bulk electronic structure, we leave it as an open question.

For the broad feature located well below $`E_F`$ at X, instead of a bonding $`\sigma `$ state as suggested in Ref. , we reassign it to a chain derived state. This assignment is concluded from the quasi one-dimensional (1D) character of this feature, i.e., the dispersion is strong along the chain direction ($`\mathrm{\Gamma }`$-Y and X-S) and weak in the perpendicular direction ($`\mathrm{\Gamma }`$-X). Further support comes from the fact that the two crossings along $`\mathrm{\Gamma }`$-Y and X-S agree very well with the 1D chain Fermi surface (FS) determined by a full FS mapping. The dispersion of this chain state and the related chain FS are consistent with the LDA band structure . It is also interesting to note that this feature looks rather similar to the chain state observed in PrBa2Cu3O7  and that the dispersion along the chain direction qualitatively agrees with the holon band predicted by the $`t`$-$`J`$ model calculation at quarter filling . However, we were not able to identify in the current data a spinon band, which may possibly be masked by the hump around 0.12 eV in Y123.

Now let us focus on the superconducting peak and the hump that are clearly evident in Fig. 1a. For better illustration, the EDC at X is replotted in Fig. 3a with fitting curves. In the fitting procedure, an integration background that smoothly merges into the high-energy hump is subtracted to extract the double peak feature, which is then fitted with a simple spectral function: a product of a Fermi function and two Lorentzians, convolved with a gaussian instrumental resolution function. After subtracting the surface state peak, the rest of the spectrum looks strikingly similar to the well-known ‘peak-dip-hump’ structure observed in the ARPES spectra of Bi2212 at ($`\pi ,0`$) below $`T_c`$ (Fig. 3b), indicating that the second peak may be a superconducting peak.

More direct evidence of the superconducting nature of this peak comes from its temperature dependence (Fig. 4a). Due to the unstable nature of the Y123 surface at high temperatures, great caution has been taken in cycling the temperature. A LEED pattern taken after the measurements (Fig. 4b) shows no visible difference from that of a freshly cleaved surface, indicating a well-ordered surface, although some minor aging effects can be detected. Following the same fitting procedure (with both peak positions fixed), we extract the superconducting peak and plot in Fig. 4c the normalized superconducting peak ratio (SPR), which is defined after Ref.  as the ratio between the extracted peak intensity and the total spectral weight (excluding the surface state peak) integrated over \[0.3 eV, -0.1 eV\]. Despite the scatter in the data and the large error bars, the general trend suggests that the appearance of this peak is related to the superconducting transition.

Furthermore, the remarkable resemblance of the spectra in Fig. 3b implies that this ‘peak-dip-hump’ structure at X is a CuO2 plane derived feature since the CuO chain is absent in Bi2212. This argument is further supported by data taken on Zn-doped ($``$2%) Y123 with the same oxygen content under identical experimental conditions (Fig. 3c). Compared with the data from Zn-free Y123 (Fig. 3a), the major change is that the superconducting peak shifts to lower BE, i.e., the gap size decreases upon Zn-doping. Since Zn atoms only occupy the Cu sites in the CuO2 plane , this change in the superconducting peak suggests that it is directly related to the CuO2 plane. Due to the weak dispersion of this peak around ($`\pi ,0`$), we can approximately quantify the superconducting gap using the peak position, giving a gap size of 25 and 29 meV at X for Zn-doped ($`T_c`$ = 79 K) and Zn-free ($`T_c`$ = 89 K) Y123, respectively. This indicates that, as expected, the gap magnitude scales with $`T_c`$.

At this point, one may speculate as to why the ‘peak-dip-hump’ feature is absent at Y in Fig. 1b. This is partly due to the presence of the intense surface state peak, which masks the underlying bulk electronic states close to $`E_F`$. In fact, the data taken at Y with the polarization perpendicular to the chain direction (Fig. 3d), where the surface state peak is suppressed, clearly show a second peak at higher BE. In analogy with the assignment of the superconducting peak at X we may reasonably interpret this additional feature as evidence for a superconducting peak at Y. The hump associated with this peak is almost indiscernible in Fig. 3d, but becomes more evident moving away from Y to S as it disperses towards $`E_F`$ (Fig. 2f). The temperature dependent measurements on this feature are very difficult due to its weak intensity and the presence of the strong surface state peak. Hence, we deliberately aged the sample by warming it to the point where the surface state peak completely vanishes (Fig. 4d). A small peak on top of a relatively large background exhibits a reproducible temperature dependence expected for a superconducting peak. As a consequence of aging, this peak shifts to lower BE ($``$33 meV), as compared with Fig. 3d. Since the Y-point is outside the 1D chain FS, the peak and the hump detected at Y should be exclusively related to the CuO2 plane.

Having identified all of the features in the low energy excitation spectra, our data provide the first observation of a ‘peak-dip-hump’ structure near ($`\pi ,0`$) in a system other than Bi2212. This finding partially resolves the longstanding issue regarding the universality of the superconducting peak, which is crucial for a complete understanding of HTSC as it appears to contain information on both pairing and phase coherence . The failure to detect this peak in previous photoemission studies on Y123 was mainly due to the presence of the intense surface state peak and insufficient energy and momentum resolution. However, it is not clear at present whether the absence of this peak in other superconducting cuprates, such as Bi2Sr2CuO6+δ and La2-xSrxCuO4, is simply a consequence of lower $`T_c`$ , i.e., low superfluid density, or is also related to the lack of the CuO2 bilayer that exists in both Bi2212 and Y123. Further experiments on other cuprate superconductors, particularly the single-plane Tl- or Hg-based compounds with higher $`T_c`$, are necessary to clarify this issue.

Now let us turn our attention to the strong $`a`$-$`b`$ anisotropy of the electronic structure as revealed by a direct comparison between Fig. 1a and 1b. For the chain state, the $`a`$-$`b`$ asymmetry is in accord with the 1D character of the CuO chain. More interesting is the $`a`$-$`b`$ anisotropy associated with the peak and the hump as they are related to the superconductivity in the CuO2 plane. The major asymmetry lies in the energy positions of the peak and the hump: 29 ($`\mathrm{\Delta }_x`$) and 120 meV ($`\omega _x`$) for X, and 44 ($`\mathrm{\Delta }_y`$) and 180 meV ($`\omega _y`$) for Y. This yields a significant ($``$50%) difference in the gap magnitude between X and Y, a substantial deviation from the ideal $`d`$-wave pairing state. This $`a`$-$`b`$ asymmetry of the superconducting gap in Y123 has also been reported by Limonov and co-workers in their Raman scattering study , from which a value of 25 and 30 meV can be derived for $`\mathrm{\Delta }_x`$ and $`\mathrm{\Delta }_y`$, respectively. It is also remarkable to see that the gap magnitude scales with the hump positions: $`\mathrm{\Delta }_x/\omega _x\mathrm{\Delta }_y/\omega _y0.24`$. The scaling of these two energy scales has also been observed by ARPES in Bi2212 as a function of doping . The observation of this scaling behavior in the same material along two different directions imposes strong constraints on theory.

In addition to the asymmetry of the gap magnitude, we also notice that there is a clear difference in the intensities of the peak and the hump between X and Y. According to Ref. , this could possibly be a signature of the in-plane anisotropy of the superfluid density. However, the polarization effects due to the presence of the CuO chains on the surface layer may complicate the situation, making a quantitative comparison difficult. We note that strong $`a`$-$`b`$ anisotropy of the penetration depth (i.e., the superfluid density) was also found in microwave and far infrared spectroscopic measurements , which was interpreted as evidence for the presence of superconductivity in the chain layer. While it is still a controversial issue whether the ground state of the CuO chain is superconducting or a charge density wave , our data offer an alternative interpretation: this anisotropy may in fact reside in the CuO2 plane itself.

Such a strong $`a`$-$`b`$ anisotropy is crucial for a comprehensive understanding of the superconducting properties of Y123 and the interpretation of related experiments. It is therefore important to understand the origin of this in-plane anisotropy. One possible explanation could be based on the orthorhombicity of the CuO2 plane in Y123. According to the LDA band calculations, the ratio between the hopping integrals along the $`a`$ and $`b`$ directions scales roughly as $`(b/a)^4`$, which, in turn, yields considerable differences in the calculated electronic structures between $`\mathrm{\Gamma }`$-X-S and $`\mathrm{\Gamma }`$-Y-S  that qualitatively agree with our data. However, no direct conclusion can be drawn for the superconducting state properties. Alternatively, a strong chain-plane coupling may be present. In particular, as recently suggested by Atkinson , a nontrivial $`d`$-wave gap structure, which seems to be consistent with our data, exists on the anti-bonding band related FS due to the hybridization with the chain band.

We thank for J. C. Davis, A. L. de Lozanne, H. Eisaki, T. Mizokawa and S. Schuppler for fruitful discussions, T. Tohyama, and S. Maekawa for unpublished results from their calculations. The data presented here were obtained at SSRL, which is operated by the Department of Energy, Office of Basic Energy Sciences. The Stanford work was also supported by ONR Grant N00014-98-1-0195 and NSF grant DMR-9705210.


text:
# Absorption Line Shape of a One-Dimensional Bose Gas

## Abstract

We discuss the line shape for an internal state transition for bosonic atoms confined as a one-dimensional gas. Typical line shape is an edge singularity due to the absence of Bose-Einstein Condensation in such systems.

PACS numbers: 03.75.Fi

Recent progress in trapped bosonic gases has stimulated a lot of activities in the study of many-body effects in such systems. Most theoretical papers concentrate on effects due to the existence of the Bose-Einstein Condensate (BEC) and the associated macroscopic wave-function. With current techniques, it is likely that one can as well study an effectively one-dimensional (1D) trapped gas . This situation occurs when the bosons are put in a trap which is tightly confining in two directions ( $`y`$ and $`z`$) but is essentially free in the third ($`x`$), and with particles occupying only the lowest quantized subband for motions in the $`y`$-$`z`$ plane. This system is unique in that, for any finite interaction among the bosons, there is no Bose-Einstein Condensation even at zero temperature. . The long-wavelength quantum fluctuations destroy the phase coherence of the system in the macroscopic limit. For a given density $`n`$, the occupation number $`N_0`$ of the lowest single particle state, though in general increases with the size and hence the total number of particles $`N`$ in the system, is nevertheless only a negligible fraction of the total ($`\genfrac{}{}{0pt}{}{\mathrm{lim}}{N\mathrm{}}\frac{N_0}{N}0`$ ). Though many properties of 1D bose gas have been studied in the past , the subject remains somewhat academic since there exists no physical realization. The possibility of realizing this 1D bose gas has already stimulated some more theoretical work on this system recently

The present 1D bosonic system is closely related to its fermionic counterpart; and they are collectively known as the “Luttinger liquids” . One-dimensional fermionic systems have been studied much more extensively in the context of quantum wires, carbon nanotubes  etc. There the one-dimensional nature leads to the absence of fermi liquid behavior, so the properties of the system is qualitatively different from its non-interacting counterpart.

The atomically trapped 1D bose gas offers us an excellent new opportunity to study the peculiar properties of quantum systems in reduced dimensions, in particular the consequence due to the absence of BEC. Not only that this is the first 1D bosonic system available, it provides advantages over the mentioned 1D fermionic cases in that the present system is intrinsically clean, and moreover numerous atomic and optical tools can now be used. It may even be possible to further study novel systems such as mixtures etc.

In this paper we consider one such example, namely the absorption line shape of an optical transition in this 1D bose gas. We imagine initially a quasi-1d system of bosons in their internal ground state (referred to as the ‘a-atoms’ below). An incident optical wave excites one of the atoms to a different internal state. We are interested in the absorption line shape of this process; i.e., the probability that the absorption takes place at a given frequency $`\mathrm{\Omega }`$. We shall consider the problem at $`T=0`$, paying particular attention to the special feature due to the absence of BEC. An analogous experiment has already been done in a three dimensional trap for the 1s $``$ 2s transition in hydrogen . For a uniform 3D bose gas at $`T=0`$, the line is expected to be narrow, Lorentzian (see, e.g. , c.f. however ) with a width governed by the small ‘gaseous parameter(s)’ $`(n_{3d}a^3)^{1/2}`$. Here $`n_{3d}`$ is the particle number density in 3d and $`a`$(’s) is (are) the s-wave scattering length(s) for scattering among the bosons. For this three dimensional case, the main weight of the transition comes from exciting a particle from momentum $`\stackrel{}{p}=0`$. There is a macroscopic number of such particles. For a weakly interacting gas, the potential felt by the rest of the bosons is only slightly modified after the transition. We shall see that the situation in 1D is very different due to the absence of BEC. Exciting a single particle necessitates a substantial rearrangement of the relative motion among the rest of the particles, i.e., emission of a large number of phonons. The basic line shape is typically an edge singularity.

We shall then consider an atomic trap with tight confinement in the $`y`$-$`z`$ directions. For definiteness we shall assume that the confinement potentials in these directions are isotropic and harmonic, with frequency $`\omega _{}`$. We shall consider the case where the atoms are limited within $`0<x<L`$ but otherwise free. A weak harmonic trap potential along the x-axis has been considered in the past , and the low dimensional effects have been shown to be significantly reduced. We shall assume that such a potential is absent here. We shall further ignore the effect due to finiteness of $`L`$, a condition which we shall return below.

The bosonic system is described by a Hamiltonian containing the kinetic energy and a short-range interaction among the bosons. Anticipating that we shall eventually study the system in the quasi-1d limit, the field operator for the original bosons, referred to hereafter as the a-bosons, is expanded as $`\psi _a(\stackrel{}{r})=\psi _a(x)\chi _0(\stackrel{}{r}_{})+_j\psi _{aj}(x)\chi _j(\stackrel{}{r}_{})`$ where $`\psi _a(x)`$ is the annihilation operator for the a-bosons in the lowest subband ($`0`$) where the transverse wavefunction is $`\chi _0(\stackrel{}{r}_{})`$. $`\psi _{aj}(x)`$’s ( $`j=1,2,\mathrm{}`$) are the corresponding operators for the higher subbands (wavefunctions $`\chi _j(\stackrel{}{r}_{})`$). In the ground state $`|G>`$, $`\psi _{aj}(x)|G>=0`$ since all particles reside in the lowest subband. We shall see that $`\psi _{aj}`$ for $`j=1,2,\mathrm{}`$ will not appear anywhere below, and can pretend that the expansion of $`\psi _a(\stackrel{}{r})`$ will consist thus only the term involving the lowest subband. The extent of the wavefunction $`\chi _0`$ in the y-z direction will be denoted as $`a_{}`$. For harmonic trapping potentials in these directions $`a_{}=\sqrt{\frac{\mathrm{}}{m\omega _{}}}`$.

We shall not write down the interaction part of the Hamiltonian involving the a-bosons explicitly since we won’t need it below. For delta-function type interactions it is possible to write down all the formal many-body wavefunctions . Unfortunately manipulations of such wavefunctions are typically very mathematically involved. Moreover we shall eventually be interested in ‘final’ states where an a-atom has been excited internally. Such wavefunctions, where effectively there is a ‘foreign’ atom present, are not known (with an exception noted below). We thus proceed rather differently and follow Haldane , concentrating only on the low energy excitations which are density oscillations in the system. The field operator $`\psi _a(x)`$ ( $`=_pa_pe^{ipx}/\sqrt{L}`$) describing the motion along $`x`$ is re-expressed in terms of density $`n(x)`$ and phase $`\varphi (x)`$ by $`\psi _a(x)=[n(x)]^{1/2}e^{i\varphi (x)}`$. The effective Hamiltonian for the 1d motion can be written as

$$H_0=\frac{\mathrm{}}{2\pi }𝑑x\left[vK(\varphi )^2+\frac{v}{K}(\theta )^2\right]$$
(1)
where $`\theta `$ is related to the number density fluctuations $`\delta n(x)`$ by $`\theta =\pi \delta n`$. $`\varphi (x)`$ and $`\theta (x)`$ must be considered as operators satisfying the commutation relation $`[\varphi (x),_x^{}\theta (x^{})]=i\pi \delta (xx^{})`$. $`v`$ is the density (sound) velocity wave of the system, and $`K`$ is the (dimensionless) Luttinger liquid interaction parameter. $`vK=\pi n_o\mathrm{}/m`$ and $`K`$ is related to the compressibility of the system via $`K=\pi v\mathrm{}(n_o/\mu )`$ where $`n_o`$ is the linear number density and $`\mu `$ is the chemical potential of the system. $`K`$ depends on the interaction among the bosons. For short range interactions, $`K`$ in principle have already been obtained  but there is no general analytic form known. We simply note here that $`K=1`$ for the strongly repulsive (impenetrable) limit while $`K\mathrm{}`$ if the interaction among the bosons is weak.

$`H_0`$ can be diagonalized easily  since it is quadratic. It is useful to introduce bosonic operators $`b_q`$ and $`b_q^{}`$ which describe the sound modes of the system:

$`\theta (x)`$ $`=`$ $`i{\displaystyle \underset{q0}{}}|{\displaystyle \frac{\pi K}{2qL}}|^{1/2}sgn(q)e^{iqx}(b_q^{}+b_q)`$ (2)
$`\varphi (x)`$ $`=`$ $`i{\displaystyle \underset{q0}{}}|{\displaystyle \frac{\pi }{2qLK}}|^{1/2}e^{iqx}(b_q^{}b_q)`$ (3)
Then, apart from a constant, $`H_0=_{q0}\mathrm{}\omega _qb_q^{}b_q`$ where the mode frequencies are $`\omega _q=v|q|`$.

Now we consider the external optical field responsible for the transition in internal state. Since the atom with a new internal state is distinguishable from the original atoms, we shall call the resulting atom the ‘c-atom’. For definiteness we consider only the ‘Doppler-free’ part of the spectrum ( no momentum transfer from the exciting laser beam(s)) and assume that the perturbation responsible for the excitation is uniform throughout the 1D bose gas. The relevant part of the Hamiltonian can be then written as $`H_{\mathrm{ex}}=wd^3\stackrel{}{r}\left(\psi _c^{}(\stackrel{}{r})\psi _a(\stackrel{}{r})+\psi _a^{}(\stackrel{}{r})\psi _c(\stackrel{}{r})\right)`$ where $`\psi _c(\stackrel{}{r})`$ is the field operator for the atom in the excited internal state. We are interested in the rate of transition as a function of the excitation frequency $`\mathrm{\Omega }`$. This rate $`I(\mathrm{\Omega })`$ is given simply by the golden rule:

$$I(\mathrm{\Omega })=2\pi \underset{F}{}|<F|H_{ex}|G>|^2\delta (\mathrm{\Omega }(E_FE_G))$$
(4)
where $`|F>`$ are the set of final states. To proceed further we need to know the fate of the c-atom. There are many possibilities and we shall simply consider two extreme limits:

(1) An atom in the internal state c is not affected at all by the trapping potential for the atoms in state a. This is possible if, e.g., this potential is due to a laser with frequency near a dipole resonance for the a-atoms but far away from any of those for c. In this case the c-atom escapes from the trap and no longer interacts with the remaining a-atoms.

(2) The c-atom also feels the strong transverse confinement potential. It remains inside the trap and continues to interact with the a-atoms.

We shall treat these two cases in turn.

Case 1: In this case the final state is simply a product state, $`|F>=|\mathrm{\Phi }>\times |\stackrel{}{p}>`$, with the $`N1`$ a-atoms in state $`|\mathrm{\Phi }>`$ (with x momentum $`p_x`$) and the c-atom in the plane wave state $`|\stackrel{}{p}>`$ with momentum $`\stackrel{}{p}`$. The rate for transition into states with given $`\stackrel{}{p}`$ is

$$I_\stackrel{}{p}(\mathrm{\Omega })=2\pi |f_\stackrel{}{p}|^2\underset{\mathrm{\Phi }}{}|<\mathrm{\Phi }|a_{px}|G>|^2\delta (\mathrm{\Omega }(E_\mathrm{\Phi }+\stackrel{}{p}^2/2m+uE_G))$$
(5)
where $`f_\stackrel{}{p}d^2\stackrel{}{r}_{}e^{i\stackrel{}{p}_{}\stackrel{}{r}_{}}\chi _0(\stackrel{}{r}_{})`$ is a form factor, and $`u`$ is the difference in internal energy between the a- and c- atoms. This rate is thus proportional to the spectral function for annihilation of a particle in the 1D bose gas (hereafter $`p_xp`$ )

$$B(p,\omega ^{})\underset{\mathrm{\Phi }}{}|<\mathrm{\Phi }|a_p|G>|^2\delta (\omega ^{}(E_\mathrm{\Phi }E_{G,N1}))$$
(6)
with $`\omega ^{}=\mathrm{\Omega }+E_GE_{G,N1}\stackrel{}{p}^2/2mu`$ and where $`E_{G,N1}`$ is the ground state energy of $`N1`$ atoms in the 1D trap. $`B(p,\omega ^{})`$ is given by $`\frac{1}{\pi }\mathrm{Im}G_1(p,\omega ^{})`$ where

$$G_1(x,t)i<\psi _a^{}(0,0)\psi _a(x,t)>h(t)$$
(7)
The required Green’s function can be evaluated by finding the long distance and time behavior of $`G_1(x,t)`$. Substituting eq (3) into eq(7) one finds $`G_1\frac{1}{(x^2+v^2t^2)^{1/4K}}h(t)`$. The Fourier transform $`G_1(p,\omega ^{})`$ and hence $`B(p,\omega ^{})`$ for a given $`p`$ has the scaling form $`h(\omega ^{}v|p|)/(\omega ^2v^2p^2)^{1\alpha }`$ with $`\alpha =\frac{1}{4K}`$. Since $`1<K<\mathrm{}`$ the line shape is of the form of an edge singularity. For the three dimensional bose gas, the corresponding $`B(\stackrel{}{p},\omega ^{})`$ is proportional to a delta function at the Bogoliubov mode frequency (with an additional small incoherent part). The lack of a delta function in the present case reflects the absence of a condensate. If $`\stackrel{}{p}`$ can be independently measured this experiment would be analogous to angle resolved photoemission spectroscopy employed extensively in recent studies of high T-c cuprates .

Case 2: In this case the final result depends on the interaction between the a and c atoms, which has the general form $`g_{ac}𝑑\stackrel{}{r}\psi _c^{}(\stackrel{}{r})\psi _c(\stackrel{}{r})\psi _a^{}(\stackrel{}{r})\psi _a(\stackrel{}{r})`$. Here $`g_{ac}`$ is defined in 3D and is related to the 3D scattering length $`a_{ac}`$ by $`g_{ac}=4\pi \mathrm{}a_{ac}/m`$ . We express $`\psi _c(\stackrel{}{r})`$ also in subbands $`\psi _c(\stackrel{}{r})=_j\psi _{cj}(x)\chi _{cj}(\stackrel{}{r}_{})`$. The final states in eq (4) in general can involve terms that are off-diagonal in the band index $`j`$ for the c-atom due to the interaction between a and c. Such processes correspond to the possibility that the transverse motion of the c-atom be modified due to the interaction with the a-atoms. For a tight trap however, these contributions are small in the dilute limit. They involve the parameter $`n_og_{ac}/\pi a_{}^2\omega _{}a_{ac}/l`$. Here $`l=1/n_o`$ is the average interparticle spacing among the a-atoms. Ignoring thus these contributions, we have $`I(\mathrm{\Omega })=_j|f_j|^2I_j(\mathrm{\Omega })`$, i.e., the transition line then becomes a superposition of ‘lines’ involving transition to final states where the c-atom is within a given $`j`$-th subband, with a weight given by the form factor $`f_j=𝑑\stackrel{}{r}_{}[\chi _0(\stackrel{}{r}_{})\chi _{cj}^{}(\stackrel{}{r}_{})]`$. The shape of each of these lines is determined by $`\mathrm{Im}[\mathrm{D}_j^R(\mathrm{\Omega })]`$ where

$$\mathrm{D}_j^R(t)=i𝑑x_1𝑑x_2<\psi _a^{}(x_1,t)\psi _{cj}(x_1,t)\psi _{cj}^{}(x_2)\psi _a(x_2)>h(t)$$
(8)
which we shall now compute. The effective Hamiltonian is given by the sum of $`H_0`$ involving only the a-atoms and $`H_{cj}`$ which describes the motion of the c-atom within its $`j`$-th subband and its interaction with the a-atoms. $`H_0`$ has already been given in eq(1) and the effective $`H_{cj}`$ is given by $`H_{cj}^{(1)}+H_{cj}^{\mathrm{int}}`$ where

$$H_{cj}^{(1)}=𝑑x\left\{\frac{\mathrm{}^2}{2m}\frac{\psi _{cj}^{}(x)}{x}\frac{\psi _{cj}(x)}{x}+(ϵ_{cj}+u)\psi _{cj}^{}(x)\psi _{cj}(x)\right\}$$
(9)
consists of the ‘one-body’ (kinetic, trap and internal energy) terms and

$$H_{cj}^{\mathrm{int}}=g_{ac,jj}𝑑x\psi _{cj}^{}(x)\psi _{cj}(x)\psi _a^{}(x)\psi _a(x)$$
(10)
is due to the interaction between the a and c atoms. Here $`ϵ_{cj}`$ is the subband energy for the c-atom in its j-th subband and $`g_{ac,jj}=g_{ac}d^3\stackrel{}{r}|\chi _{cj}\chi _0(\stackrel{}{r}_{})|^2`$ depends on $`j`$. For small $`j`$ it is of order $`g_{ac}/\pi a_{}^2`$ if the trap potentials for a and c atoms are similar.

It is convenient to express $`\mathrm{D}_j^R(t)`$ in the ‘interaction’ picture:

$$\mathrm{D}_j^R(t)=i𝑑x_1𝑑x_2<\psi _a^{}(x_1,t)\psi _{cj}(x_1,t)\widehat{T}e^{i_0^tH_{cj}^{\mathrm{int}}(t^{})𝑑t^{}}\psi _{cj}^{}(x_2)\psi _a(x_2)>h(t)$$
(11)
where $`\widehat{T}`$ is the time-ordering operator and the expectation value is to be calculated with the Hamiltonion $`H_0+H_{cj}^{(1)}`$. One can rewrite eq(11) in momentum space and expand the exponent involving the interaction. $`H_{cj}^{\mathrm{int}}`$ then produces ‘scattering’ terms where the momentum of the c-atom is changed from one value to the other. The resulting calculation does not seem to be tractable analytically in general. The problem however can be simplified significantly if we assume that $`g_{ac,jj}`$ is small (condition given below) and concentrate again on the low frequency limit, i.e., just above the threshold for transition into the $`j`$-th subband (given by $`\xi _{cj}=ϵ_{cj}ϵ_0\mu +u+g_{ac,jj}n_{1d}`$ to lowest order in $`g_{ac}`$). In this case transitions must be made to states where the momenta $`p`$ of the c-atom are small. Provided that $`p<<mv`$, then the kinetic energy (of motion along x) transferred to the c-atom $`p^2/2m`$ is much less than that to the density oscillations of a-atoms $`vp`$. Thus so long as we are interested in frequency deviations $`\mathrm{\Delta }\mathrm{\Omega }`$ from the threshold which satisfy $`\mathrm{\Delta }\mathrm{\Omega }<<mv^2`$, one can ignore the kinetic energy of the motion of the c-atoms along $`x`$. For an impenetrable bose gas this amounts to limiting ourselves to $`\mathrm{\Delta }\mathrm{\Omega }<<\pi ^2\mathrm{}^2/ml^2`$. Ignoring thus the first term in $`H_{cj}^{(1)}`$ and rewriting the result back in real space, we obtain

$$\mathrm{D}_j^R(t)=iLe^{i\xi _{cj}t}<\psi _a^{}(0,t)\widehat{T}e^{i_0^t\stackrel{~}{H}_{a,int}(t^{})𝑑t^{}}\psi _a(0)>h(t)$$
(12)
where the effective interaction Hamiltonian $`\stackrel{~}{H}_{a,int}`$ acts only on the a-atoms and is given by $`g_{ac,jj}\psi _a^{}(0)\psi _a(0)`$. This result can be understood physically as follows: the external optical field produces an a $``$ cj transition at $`t=0`$ and at a general location say $`x_0`$, annihilating an a-atom while producing a c-atom there in the j-th subband. Since the $`c`$ atom moves with velocity $`p/m`$ whereas the sound waves move with the speed $`v`$, at small $`p`$ the c-atom is essentially stationary. For $`t>0`$, due to the c-atom created, a delta-function interaction at $`x_0`$ acts on the a-atoms. We are interested in the overlap between the initial state with an a-atom destroyed at $`x_0`$ and the final states with this extra interaction potential. This overlap is independent of the location of the transition $`x_0`$ which then has been set to $`0`$. The problem has thus become similar to that of X-ray absorption in solid state physics. There initially ( $`t<0`$) one has an electron gas in its ground state. At $`t=0`$ an X-ray photon is absorbed which creates a charged nuclei with an extra electron added to the electron gas, which at $`t>0`$ feels an extra local potential due to the charged nuclei. Notice however there are slight differences. Here for $`t<0`$ we have an equilibrium 1D interacting bose gas and at $`t>0`$ there is one less boson than initially.

The correlation function in (12) and thus the line shape can be found analogous to the X-ray problem . The expectation value needed can be rewritten as $`<\psi ^{}(0)e^{i\stackrel{~}{H}t}\psi (0)>`$ where $`\stackrel{~}{H}=H_0+\stackrel{~}{H}_{int}`$. $`H_0=\mathrm{}\omega _qb_q^{}b_q`$ and $`\stackrel{~}{H}_{int}=g_{ac,jj}\delta n(0)`$ can both be written in terms of $`b_q`$ and $`b_q^{}`$. The annihilation operator $`\psi e^{i\varphi }`$ acts like a displacement operator for these bosons: $`e^{i\varphi }b_qe^{i\varphi }=b_q|\frac{\pi }{2qLK}|^{1/2}`$. Thus the expectation value is the same as $`<\widehat{T}e^{i_0^t(g_{ac,jj}\pi v_N)\delta n(0,t^{})𝑑t^{}}>`$ calculated with the Hamiltonian $`H_0`$. Here we have defined $`v_Nv/K`$. The calculation can be easily done since $`H_0`$ is quadratic in the boson operators $`b_q`$ and $`b_q^{}`$. The long time dependence is given by, apart from a part oscillating sinusoidally with $`t`$ which contributes to the line shift, $`\mathrm{D}_j^R(t)1/t^{\alpha _j^{}}`$ and thus the line shape $`1/[\mathrm{\Delta }\mathrm{\Omega }]^{1\alpha _j^{}}`$ with $`\alpha _j^{}=\frac{1}{2K}[1\frac{g_{ac,jj}}{\pi v_N\mathrm{}}]^2`$. Notice that the line shape is in general different for transition to different subbands of the ‘c-atom’.

It is instructive to compare this result with that of the X-ray problem when the interaction between the nucleus after the X-ray and the electrons is modelled by a delta-function interaction of strength $`V`$. There the long time behavior of the corresponding correlation function is given by $`1/t^{\alpha _X}`$ and the line shape is $`1/[\mathrm{\Delta }\mathrm{\Omega }]^{1\alpha _X}`$ with $`\alpha _X=[1+N(0)V]^2`$, where $`N(0)`$ is the density of states at the fermi level for the electron gas . Here, the interaction strength is replaced by $`g_{ac,jj}`$ and the density of states is replaced by $`\frac{n}{\mu }=\frac{1}{\pi v_N\mathrm{}}`$. In both cases the power law line shapes (for $`\alpha _X0`$ and $`\alpha _j^{}0`$) are due to ‘orthogonal catastrophe’, that the overlap between the state just after the excitation and the new ground state of the system vanishes. There are however two differences: the $`+`$ sign in $`\alpha _X`$ is replaced by a $``$ sign since now there is one less boson at $`t>0`$ rather than one more electron in the X-ray absorption problem. The decay of $`D_j^R(t)`$ in time is slower ($`\alpha _j^{}`$ is smaller) if $`g_{ac}>0`$: near the location of the excitation, the reduction in local density of the a-atoms is also what an interaction with $`g_{ac}>0`$ prefers. The factor $`1/2K`$ arises from the fact that we have a one-dimensional interacting bose gas rather than a three-dimensional fermi liquid.

As the repulsive interaction $`g_{ac,jj}`$ between the a and c atoms increases, $`\alpha _j^{}`$ decreases and the line sharpens up. The above quantitive result for the exponents $`\alpha _j^{}`$ however holds only for small $`g_{ac,jj}`$ ( $`<<\pi v_N\mathrm{}`$ ). In the limit where interactions among all the bosons are strongly repulsive, one reaches the impenetrable limit where the bosons cannot pass each other. In this limit the statistics, i.e., the indistinguishability among the a-bosons and the distinguishability between the a- and c- bosons become irrelevant. One can check that the wavefunctions written down in  for a system of identical bosons only can be generalized to the present case if we replace one of the bosons by the foreign c-atom. It can be further verified that the absorption line shape becomes delta functions (one for each subband cj) in this impenetrable limit.

The experimental observation of the line shape discussed above should be feasible. The possibility of obtaining a quasi-1D bose gas has already been discussed in ref . The above line-shape is applicable so long as $`\mathrm{\Delta }\mathrm{\Omega }\genfrac{}{}{0pt}{}{<}{}\frac{\pi ^2\mathrm{}}{ml^2}\omega _l`$. For hydrogen atoms and with $`l1\mu m`$, $`\omega _l100\mathrm{k}\mathrm{H}\mathrm{z}`$, much larger than the present available experimental resolutions. It is sufficient for the temperature $`T`$ to be small compared with the line-width. The condition $`T<<\omega _l`$ is satisfied for $`T<\mu K`$, a temperature readily achievable. A finite length $`L`$ of the system will break the line into a set of discrete sub-lines with separations of order $`\pi ^2\mathrm{}/mL^2`$, but the shape of the ‘line’ can still be observable as long as $`L>>l`$.

I thank Ite Yu for a useful correspondence.


text:
# Full Current Statistics in Diffusive Normal-Superconductor Structures

## Abstract

We study the current statistics in normal diffusive conductors in contact with a superconductor. Using an extension of the Keldysh Green’s function method we are able to find the full distribution of charge transfers for all temperatures and voltages. For the non-Gaussian regime, we show that the equilibrium current fluctuations are enhanced by the presence of the superconductor. We predict an enhancement of the nonequilibrium current noise for temperatures below and voltages of the order of the Thouless energy $`E_{\text{Th}}=D/L^2`$. Our calculation fully accounts for the proximity effect in the normal metal and agrees with experimental data. We demonstrate that the calculation of the full current statistics is in fact simpler than a concrete calculation of the noise.

The electric current in conductors in general fluctuates. The full statistics of these fluctuations can be accessed in some cases. However, the full statistics is not easily accessible experimentally. The current experiments mainly concentrate on noise power measurements. This has become an important tool to extract information about the relevant charge transport mechanism. One can, for example, extract the effective charge of carriers responsible for the transport. Heterostructures with normal metals (N) and superconductors (S) were intensively investigated. Generally, the interest in these stems from the interplay between transport of doubly charged Cooper pairs and singly charged normal quasiparticles. The ratio between noise power in the superconducting state and in the normal state can roughly be interpreted as effective charge. Experimental indications of a doubled shot noise due to Andreev reflection in diffusive wires with one superconducting lead have been reported in Refs. . Other contributions report an enhancement of the current noise in SNS structures  with effective charges much larger that unity, possibly originating from multiple Andreev reflections.

On the theoretical side first calculations of the noise  and the full statistics  for short contacts predicted an enhanced shot noise with respect to the normal state value. For certain cases, such as tunnel junctions or disordered contacts, a doubled effective charge was predicted. The drawback of these calculations is the limitation to short contacts or small energies, since fully coherent propagation of electrons and holes is assumed and the energy dependence of the scattering amplitudes has been disregarded. For the opposite regime, where coherence of electrons and holes in the normal region plays no role, a modification of the Boltzmann-Langevin approach has been put forward recently. Interestingly, it also gives a doubling of the shot noise in the incoherent regime. Similar limitations one finds in the available calculations of the noise of SNS structures. It is tempting to “interpolate” between these two limits of coherent and incoherent propagation and to conclude, that nothing interesting happens in the intermediate regime. We will show below, that this is not the case. As analogy let us note the similarity to the linear conductance of a diffusive normal wire with one superconducting lead. Using the previously mentioned approaches one finds that the conductance has exactly the same value as in the normal state independent of temperature. Only a full calculation using the Keldysh Green’s function technique revealed, that the conductance is significantly enhanced at energies of the order of the Thouless energy $`E_{\text{Th}}`$, displaying a reentrant behaviour. This behaviour of the conductance has been experimentally verified. At present no calculation of current noise is available, which fully accounts for the peculiarities of the proximity effect.

Another recent development is to study theoretically the full statistical properties of current fluctuations. The method of choice to do so is the so called full counting statistics, pioneered in Ref. . Introducing a counting field $`\chi `$ that couples to the current operator, one can access the full current distribution. Derivatives of the current with respect to $`\chi `$ immediately generate all moments of the distribution. Thus, one circumvents the cumbersome calculation of all the moments.

In this Letter we introduce a theoretical method that allows to obtain the full current statistics for a wide class of SN-structures. We further present concrete results for the equilibrium current distribution and the nonequilibrium current noise with full account of the proximity effect. For simplicity, we restrict ourselves here to a diffusive wire mounted between a normal reservoir and a superconducting one. First we determine the equilibrium current distribution. It turns out that the large current fluctuations are more probable in the superconducting case. Second we calculate the temperature and voltage dependent current noise. We find that the noise is enhanced (similar to the conductance) at voltages of the order of and for temperatures below the Thouless energy. In the respective limits our results agree with previous calculation.

Let us first present the theoretical framework. We use a recently developed extension  of the Keldysh technique to compute the statistics of our proximity structure. There it was shown that the current statistics can be obtained by imposing the modified boundary condition on one reservoir’s(L) Green’s function

$$\stackrel{ˇ}{G}_\text{L}(\chi )=e^{\frac{i}{2}\chi \stackrel{ˇ}{\tau }_\mathrm{K}}\stackrel{ˇ}{G}_\text{L}e^{\frac{i}{2}\chi \stackrel{ˇ}{\tau }_\mathrm{K}}.$$
(1)
Here $`\stackrel{ˇ}{G}_\text{L}`$ is the standard Green’s function of a reservoir and $`\stackrel{ˇ}{\tau }_\mathrm{K}=\widehat{\sigma }_1\overline{\tau }_3`$ is a matrix in Keldysh($`\widehat{}`$)-Nambu($`\overline{}`$) space. Statistical properties are encoded in the dependence on the counting field $`\chi `$. The current probability distribution is found from

$$P(I)=_\pi ^\pi 𝑑\chi e^{S(\chi )i\chi It/e}.$$
(2)
Here $`t`$ denotes the time of observation. The action $`S(\chi )`$ can be found from integration of the $`\chi `$-dependent current

$$i\frac{e}{t}\frac{S(\chi )}{\chi }=I(\chi )=\frac{1}{8e}𝑑E\text{Tr}\left[\stackrel{ˇ}{\tau }_K\stackrel{ˇ}{I}(\chi )\right].$$
(3)
The spectral matrix current in a diffusive wire is given by

$$\stackrel{ˇ}{I}(\chi )=\sigma \stackrel{ˇ}{G}(\chi )\frac{}{x}\stackrel{ˇ}{G}(\chi ),$$
(4)
where $`\sigma `$ is the conductivity. This matrix current has to be found by defining an appropriate circuit. The manual on how to do this can be found in . Note, that the only change in comparison to the calculation of the conductance is the modified boundary condition (1). All other relations defining the actual circuit remain unchanged, as long as they are respecting the full matrix structure.

We stress that our approach is not in contradiction with the general scattering matrix approach of Ref. . If one would know the electron-hole scattering amplitudes for the system under consideration, and would not disregard their energy dependence, one could obtain the same result. The characterization of these amplitudes would have to be performed along the lines of Ref. . We also emphasize, that the calculation of the full current statistics can be done in a finer, simpler and more compact way than a separate calculation of it’s second order perturbations series, i. e. the noise.

Let us now specify our system. A diffusive metal is connected to a normal terminal on one end and to a superconducting terminal at the other end. Inside the mesoscopic wire the quasiclassical transport equations are obeyed. In the normal metal they read

$$D\frac{}{x}\left(\stackrel{ˇ}{G}(x,\chi )\frac{}{x}\stackrel{ˇ}{G}(x,\chi )\right)=[iE\overline{\tau }_3,\stackrel{ˇ}{G}(x,\chi )].$$
(5)
Here $`D`$ is the diffusion constant and $`x`$ the coordinate along the wire, which has a length $`L`$. It’s conductance is $`G_\text{N}=\sigma A/L`$ (cross section $`A`$). At both ends boundary conditions to reservoirs have to be supplied. At the normal end with ideal connection the Green’s function is continuous: $`\stackrel{ˇ}{G}(0,\chi )=\stackrel{ˇ}{G}_\text{L}(\chi )`$. The other end is connected to a superconducting reservoir by a contact of negligible resistance, which leads to the boundary condition $`\stackrel{ˇ}{G}(L,\chi )=\stackrel{ˇ}{G}_\text{R}`$. A circuit representation of the system is depicted in the inset of Fig. 1.

In a normal reservoir (which we will consider in the rest of the paper) $`\stackrel{ˇ}{G}_\text{L}`$ is given by

$$\stackrel{ˇ}{G}_\text{L}=\left(\begin{array}{cc}\overline{\tau }_3\hfill & \hfill \overline{K}\\ 0\hfill & \hfill \overline{\tau }_3\end{array}\right),\overline{K}=2\left(\begin{array}{c}12f(E)0\hfill \\ 012f(E)\hfill \end{array}\right).$$
(6)
The distribution function at voltage $`V`$ and temperature $`T`$ is given by $`f(E)=(\mathrm{exp}((E+eV)/T)+1)^1`$. A superconducting reservoir at zero voltage is described by

$$\stackrel{ˇ}{G}_\text{R}=\left(\begin{array}{cc}\overline{R}& (\overline{R}\overline{A})\mathrm{tanh}(E/2T)\\ 0& \overline{A}\end{array}\right).$$
(7)
Advanced and retarded Green’s functions in (7) possess the structure $`\overline{R}(\overline{A})=g_{\text{R,A}}\overline{\tau }_3+f_{\text{R,A}}\overline{\tau }_1`$, $`f_{\text{R,A}}=i\mathrm{\Delta }/((E\pm i\delta )^2\mathrm{\Delta }^2)^{1/2}`$ and $`g_{\text{R,A}}`$ following from the normalization condition $`f^2+g^2=1`$ in a standard BCS-superconductor.

Let us briefly comment on the numerical procedure of the solution. It is most convenient to solve the matrix equation (5) directly. For this purpose the diffusive wire is represented by a discrete set of $`n`$ nodes, each represented by a Green’s function $`\stackrel{ˇ}{G}_k`$ connected in series by tunnel junctions of conductance $`g=(n+1)G_\text{N}`$. The matrix current between two neighboring nodes is then given by $`\stackrel{ˇ}{I}_{k,k+1}=\frac{g}{2}[\stackrel{ˇ}{G}_k,\stackrel{ˇ}{G}_{k+1}]`$. The right hand side of (5) has a similar form and is represented as a ’decoherence’ terminal (with a Green’s function $`\stackrel{ˇ}{G}_{\text{dec}}=i2(E/nE_{\text{Th}})\overline{\tau }_3`$) connected to each node. The matrix current conservation for node $`k`$ follows from discretizing Eq. (5) and reads $`[g(\stackrel{ˇ}{G}_{k+1}+\stackrel{ˇ}{G}_{k1})+\stackrel{ˇ}{G}_{\text{dec}},\stackrel{ˇ}{G}_k]=0`$. The resulting set of equations for the nodes Green’s functions is then solved by iteration.

Because of its importance in the following discussion, we have depicted the differential resistance of the proximity wire in the inset of Fig.2. It displays the famous reentrance behaviour , i.e. starting from the normal state resistance at zero voltage a minimum at energies of the order of several $`E_{\text{Th}}`$ occurs. Above the minimum it decays slowly to the normal state value as $`(E_{\text{Th}}/E)^{1/2}`$. Therefore, the Thouless energy is the central quantity in the physics of the proximity effect. The disagreement of the theory and the experimental data from Kozhevnikov et al.  may possibly result from heating effects, not accounted for in the theoretical calculation.

As a first application, we study the distribution of current fluctuations in equilibrium. We therefore put $`V=0`$ and find the solution of the above equations for different values of $`\chi `$. We evaluate the integral over $`\chi `$ in (2) in the saddle point approximation, we take $`\chi `$ as complex and expand the exponent in (2) around maximum. The integral yields then $`P(I)\mathrm{exp}(S(\chi )\chi It/e)`$, which we plot implicitly as a function of $`I(\chi )`$. To extract generic fluctuation properties of the proximity effect we set here $`\mathrm{\Delta }T,E_{\text{Th}}`$.

Results of this calculation are displayed in Fig. 1. The current $`I`$ is normalized by $`G_\text{N}T/e`$ and $`\mathrm{ln}P(I)`$ is plotted in units of $`G_\text{N}Tt/e^2`$. The solid line shows the distribution in the normal state, which does not depend on the Thouless energy. In our units, this curve is consequently independent on temperature. In the superconducting state the Thouless energy does matter, and the distributions depend on the ratio $`E_{\text{Th}}/T`$. We observe that large fluctuations of the current in the superconducting case are enhanced in comparison to the normal case, and in both cases are enhanced in comparison to Gaussian noise. For comparison we plotted the Gaussian distribution $`\mathrm{exp}(tI^2/4G_\text{N}T)`$ of the current measurements, owing to the fact that the conductance is the same in all cases. The differences between the normal and the superconducting state occur in the regime of non-Gaussian fluctuations.

Let us now turn to the nonequilibrium current noise as second application. Previous results can be summarized in two statements. At temperatures larger than an applied voltage the wire displays the usual thermal noise $`S_\text{I}(TeV)=4G_\text{N}T`$ in accordance with the fluctuation-dissipation theorem. Note that in a calculation, which neglects the proximity effect, the conductance of the normal wire is equal to the normal state conductance $`G_\text{N}`$ independent of voltage or temperature. For large voltages $`eVT`$, on the other hand, the shot noise $`S_\text{I}(TV)=(4/3)eG_\text{N}V`$ is doubled in comparison to the normal state value $`(2/3)eG_\text{N}V`$. As discussed previously these results do not depend on coherence between electrons and holes. Consequently the Thouless energy plays no role. However, this can not be true, since the conductance is enhanced by approximately $`10\%`$ at energies of the order of $`E_{\text{Th}}`$. In the following we will show, that the proximity effect indeed changes the nonequilibrium noise of the system.

Using the method developed in , we are able to take the coherence between electrons and holes fully into account. We have to solve the Usadel equation (5) taking into account the boundary condition (1). Then, the current noise power is given by

$$S_\text{I}=2ei\frac{I(\chi )}{\chi }|_{\chi =0}.$$
(8)
There are two ways to attack this problem. One way to determine the noise is to expand the Green’s functions and Eq. (5) to first order in $`\chi `$ and thus to obtain an equation for the noise. The other way is to solve the full matrix equations, and perform the differentiation in (8) numerically. Below we will use the second way to find the nonequilibrium noise. Nevertheless, let us sketch the derivation of an equation for the noise. We define $`\stackrel{ˇ}{G}(x,\chi )=\stackrel{ˇ}{G}_0(x)i(\chi /2)\stackrel{ˇ}{G}_1(x)`$ and $`\stackrel{ˇ}{I}(x,\chi )=\stackrel{ˇ}{I}_0(x)i(\chi /2)\stackrel{ˇ}{I}_1(x)`$. As result we find

$`\stackrel{ˇ}{I}_1(x)`$ $`=`$ $`\sigma \left(\stackrel{ˇ}{G}_0(x){\displaystyle \frac{}{x}}\stackrel{ˇ}{G}_1(x)+\stackrel{ˇ}{G}_1(x){\displaystyle \frac{}{x}}\stackrel{ˇ}{G}_0(x)\right),`$ (9)
$`{\displaystyle \frac{D}{\sigma }}{\displaystyle \frac{}{x}}\stackrel{ˇ}{I}_1(x)`$ $`=`$ $`[iE\overline{\tau }_3,\stackrel{ˇ}{G}_1(x)].`$ (10)
From this equation the generalization of the Boltzmann-Langevin equation to superconductors can in principle be derived. The boundary conditions at the reservoirs read $`\stackrel{ˇ}{G}_1(0)=[\stackrel{ˇ}{\tau }_\text{K},\stackrel{ˇ}{G}_\text{L}]`$ at the left end and $`\stackrel{ˇ}{G}_1(L)=0`$ at the right end. Finally the noise is $`S_\text{I}=e𝑑E\text{Tr}\stackrel{ˇ}{\tau }_\text{K}\stackrel{ˇ}{I}_1(x)`$. By taking the trace of Eq. (10) multiplied with $`\stackrel{ˇ}{\tau }_\text{K}`$ it follows that it does not matter, where the noise is evaluated, which is as it should be.

We turn now to concrete results for the noise power. The influence of the proximity effect is most easily seen in the differential noise $`dS/dV`$. Numerical results for different temperatures are displayed in Fig. 2. The inset shows the differential conductance for the same parameters. The differential noise shows a remarkable enhancement at energies of the order of the Thouless energy. Following a linear increase at low voltage the differential noise overshoots the doubled normal differential noise, which is recovered at large voltages. The maximal differential noise occurs, if the voltage is of the order of $`4E_{\text{Th}}`$. The largest enhancement is found for $`TE_{\text{Th}}`$ and is roughly $`10\%`$. At higher temperature the differential noise approaches the Boltzmann-Langevin result , shown as grey lines. At zero temperature the reentrant behaviour makes the connection to the result obtained within random matrix theory. The nontrivial behaviour in the regime between these two approaches shows the importance of phase coherence.

We can compare our results with experimental data on the noise power obtained by Kozhevnikov et al..  The experimental parameters correspond to $`T/E_{\text{Th}}=0.2`$. Given that our approach contains no adjustable parameters, the agreement is very good. A possible explanation for the difference between theory and experiment is that heating effects in the experiments may be important. These have been completely disregarded in the theoretical calculation. This would also explain the smaller reentrance of the differential resistance seen in the experiment. Note, however, that the energy scale, at which the influence of the proximity effect is seen, is unchanged.

In conclusion, we have developed a method to calculate statistical properties (with emphasis on the current noise power) of normal-metal–superconductor heterostructures. The method is embedded in a matrix circuit theory , which allows to find the full current statistics of a wide class of systems. We applied the method to a normal diffusive wire with one normal and one superconducting reservoir. In equilibrium it turns out that the current fluctuation distribution differs from that of a purely normal system in the non-Gaussian regime. Large fluctuations of the current are enhanced by the proximity effect. For temperatures below $`E_{\text{Th}}`$ we found that the nonequilibrium current noise shows a reentrant behaviour with a maximum for voltages of the order of $`E_{\text{Th}}`$. This is in qualitative agreement with experiment.

We thank A. A. Kozhevnikov and R. J. Schoelkopf for discussions and for providing us with the original data from Ref. . W.B. was financially supported by the “Stichting voor Fundamenteel Onderzoek der Materie” (FOM) and the “Alexander von Humboldt-Stiftung”.


text:
# 1 Calorimeters geometry considered

## 1 Calorimeters geometry considered

The pion responses are obtained with GEANT simulations with a detailed description of CMS calorimeter geometry (CMSIM 115), version TDR-2/1/.

The ECAL is the PbWO4 one (readout number 1).

An additional 1 cm scintillator layer is placed in front of HB to compensate for the energy loss in the cables, electronics and cooling system on the back side of the ECAL, which is simulated by a uniform slab of 0.2 $`\lambda `$ of dead material (readout number 2).

The sampling thicknesses of copper alloy (90%Cu + 10%Zn) are 5 cm in the barrel segments, except the inner and outer plates which are 7 cm stainless steel. The 9 mm gaps between the absorber plates are filed with 4 mm scintillator planes (readout number 3).

To improve the energy measurement for late developing hadron showers, tail catcher layers are inserted in front and behind the first return yoke iron layer (RY1) (readout number 4).

The regions behind the ECAL are empty except for the scintillator layer and equivalent dead material uniform absorber layer. The effects of the CMS 4T field are fully included.

## 2 Data sample

We have performed a GEANT simulation of the response to pions for a set of incident energies from 10 to 3000 GeV at pseudorapidity set to 0. GHEISHA was used as a hadron interaction simulator. The sum of the energies deposited by the showers in the active elements of each readout has been stored in a disk file so that the showers could be analysed later.

## 3 Standard calibration

The standard calibration used just sampling coefficients for calorimeter parts with different sampling ratio. The reconstructed energy of simple shower is given by the weighted sum of the energies deposited in the readouts:

$$E^{rec}=\underset{i=1,4}{}C_iE_i,$$

where:

$`E_i`$ \- amplitude of the signal from the calorimeter longitudinal segmentation (readouts);

$`C_i`$ \- calibration coefficients.

Coefficients $`C_i`$, i = 1, 4 are determined by the minimisation of the width of the energy distributions.

The Gaussian part of reconstructed energy distributions at various incident energies are then fitted to obtain the calorimeter energy resolution. The energy resolution is parametrised by the expression:

$$\sigma /E=a/\sqrt{E}b.$$

The obtained energy resolution by the standard calibration is shown on fig. 1 and the residuals of the reconstructed energy on fig. 2.

## 4 Non-linear technique

The nonlinear behaviour of the calorimeter response could be solved by application of the non-linear method which improve the linearity of the calorimeter response and the energy resolution in the broad energy range. Non-linear technique is the selection of some additional parameters which provide correct energy reconstruction of the hadron shower.

Thus the reconstructed energy $`E_{nl}^{rec}`$ is parametrised as:

$$E_{nl}^{rec}=\underset{i=1,4}{}f_i(\stackrel{}{A},E_i)E_i,$$

where $`f_i`$ are non-linear functions of unknown parameters $`\stackrel{}{A}`$.

Large number of free parameters (more than 50) are needed to better describe the non-linearity. For such a big number of parameters minimisation performed by MINUIT is inefficient. To this end we used the specialised program REGN based on a autoregularized Newton type method/2/. In the REGN computer code/3/ the $`\chi ^2`$ is one of the criteria available for solving the system and for testing the mathematical model. Other criteria permit one to chose uniquely between several model functions the best one/4,5,6/.

The energy resolution obtained applying non-linear calibration method is shown in fig. 1 together with results from the standard calibration. In Figure 2 the residuals of the reconstructed energies are shown for both methods. The comparison of the standard calibration method to the non-linear technique shows clear improvement in the resolution and linearity of the reconstructed energy.

## 5 Summary

The CMS calorimeters pion calibration was done using two different approaches. Calibration using non-linear technique improves the resolution and linearity of the reconstructed energy.