text:
# X-ray Binaries and Globular Clusters in Elliptical Galaxies

## 1 Introduction

The X-ray emission from normal elliptical galaxies has two major components: soft ($`kT0.21`$ keV) emission from diffuse gas and harder ($`kT6`$ keV) emission from populations of accreting (low-mass) stellar X-ray binaries (LMXB). The X-ray properties of the LMXB component have been difficult to determine, due to their spatial confusion with diffuse gaseous emission and to spectral hardness which places much of the LMXB emission outside the effective bandpasses of most imaging X-ray satellite spectrometers. The presence of the LMXB component has been inferred in part through observations of spectral hardening in ellipticals with progressively smaller X-ray to optical luminosity ratios (Kim, Fabbiano, & Trinchieri 1992), indicating that they have relatively little gas, exposing the harder LMXB component. Populations of LMXBs are also expected in ellipticals simply by analogy with detections of discrete LMXB sources in nearby spheroids such as the bulges and halos of our Galaxy and M31 (Forman, Jones & Tucker 1985; Canizares, Fabbiano & Trinchieri 1987), as well as in the radio galaxy Centaurus A (Turner et al. 1997).

A simple argument suggests that the total X-ray luminosities of LMXB populations in ellipticals might be proportional to the stellar luminosities of the galaxies: if the properties of low mass binary stellar systems (such as the fraction of stars in binaries, the distributions of binary separations and mass ratios, etc.) are largely independent of their galactic context, the number of LMXBs (hence their aggregate X-ray luminosity) should be simply proportional to the number of stars in the galaxy (and thus their total optical luminosity).

High angular resolution Chandra observations are now allowing individual LMXBs to be resolved out of the diffuse gaseous X-ray emission in nearby ellipticals (Kraft et al. 2000, 2001; Sarazin, Irwin, & Bregman 2000, 2001; Angelini, Loewenstein, & Mushotzky 2001; Finoguenov & Jones 2001), which makes their composite spectral analysis much easier. The bulk of the hard emission in normal ellipticals indeed comes from LMXBs, rather than from advection-dominated accretion flows onto massive, central black holes, proposed by Allen, di Matteo & Fabian (2000).

Until Chandra observes more nearby ellipticals, the strongest spectral constraints to date on the hard stellar LMXB component in a large sample of ellipticals will still come from ASCA spectra (Matsumoto et al. 1997; White 2000, 2001, 2002). Since the effective angular resolution of ASCA imaging spectrometers is $`23^{}`$ (half-power diameter), confusion prevents individual LMXBs from being easily resolved out of the diffuse gas in ellipticals. The hard LMXB component can be spectrally distinguished from the softer gaseous component, however. Matsumoto et al. (1997) separated the hard LMXB component in ellipticals from softer gaseous emission by considering ASCA GIS spectral energies above 4 keV. Stacking 4-10 keV GIS spectra for 12 ellipticals, they found a best-fit thermal bremsstrahlung model with $`kT=12_{5.5}^{+29.3}`$ keV (where errors are 90% confidence limits); a power-law with photon index $`\alpha =1.8_{0.4}^{+0.4}`$ fit equally well. Matsumoto et al. (1997) found that the X-ray luminosities of the LMXB components were proportional to the optical luminosities of the ellipticals, but with a surprisingly large scatter (a factor of $`4`$). Some of the ellipticals included in the sample have significant X-ray emission from active galactic nuclei (AGN), which may account for some of the scatter.

White (2000, 2001, 2002) performed a similar ASCA analysis, but with a larger spectral bandwidth (0.7-10 keV), on six normal ellipticals (i.e., ellipticals without significant AGN emission). Spectra were extracted from within metric radii of 6$`r_e`$ from the galactic centers, where $`r_e`$ is the optical effective radius of a galaxy. The GIS spectra of the six ellipticals were simultaneously fit with both soft (gaseous) and hard (LMXB) emission components. (Only GIS data were used because the GIS detectors have twice the effective area of the SIS detectors above 7 keV.) The temperatures (or power-law indices) of the hard components in the galaxies were tied together, while the temperatures of the soft components (if present) were allowed to vary individually. Much tighter spectral model constraints were provided by the increased spectral bandwidth compared to the 4-10 keV bandwidth in the stacked spectral study of Matsumoto et al. (1997). The spectra of the LMXB components were fit equally well by a bremsstrahlung model with $`kT=6.4_{1.1}^{+1.7}`$ keV or a power-law model with photon index $`\alpha =1.82_{0.09}^{+0.10}`$ (errors are 90% confidence limits). Individual fits to each galaxy in the set were consistent with the results of the joint fits and fluxes were obtained by adopting the best jointly fit temperature. These are the tightest constraints to date on the global spectral properties of the stellar LMXB component in ellipticals. GIS X-ray fluxes were determined for the LMXB components in an additional 8 ellipticals which had poorer photon statistics by fixing the temperature of the hard component to 6.4 keV. The resulting X-ray luminosities of the LMXB components in the 14 galaxies were found to be proportional to the optical luminosities of the galaxies, with a factor of $`410`$ scatter.

Recent Chandra observations (Sarazin et al. 2000, 2001: Angelini et al. 2001; Blanton, Sarazin, & Irwin 2001) also show that the X-ray luminosities of the resolved LMXB components in ellipticals exhibit significant scatter at a given optical luminosity. (The composite spectra of these LMXB components are also consistent with the hard ASCA spectral results described above.) Although their scatter in luminosity is much smaller than that of the softer gaseous component (which ranges a factor of 100 in X-ray/optical flux ratio), it is still larger than expected if the LMXB component is strictly proportional to the stellar component. Following a preliminary study by White (2000, 2001), we next consider two possible sources of variance in the X-ray luminosities of LMXB populations in elliptical galaxies: age differences among ellipticals and differences in their globular cluster populations.

## 2 Possible Sources of Variance in LMXB Populations in Ellipticals

### 2.1 Variations in Galaxy Ages

If LMXB populations are a subset of the primordial binary systems in a galactic stellar field, then the evolution of an LMXB population should be tied to the evolution of visible stars. The evolution of the aggregate X-ray luminosity of an LMXB population is best assessed with a population synthesis calculation (Wu 2001). A simpler, more phenomenological approach was adopted by White & Ghosh (1998). The luminosity of an LMXB population (neglecting burst sources) is expected to rise for about a Gyr (White & Ghosh 1998), driven by the time it takes the less massive secondaries in binaries to evolve off the main sequence, overflow their Roche lobes, and start dumping mass onto the more massive (primary) compact stellar remnants in the binaries. After this Gyr ramp-up, the luminosity of an LMXB population is expected to slowly decline (White & Ghosh 1998; Wu 2001).

The optical stellar luminosities of galaxies are also slowly declining with time. Figure 1 of Bruzual & Charlot (1993) shows that the (visual) optical luminosity of an elliptical galaxy declines as $`L_Vt^{0.95}`$. This time dependence is close enough to estimates for the temporal decline in the X-ray luminosity of LMXB populations that it is difficult to assess how the X-ray to optical ratio of the LMXB component is evolving at the present epoch.

For our purposes, ambiguity in the theoretical expectation for the temporal evolution of $`L_{\mathrm{LMXB}}/L_{\mathrm{opt}}`$ does not matter — we can deduce the evolution empirically. If variations in the X-ray properties of the LMXB component for galaxies of a given optical luminosity are correlated with independent age estimates, the $`L_{\mathrm{LMXB}}/L_{\mathrm{opt}}`$ ratio itself might then be used as an indicator of galactic age.

To test this possibility, we compare the $`L_{\mathrm{LMXB}}/L_{\mathrm{opt}}`$ ratio for the LMXB components from the elliptical sample of White (2002) to recent stellar age estimates for elliptical galaxies (Trager et al. 2000). The stellar age estimates were derived by comparing the strengths of optical stellar spectral line indices of $`<`$Fe$`>`$, Mg b and H$`\beta `$ (which were corrected for the effects of non-solar abundance ratios), to the single-burst stellar population models of Worthey (1994). For their sample of 50 ellipticals, Trager et al. (2000) find ages ranging from $`1.518`$ Gyr. In Figure 2.1, the relevant optical stellar age estimates (ranging from $`414`$ Gyr) are compared to the $`L_{\mathrm{LMXB}}/L_{\mathrm{opt}}`$ ratios for the LMXB populations in the eight ellipticals of (the fourteen in) White (2002) that were included in the Trager et al. (2000) sample. The X-ray luminosities were calculated in the 0.5-4.5 keV range and their errors are 90% confidence limits. There is apparently no correlation, so variations in the $`L_{\mathrm{LMXB}}/L_{\mathrm{opt}}`$ ratios for LMXB populations are not likely to be due to galactic age differences.

### 2.2 Globular Cluster Population Variations

Alternatively, most LMXBs may be formed in globular clusters, rather than in the stellar field of a galaxy. Observations of LMXBs in our own Galaxy show that it is very difficult to make LMXBs in the field. There are comparable numbers of high mass X-ray binaries (HMXBs) and LMXBs, but HMXB lifetimes are only $`10^5`$ yr, which are $`10^{25}`$ times shorter than LMXB lifetimes. Thus, the formation of LMXBs in our Galaxy is $`10^{25}`$ times less efficient than the formation of HMXBs.

Meanwhile, globular clusters are the most efficient formation site for LMXBs in our Galaxy: $`20\%`$ of all LMXBs in our Galaxy reside in globular clusters, yet globular clusters contain $`0.1\%`$ of the stars in our Galaxy (Katz 1975). Thus, globular clusters form LMXBs $`>200`$ times more efficiently than the rest of the Galaxy, through tidal capture during close stellar encounters (Clark 1975).

Given this, it is reasonable to consider whether $`most`$ LMXBs are formed in globular clusters. Grindlay (1984) suggested that X-ray bursting binaries in the bulge of our Galaxy (not residing in globular clusters) were formed in globular clusters which were later destroyed by tidal shocks while passing through the Galaxy’s disk. This disruption mechanism is not relevant to globulars residing in elliptical galaxies, but if globulars are on eccentric orbits, they may still lose stars due to galactic tides at perigalacticon. Vesperini (2000) estimates that only $`10\%`$ of globular clusters are disrupted in the most luminous ellipticals, but the fraction disrupted increases considerably for progressively smaller elliptical luminosities. Also, some LMXBs formed in globular clusters may be ejected by the velocity kicks imparted by supernovae during the formation of neutron stars. Typical supernovae kick velocities are thought to range from $`100500`$ km/s (Kalogera 1996; Terman, Taam, & Savage 1996; Fryer, Burrows, & Benz 1998), while the escape velocities of globular clusters range from $`560`$ km/s (Webbink 1985). Since LMXB lifetimes can be quite long ($`10^9`$ yr), LMXBs ejected from globular clusters can be visible outside the clusters for a long time. If LMXBs are made primarily in globular clusters, we would expect their number to be proportional to the number of globular clusters, regardless of the fraction ejected from clusters.

With the excellent angular resolution of Chandra, it is possible to determine positions of LMXBs with sufficient accuracy to determine if they coincide with globular clusters. Observations of ellipticals show that a significant fraction are located in globulars (Sarazin et al. 2000, 2001; Angelini et al. 2001; Randall, Sarazin, & Irwin 2001; White, Davis & Hanes 2002). The fraction ranges from $``$20 to 70%. These are lower limits, since deeper Chandra exposures may detect fainter X-ray sources in additional globulars. In the same region of these galaxies, the globular clusters typically provide only $`0.1`$% of the optical light, which implies that optical stars in globular clusters are much more likely to be the donor stars in LMXBs than in the field. The primary limitation with identifying LMXBs with globulars in ellipticals at present is the lack of Hubble Space Telescope (HST) globular cluster lists for many of the Chandra observed ellipticals.

Since elliptical galaxies exhibit an order of magnitude range of globular cluster populations for a given galaxian luminosity, we can assess whether these variations are correlated with variations in $`L_{\mathrm{LMXB}}/L_{\mathrm{opt}}`$. In Figure 2.2, we compare the specific frequency of globular clusters $`S`$ (the number of globular clusters per galaxy visual luminosity in units of the luminosity of a $`M_V=15`$ galaxy) from the compilation of Kissler-Patig (1997) to the $`L_{\mathrm{LMXB}}/L_{\mathrm{opt}}`$ ratio for the LMXB populations in the elliptical sample described above. The Kissler-Patig (1997) compilation includes error estimates for $`S`$, as indicated in the figure, for all but one elliptical in this sample. There appears to be a strong correlation (White 2000, 2001), and the relationship is consistent with direct proportionality: $`L_{\mathrm{LMXB}}/L_{\mathrm{opt}}S^{1.2\pm 0.4}`$. Thus, LMXB populations may indeed be controlled by globular cluster populations.

## 3 Discussion and Conclusions

This study suggests that globular cluster populations in ellipticals control the evolution of their populations of low mass stellar X-ray binaries. At the present epoch, age differences among the ellipticals in the small sample described above do not seem to be correlated with the variations in their X-ray binary populations. It would be useful to compare the LMXBs in early-type galaxies with more detailed properties (such as color and magnitude) of their globulars (other than just the number), to see if the history of star formation in the globular affects the LMXB population. The Chandra observations of NGC 1399 suggested that LMXBs are more likely to be found in redder globular clusters (Angelini et al. 2001). There also is evidence that LMXBs are more likely to be found in brighter globulars (Angelini et al. 2001), although it is unclear whether this is just due to the larger number of stars in brighter globulars, or indicates a higher probability per star in more massive clusters.

It is possible that most or all LMXBs in early-type galaxies are formed in globular clusters. In this scenario, any field LMXBs in ellipticals would have escaped from globulars. The field LMXBs might have been ejected by the kick velocities resulting from supernovae, by stellar dynamical processes such as superelastic encounters in which internal binding energy in the binaries is converted into kinetic energy of motion of the binary, or by the dissolution of the globular due to tidal effects.

There is some evidence that the LMXBs in globular clusters in elliptical galaxies are brighter in X-rays than those in the field (Angelini et al. 2001). If all of the LMXBs were formed in globular clusters, this would require that the less luminous and presumably less massive LMXBs systems were preferentially ejected from globular clusters. Alternatively, if a significant fraction of the field LMXBs were formed in the field, this would imply differing stellar and binary evolution histories in the field and in globular clusters.

Many sources associated with globulars have X-ray luminosities which exceed the Eddington luminosity of a 1.4 $`M_{}`$ neutron star, $`L_{\mathrm{Edd},\mathrm{NS}}2\times 10^{38}`$ ergs s<sup>-1</sup> (Sarazin et al. 2000, 2001; Angelini et al. 2001; Blanton et al. 2001; Randall et al. 2001; White et al. 2002). In general, LMXB luminosity functions within early-type galaxies have a break at a luminosity very close to $`L_{\mathrm{Edd},\mathrm{NS}}`$ (Sarazin et al. 2000, 2001; Blanton et al. 2001; Randall et al. 2001). This break may separate accreting black holes, at higher luminosities, from accreting neutron stars at lower luminosities (Sarazin et al. 2000). If the sources are Eddington-limited single binaries, then the brightest sources associated with globular clusters must contain quite massive ($`20M_{}`$) black holes (Sarazin et al. 2000). It is not clear how globular clusters would form and retain such massive black holes in LMXBs and also not have them form binary black hole systems (e.g., Portegies Zwart & McMillan 2000).

On the other hand, these super-Eddington sources might be globular clusters containing multiple X-ray sources (Angelini et al. 2001). The difficulty with this suggestion is that the fraction of globulars which contain at least one X-ray source is low; the fraction is about 3-4% in NGC 1399 (Angelini et al. 2001), NGC 4697 (Sarazin et al. 2000, 2001), and in NGC 4649. (It is interesting that this fraction is almost the same in these three ellipticals.) The luminosity function of the LMXBs isn’t steep enough that multiple sources should affect the luminosity function at the high end, given this low fraction. On the other hand, if the LMXBs are mainly confined to a smaller subset of globular clusters (e.g., very bright and/or very red and/or core–collapsed systems), then the fraction of the relevant subset of globulars with X-ray sources and the probability of multiple sources might be much higher.

There is some evidence that these luminous super-Eddington X-ray sources are associated with globular clusters in elliptical and S0 galaxies, but not with the bulges of spiral galaxies. There are no super-Eddington LMXBs associated with globular clusters in the bulge of our Galaxy (e.g., Hut et al. 1992). The same is true of the central 5′ of the bulge of M31 (Shirey et al. 2001). The luminosity of LMXBs in the nearby large spiral bulge NGC 1291 has a cut-off at a luminosity which is very close to $`L_{\mathrm{Edd},\mathrm{NS}}`$. Unless this is the result of statistical fluctuations and the large populations of globular clusters in elliptical galaxies, this suggests that the binaries in globulars in ellipticals have different histories than those in spiral bulges.

In general, the specific frequency of globular clusters $`S`$ increases along the Hubble sequence from late to early-type spirals, from spirals to S0s, from S0s to ellipticals, and from normal giant ellipticals to cDs (e.g., Harris 1991). If LMXBs all are born in globular clusters, then one would expect the specific frequency of LMXBs and their X-ray to optical ratio to increase along the same sequence. In particular, a number of central dominant ellipticals in galaxy clusters are notable for having enormous populations of globular clusters. White (1987) suggested that these unusually populous systems belong to the galaxy clusters, rather than individual galaxies: globular clusters will have been tidally stripped from individual galaxies during the collapse of a galaxy cluster and will then virialize along with other cluster constituents. If LMXBs are formed primarily in globular clusters, then central dominant galaxies with high specific globular cluster frequencies (such as M87, NGC 1399, NGC 3311 and NGC 4874) should be particularly rich in LMXBs, as well. Chandra observations have already confirmed the large populations of LMXBs in NGC 1399 (Angelini et al. 2001), and in the region of this galaxy with a globular cluster list based on HST observations, 70% of the LMXBs are located in globular clusters.

R. E. W. was supported in part by a National Research Council Senior Research Associateship at NASA/Goddard and by the National Aeronautics and Space Administration through $`Chandra`$ Award Numbers GO 0-1144X and GO 1-2095X issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-39073. C. L. S. was supported in part by $`Chandra`$ Awards GO0-1141X and GO1-2078X and by XMM/Newton award NAG5-10074.


format: text

text:
# Bounding The Degrees of Generators of a Homogeneous Dimension 2 Toric Ideal

## Introduction

Let $`𝒜=(a_{ij})`$ be a $`d\times n`$ matrix of integers. Let $`k`$ be an arbitrary ground field. Let $`R=k[z_1,\mathrm{},z_n]`$ and $`S=k[y_1,\mathrm{},y_d]`$. From $`𝒜`$ we get a ring homomorphism $`\psi `$ from $`R`$ to $`S`$ by sending $`z_j`$ to $`y_1^{a_{1j}}y_2^{a_{2j}}\mathrm{}y_d^{a_{dj}}`$. Let $`I`$ be the kernel of this map. Ideals which arise in this way are called toric ideals. See  for a thorough introduction to the subject.

It is natural to try to determine the syzygies of such an ideal $`I`$, a problem pursued in , or, more restrictedly, to ask for a minimal set of generators for such an ideal, the approach taken in . In our case, as in , we shall be interested in determining an upper bound for the degrees of a minimal generating set, in a special case, also singled out for consideration in , as follows: let $`d=2`$, and let all the $`a_{1j}=1`$. (It follows that the ideal $`I`$ is the homogeneous ideal of a monomial curve in projective space, but we shall not adopt that point of view here.) For simplicity, we refer to $`a_{2j}`$ as $`a_j`$. Without loss of generality let $`a_1<a_2<\mathrm{}<a_n`$.

The ring $`R`$ has a $`\times `$ grading, where $`z_j`$ has degree $`(1,a_j)`$. It is easily seen that the ideal $`I`$ is homogeneous with respect to this grading. Forgetting the $``$ component of the grading, we recover the usual $``$ grading on $`R`$. We use the word “bidegree” to refer to degree in the $`\times `$ grading, and “degree” to refer to degree in the usual $``$ grading.

This paper consists of a proof of the following theorem:

###### Theorem 1 (Main Theorem).

Let $`rs`$ be the two largest successive differences $`a_ia_{i1}`$. Then $`I`$ is generated by elements of degree no more than $`r+s`$.

Simple examples show this bound is tight: let $`r`$ and $`s`$ be relatively prime integers, let $`n=3`$, $`a_1=r`$, $`a_2=0`$, $`a_3=s`$. Then $`I`$ is the principal ideal generated by $`z_2^{r+s}z_1^sz_3^r`$.

Note that a given toric ideal $`I`$ will arise from more than one choice of $`𝒜`$. Given a toric ideal $`I`$ in $`R`$, to use the Main Theorem to obtain the best possible bound for the degrees of a minimal generating set for $`I`$, one should choose $`𝒜`$ in such a way that the greatest common divisor of the successive differences $`a_ia_{i1}`$ is 1. Such a choice is always possible, and any such choice will yield the same bound.

In , L’vovsky used cohomological results from  to prove a stronger statement than our Main Theorem, bounding the regularity of $`I`$, which is to say, bounding the degrees of the $`i`$-syzygies of $`I`$ for all $`i`$. Restricted to generators (0-syzygies), his result coincides with ours. It is not clear how the bound given by our Main Theorem compares with the bound obtained in .

In justification for this paper, aside from the intrinsic interest of a combinatorial proof of L’vovsky’s bound, we hope that the techniques of this proof may extend to higher dimensional cases.

The Main Theorem follows easily from the following combinatorial result.

###### Theorem 2 (Connectedness Theorem).

Let $`V`$, not necessarily finite, with the sizes of gaps between successive elements bounded above. Let $`rs`$ be the sizes of the two largest gaps between successive elements of $`V`$. For $`(q,c)\times `$, let $`\mathrm{\Pi }_{(q,c)}`$ be the collection of multisets with support in $`V`$, of cardinality $`q`$ and sum $`c`$. Let $`\mathrm{\Delta }_{(q,c)}`$ be the simplicial complex generated by the supports of the multisets in $`\mathrm{\Pi }_{(q,c)}`$. Then if $`q>r+s`$, $`\mathrm{\Delta }_{(q,c)}`$ is connected.

## Translation to Combinatorics

A multiset is an unordered collection of elements, in which some elements may appear with multiplicity greater than one. We use $`+`$ and $``$ for addition and subtraction of multisets, and write $`P=\{x_1,\mathrm{},x_n\}_{}`$ to indicate that the elements of $`P`$ are listed in non-decreasing order. $`P`$ is the sum of the elements of $`P`$. We now give the proof of the Main Theorem assuming the Connectedness Theorem:

###### Proof of Main Theorem.

Let $`V=\{a_1,\mathrm{},a_n\}_<`$. Let $`rs`$ be the two largest successive differences $`a_ia_{i1}`$. By the Connectedness Theorem, $`\mathrm{\Delta }_{(q,c)}`$ is connected for $`q>r+s`$. The Main Theorem now follows from the following lemma:

###### Lemma 1 (Translation Lemma).

No minimal generating set for $`I`$ has generators in bidegree $`(q,c)`$ iff $`\mathrm{\Delta }_{(q,c)}`$ is connected.

Remark: This is a special case of a result of , which gives information about the degrees of minimal generators of $`I`$ and also of all its $`i`$-syzygies, based on the homology of $`\mathrm{\Delta }_{(q,c)}`$. In the interest of self-containedness, we give an elementary proof of the result we need.

###### Proof.

($``$) Let $`I^<`$ denote the ideal of $`R`$ generated by the elements of $`I`$ of bidegree $`(q^{},c^{})`$ with $`q^{}<q`$. Then $`I_{(q,c)}^<`$, the $`(q,c)`$-bigraded part of $`I^<`$, is a sub-vector space of $`I_{(q,c)}`$. We wish to show that $`I_{(q,c)}=I_{(q,c)}^<`$.

For $`P`$ a multiset with support in $`V`$, let $`z^P`$ denote the monomial in $`R`$ where the exponent of $`z_i`$ is the multiplicity of $`a_i`$ in $`P`$. $`I_{(q,c)}`$ is spanned as a $`k`$-vector space by elements of the form $`z^Pz^P^{}`$, for $`P`$ and $`P^{}`$ in $`\mathrm{\Pi }_{(q,c)}`$.

Let $`P`$ and $`P^{}`$ be two elements of $`\mathrm{\Pi }_{(q,c)}`$ with non-empty intersection, say $`A`$. Then $`z^Pz^P^{}=z^A(z^{PA}z^{P^{}A})`$. But $`z^{PA}z^{P^{}A}`$ is in $`I`$, in bidegree $`(|PA|,(PA))`$. Thus $`z^Pz^P^{}I_{(q,c)}^<`$.

Now suppose $`P`$ and $`P^{}`$ are arbitrary elements of $`\mathrm{\Pi }_{(q,c)}`$. Since $`\mathrm{\Delta }_{(q,c)}`$ is connected, it follows that there exist $`P=P_0,P_1,\mathrm{},P_t=P^{}`$ with each $`P_i\mathrm{\Pi }_{(q,c)}`$, such that for each $`i`$, $`P_i`$ and $`P_{i+1}`$ have at least one element in common. By the previous argument, then, $`z^{P_i}z^{P_{i+1}}I_{(q,c)}^<`$, from which it follows that $`z^Pz^P^{}I_{(q,c)}^<`$, as desired.

($``$) $`I_{(q,c)}^<`$ is spanned by elements of the form $`z^A(z^Bz^B^{})`$ with $`A`$ non-empty, or equivalently of the form $`z^Pz^P^{}`$ with $`P`$ and $`P^{}`$ having non-empty intersection. It follows that for any $`z^Pz^P^{}I_{(q,c)}^<`$, the supports of $`P`$ and $`P^{}`$ are in the same component of $`\mathrm{\Delta }_{(q,c)}`$. Thus, if $`\mathrm{\Delta }_{(q,c)}`$ has more than one component, $`I_{(q,c)}^<I_{(q,c)}`$, as desired. ∎

## Combinatorial Lemmas

We now develop the combinatorial tools to prove the Connectedness Theorem. Fix $`V`$, with $`r`$ and $`s`$ the sizes of the two largest gaps between successive elements of $`V`$. Fix $`q>r+s`$ and fix $`c`$. Let $`\mathrm{\Delta }=\mathrm{\Delta }_{(q,c)}`$ and $`\mathrm{\Pi }=\mathrm{\Pi }_{(q,c)}`$.

Given a multiset $`P=\{x_1,\mathrm{},x_p\}_{}`$, define $`m(P)=_iix_i`$. Among multisets of the same sum and cardinality, the intuition is that $`m`$ is a measure of how spread out $`P`$ is — the more spread out $`P`$ is, the greater $`m(P)`$ will be.

We now show that, for a suitable class of $`P\mathrm{\Pi }`$, we can find another multiset $`P^{}\mathrm{\Pi }`$ which is more spread out than $`P`$.

###### Lemma 2 (Expansion Lemma).

Let $`P\mathrm{\Pi }`$, $`P=A+C`$. $`|C|=r+s`$, and $`C`$ contains neither the greatest element nor the least element of $`V`$. Then there exists some $`P^{}\mathrm{\Pi }`$, $`P^{}=A+C^{}`$, such that $`m(P^{})>m(P)`$.

Note: the following proof owes its basic approach to the proof of Theorem 6.1 of .

###### Proof.

Consider the following algorithm, which obtains a sequence of multisets $`C_i`$, with $`C_0=C`$, where $`C_i`$ is obtained from $`C_{i1}`$ by replacing one of the original elements of $`C`$ by either the next larger or the next smaller element of $`V`$. When thinking about this algorithm, it’s helpful to think of the elements of $`C`$ as stones sitting on a number line, where the allowed positions for the stones are the numbers in $`V`$. $`C_i`$ is obtained from $`C_{i1}`$ by jumping one stone which hasn’t been moved yet to the next higher or lower allowed position.

###### Algorithm 1 (Expansion Algorithm).

$`C_0:=C`$.
$`Active:=C`$. (These are element of $`C`$ which haven’t moved yet.)
$`s_0:=0`$.
For $`i:=1`$ to $`r+s`$ do:
If $`s_{i1}0`$,
Remove the largest element from Active, and call it $`x_i`$.
$`C_i:=C_{i1}\{x_i\}+\{\text{the least element of }V\text{ greater than }x_i\}`$.
If $`s_{i1}>0`$,
Remove the smallest element from Active, and call it $`x_i`$.
$`C_i:=C_{i1}\{x_i\}+\{\text{the greatest element of }V\text{ less than }x_i\}`$.
Let $`s_i:=C_iC`$.

We would now like to bound $`s_i`$. Observe that if $`s_{i1}0`$ then $`s_i>s_{i1}`$, and if $`s_{i1}>0`$, then $`s_i<s_{i1}`$. Thus, the absolute value of $`s_i`$ can be no greater than the largest jump possible on a single step, which is $`r`$. However, we can be a little more precise.

The elements of $`C`$ which are increased by the algorithm are greater than or equal to all the elements which are decreased by the algorithm. Thus, any gap between successive elements of $`V`$ can be jumped in only one direction in the course of running this algorithm (though it may be jumped more than one time).

Suppose there is a unique gap of size $`r`$, and it is jumped only in the increasing direction. Then it follows that $`s+1s_ir`$ for all $`i`$. Symmetrically, suppose there is a unique gap of size $`r`$, and it is jumped only in the decreasing direction. Then $`r+1s_is`$. If there is a unique gap of size $`r`$ which is not jumped, or there is more than one gap of size $`r`$ (in which case $`r=s`$), then $`s+1s_is`$. In any case, we deduce that there are at most $`r+s`$ possible values for $`s_i`$. But there are $`r+s+1`$ of the $`s_i`$. Thus, at least two of the $`s_i`$ must be equal, say $`s_j=s_l`$, with $`j<l`$. Obtain $`C^{}`$ from $`C`$ by making the same jumps as were made in the algorithm on steps $`j+1`$ through $`l`$. Then $`C^{}=C+s_ls_j=C`$. Let $`P^{}=A+C^{}`$. In going from $`P`$ to $`P^{}`$, the elements which have been increased are all greater than or equal to the elements which have been decreased, and thus $`m(P^{})>m(P)`$. ∎

Using only this lemma, we can prove the Connectedness Theorem with the additional assumption that $`V`$ is not bounded below.

###### Proof of Connectedness Theorem assuming $`V`$ is not bounded below.

Let $`x`$ and $`y`$ be vertices of $`\mathrm{\Delta }`$. We want to show that they lie in the same component. Choose some $`P\mathrm{\Pi }`$ with $`xP`$, and some $`P^{}\mathrm{\Pi }`$, with $`yP^{}`$. Let $`x^{}`$ be the maximum element of $`P`$, $`y^{}`$ the maximum element of $`P^{}`$. If $`x^{}=y^{}`$ then we are done. So assume without loss of generality that $`x^{}<y^{}`$. Choose a subset $`A`$ of $`P`$, whose size is $`q(r+s)`$, and which contains $`x`$. Thus, we can write $`P=A+C`$, satisfying the hypotheses of the Expansion Lemma.

Now apply the Expansion Lemma recursively, with this choice of $`A`$ fixed, obtaining a sequence of multisets $`P=P_0,P_1,P_2,\mathrm{},P_t`$ until either $`m(P_t)(1+2+\mathrm{}+q)y^{}`$ or $`P_t`$ contains the largest element of $`V`$. In either case, $`P_t`$ clearly contains an element which is greater than or equal to $`y^{}`$. Thus, some previous $`P_i`$ contains $`y^{}`$. Now by construction $`P_i\mathrm{\Pi }`$, and $`x`$ and $`y^{}`$ are both contained in $`P_i`$, which finishes the proof. ∎

###### Lemma 3 (Multiple Expansion Lemma).

Let $`P,P^{}\mathrm{\Pi }`$. Suppose the largest and smallest elements of $`P+P^{}`$ occur in $`P^{}`$. Then there is a $`Q\mathrm{\Pi }`$ containing at least one element of each of $`P`$ and $`P^{}`$.

###### Proof.

The proof is essentially the argument given above, proving the Connectedness Theorem in the case where $`V`$ is not bounded below (but we do not make the assumption that $`V`$ is not bounded below). Fix some set $`A`$ in $`P`$ of size $`q(r+s)`$. Apply the Expansion Lemma recursively with $`A`$ fixed, producing a sequence of multisets $`P=P_0,P_1,\mathrm{}`$, until some $`P_t`$ includes some element of $`P^{}`$. The Expansion Lemma can always be applied because at no stage before halting does $`P_i`$ include the largest or smallest element of $`V`$, since prior to including one of these elements, it would include an element of $`P^{}`$. Then $`Q=P_t`$ satisfies the conditions in the statement of the lemma. ∎

We now prove another lemma, similar to the Multiple Expansion Lemma, which we will need to prove the Connectedness Theorem in full generality.

###### Lemma 4 (Criss-Cross Lemma).

Let $`P,P^{}\mathrm{\Pi }`$. Suppose the largest element of $`P+P^{}`$ occurs in $`P^{}`$, while the smallest element of $`P+P^{}`$ occurs in $`P`$. Then there exists a $`Q\mathrm{\Pi }`$ which contains at least one element of each of $`P`$ and $`P^{}`$.

###### Proof.

Assume, without loss of generality, that $`P`$ and $`P^{}`$ are disjoint. Pick $`fP`$, $`f^{}P^{}`$, with $`f<f^{}`$. Let $`B=P\{f\}`$, $`B^{}=P^{}\{f^{}\}`$.

Split $`B`$ into two multisets, $`X`$ and $`Y`$, where $`X`$ consists of the elements of $`B`$ less than all elements of $`B^{}`$, and $`Y`$ is the remainder. Similarly, split $`B^{}`$ into $`X^{}`$, the elements greater than all elements of $`B`$, and $`Y^{}`$, the remainder.

###### Lemma 5 (Size Lemma).

Either $`|Y|`$ is greater than the longest gap below $`B^{}`$ or $`|Y^{}|`$ is greater than the longest gap above $`B`$.

###### Proof.

First, I claim that $`|Y|+|Y^{}|>r+s`$. Let $`B=\{b_1,\mathrm{},b_{q1}\}_{}`$, $`B^{}=\{b_1^{},\mathrm{},b_{q1}^{}\}_{}`$. For $`1i|X|`$, and for $`q|X^{}|iq1`$, $`b_i<b_i^{}`$. So, if $`|X|+|X^{}|q1`$, then $`B<B^{}`$. But $`B=Pf>P^{}f^{}=B^{}`$, which is a contradiction. Thus, $`|X|+|X^{}|<q1`$, so $`|Y|+|Y^{}|>q1r+s`$.

Now, we can see that at least one of the following four cases holds:

1. There is a unique gap of size $`r`$ below $`B^{}`$ and $`|Y|r+1`$.
2. All the gaps below $`B^{}`$ are of size no more than $`s`$ and $`|Y|s+1`$.
3. There is a unique gap of size $`r`$ above $`B`$ and $`|Y^{}|r+1`$.
4. All the gaps above $`B`$ are of size no more than $`s`$ and $`|Y^{}|s+1`$.

If there is a unique gap of size $`r`$ below $`B^{}`$ then either $`|Y|r+1`$ (case (1)) or $`|Y^{}|s+1`$ (case (4)). Similarly, if there is a unique gap of size $`r`$ above $`B`$, we are in case (2) or case (3). Otherwise, we are in case (2) or (4).

It is now easy to see that in cases (1) or (2), $`|Y|`$ is greater than the largest gap below $`B^{}`$, while in cases (3) or (4), $`|Y^{}|`$ is greater than the largest gap above $`B`$, proving the lemma. ∎

Now, using the Size Lemma, by symmetry, we may assume without loss of generality that $`|Y|`$ is greater than the largest gap below $`B^{}`$. Let the size of this gap be $`g`$. Now consider the following algorithm.

###### Algorithm 2 (Criss-Cross Algorithm).

$`B_0:=B`$.
$`ActiveX:=X`$. (These are the elements of $`X`$ that haven’t moved yet.)
$`ActiveY:=Y`$. (And similarly for $`Y`$.)
$`i:=0`$.
$`s_0:=0`$.
Repeat until $`s_i=s_j`$ for some $`j<i`$:
$`i:=i+1`$.
If $`s_{i1}0`$, then
If ActiveX is non-empty,
Remove the largest element of ActiveX and call it $`x_i`$.
Otherwise,
Remove the largest element of ActiveY and call it $`x_i`$.
$`B_i:=B_{i1}\{x_i\}+\{\text{the least element of }V\text{ greater than }x_i\}`$.
If $`s_{i1}>0`$,
Remove the smallest element of ActiveY and call it $`x_i`$.
$`B_i:=B_{i1}\{x_i\}+\{\text{the greatest element of }V\text{ less than }x_i\}`$.
Let $`s_i=B_iB`$.

I first argue that there are always enough elements of ActiveX and ActiveY to run this algorithm. There are two potential problems:

First, suppose $`s_{i1}0`$, $`\text{ActiveX}=\mathrm{}`$, $`\text{ActiveY}=\mathrm{}`$. Then $`iq`$. So at least $`r+s+1`$ values of $`s_j`$ have been defined. But, as in the proof of the Expansion Lemma, no gap can get jumped in both directions, so there are at most $`r+s`$ possible values for the $`s_j`$, so the algorithm could never have reached this step, since there must already have been some value $`s_j`$ which has occured more than once.

Second, suppose $`s_{i1}>0`$, and $`\text{ActiveY}=\mathrm{}`$. If $`\text{ActiveX}=\mathrm{}`$, then $`iq`$, and we have a contradiction as in the previous case. So assume that $`\text{ActiveX}\mathrm{}`$.

Thus, every element of $`Y`$ that have been removed from ActiveY was removed on a turn $`j`$ with $`s_j>0`$. Further, we can say that $`gs_j`$ for all these $`j`$, since we have not yet jumped any elements of $`Y`$ in the positive direction, so all the gaps that we have jumped in the positive direction are no longer than $`g`$. Since the algorithm has not yet terminated, all the $`s_j`$ are different, and they have only $`g`$ different possible positive values. Thus, the case $`s_j>0`$ has been encountered at most $`g`$ times. But since $`|Y|g+1`$, it follows that $`\text{ActiveY}\mathrm{}`$, which contradicts the assumption that we had run out of ActiveY.

Thus, the algorithm runs successfully, and we obtain $`i>j`$ with $`s_i=s_j`$. Obtain $`\stackrel{~}{B}`$ from $`B`$, by making the same jumps as were made on steps $`j+1`$ through $`i`$ of the algorithm. As in the proof of the Expansion Lemma, $`\stackrel{~}{B}=B`$.

Let $`P_0:=P`$, $`P_1:=\{f\}+\stackrel{~}{B}`$. If $`P_1`$ contains an element of $`B^{}`$, we are done with $`Q=P_1`$. Otherwise, we apply the algorithm recursively, keeping $`B^{}`$, $`f`$, and $`f^{}`$ fixed. I claim that we eventually obtain a $`P_i`$ which contains an element of $`B^{}`$.

Note first that, until we obtain a $`P_i`$ which contains an element of $`B^{}`$, we stay in the same case of the Size Lemma, because no element of $`B`$ can jump past an element of $`B^{}`$ without hitting it, so $`X`$ and $`Y`$ will stay the same size. Thus, we can succesfully repeat the algorithm until some $`P_i`$ intersects $`B^{}`$.

For a multiset $`D=\{x_1\mathrm{},x_p\}_{}`$, define $`\stackrel{~}{m}(D)=_i(p+1i)x_i`$. This function plays a similar role to $`m`$, except that as a set gets more spread out, $`\stackrel{~}{m}`$ decreases.

Now, I claim that $`\stackrel{~}{m}(Y)`$ is decreased each time we run the algorithm, since the decreases to $`Y`$ are at least as great as its increases, and the elements that are decreased precede the elements that are increased.

Thus, at some point, if the algorithm could run forever, $`\stackrel{~}{m}(Y)`$ would be smaller than $`(1+2+\mathrm{}+|Y|)`$ times the smallest element of $`B^{}`$, but this would mean that $`Y`$ contained elements smaller than every element of $`B^{}`$, which would be a contradiction.

So eventually, some $`P_i`$ contains an element of $`B^{}`$. Then $`Q=P_i`$ satisfies the conditions of the statement of the lemma. ∎

## Proof of the Connectedness Theorem

Now, we prove the connectedness of $`\mathrm{\Delta }`$ in full generality.

###### Proof of Connectedness Theorem.

We begin as in the case of $`V`$ not bounded below. Let $`x`$ and $`y`$ be vertices of $`\mathrm{\Delta }`$. We want to show that $`x`$ and $`y`$ are in the same component. Choose a set $`P`$ in $`\mathrm{\Pi }`$ which contains $`x`$, and a set $`P^{}`$ in $`\mathrm{\Pi }`$ which contains $`y`$.

If $`P`$ and $`P^{}`$ intersect, then we are done. So assume they do not. Suppose that the greatest and the least elements of $`P+P^{}`$ both occur in one of $`P`$ or $`P^{}`$, without loss of generality, say $`P^{}`$. Now apply the Multiple Expansion Lemma to obtain a multiset $`Q\mathrm{\Pi }`$ intersecting both $`P`$ and $`P^{}`$. As in the proof of the Connectedness Theorem with $`V`$ not bounded below, this shows that $`x`$ and $`y`$ are connected in $`\mathrm{\Delta }`$.

Now, suppose that the greatest and least elements of $`P+P^{}`$ do not both occur in either $`P`$ or $`P^{}`$. Without loss of generality, let $`P^{}`$ contain the greatest. This puts us in the position to apply the Criss-Cross Lemma, and again, the multiset $`Q`$ which we obtain from it shows that $`x`$ and $`y`$ are connected in $`\mathrm{\Delta }`$. ∎

## Further Directions

First, it would be good to give a combinatorial proof of the entire result of L’vovsky (bounding the degrees of all $`i`$-syzygies of $`I`$, not just the generators). An argument might use the full strength of the result mentioned in the proof of the Translation Lemma from  to translate the problem into a combinatorial framework. The necessary combinatorial result might then be established by an induction argument with a version of the Connectedness Theorem as a base case, but so far we have been unable to accomplish this except in the case where $`V`$ is not bounded below, which is of limited interest.

Also, as mentioned in the introduction, one might apply the techniques of this paper to prove degree bounds or regularity bounds in higher dimensions. However, this is considerably trickier. Our strategy for modifying a multiset $`P`$ until it hits an element of $`P^{}`$ is akin to a game of hide-and-seek in one dimension; the interested reader will see why the game is usually played in two dimensions.

## Acknowledgments

I would like to thank Hal Schenck for suggesting this problem to me, and for fruitful initial discussions. I would also like to thank Philippe Giménez and the referee for their helpful comments.


text:
# Condensate of a charged boson fluid at non-integer dimensions

## I introduction

Charged boson fluid(CBF) has been attracting attention as an interesting model of many physical systems at low temperatures. Several theoretical and computational studies have been carried out on the 2D and the 3D CBF at zero and finite temperatures. Also, a general expression for the superfluid density of dilute Bose gas for $`D>2`$ was given by Fisher and Hohenberg

It has been, recently, noted that porous media and films can be studied using fractal dimensionality between $`2<D<3`$. This study also revealed that several prominent features of porous media originates from dimensionality and detailed nature of mutual interaction plays a much less prominent role.

In this paper, we study the condensate and the superfluid density of CBF at finite temperatures in the dimensionality between $`2<D<3`$. Then, the result will be compared with the experimental data from superconducting thin films. It will be shown that the nature of gap states and detailed property of superconducting mechanism plays a much less prominent role than that of the geometrical factors.

We will consider a model of point-like spinless charged bosons embedded in a uniform neutralizing background and coupled by a Coulomb type interaction. The Bogoliubov approximation will be applied for the model Hamiltonian in weakly interacting regime to calculate the condensate. The condensate fraction will be calculated as a function of temperature and dimensionality between 2 and 3. Then, the result is utilized to calculate the superfluid fraction with the help of the $`D`$-dimensional Fisher and Hohenberg formula.

Finally, the theoretical results on the charged bosons in fractal dimensions between 2 and 3 will be compared to experimental data on superfluid density of superconducting films and the physical implications will be discussed.

## II Bogoliubov approach to a charged boson fluid

We begin our theoretical scheme from the well-established Bogoliubov approach for the CBF on a uniform neutralizing background and generalize to non-integer dimensions. The model Hamiltonian is given by

$$H=𝑑𝐫\psi ^{}(𝐫)\left(\frac{^2}{2m}\mu \right)\psi (𝐫)+\frac{1}{2}𝑑𝐫𝑑𝐫^{}\psi ^{}(𝐫)\psi ^{}(𝐫^{})V(|𝐫𝐫^{}|)\psi (𝐫^{})\psi (𝐫),$$
(1)
where $`V(|𝐫𝐫^{}|)=Q^2/|𝐫𝐫^{}|`$, $`m`$ is the effective mass, and $`Q`$ is the effective charge of the CBF. $`\psi (𝐫)`$ is the boson field operator and $`\mu `$ is the chemical potential. $`\mathrm{}=c=1`$ is assumed for convenience.

D-dimensional coupling strength is measured by the dimensionless parameter $`r_s`$ defined by $`Ar_s^D=(na_B^D)^1`$. A large $`r_s`$ means a long range of interaction and a small density of particles. $`A`$ is the D-dimensional volume parameter given by $`A=2\pi ^{D/2}D^1\mathrm{\Gamma }(D/2)^1`$ and $`a_B`$ is the effective Bohr radius of the system defined by $`a_B=1/mQ^2`$.

We first apply the Bogoliubov transformation at T=0. The field operator is written as

$$\psi (𝐫,t)=\psi _0+\stackrel{~}{\psi }(𝐫,t),$$
(2)
where $`\psi _0`$ is a macroscopic order parameter representing the condensate. $`\stackrel{~}{\psi }`$ describes the particles promoted out of the condensate and can be expressed as a linear transformation

$$\stackrel{~}{\psi }(𝐫,t)=\underset{k}{}\left[u_k(𝐫,t)a_k+v_k^{}(𝐫,t)a_k^{}\right],$$
(3)
where $`a_k(a_k^{})`$ is the bosonic quasi-particle operator. Also, we have $`\mu =0`$, and $`n_0=\psi _0^2`$ for the uniform condensate.

A straightforward calculation produces the following results in D-dimensions

$`u_k^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[ϵ_k^1\left\{n_0U_D(k)+{\displaystyle \frac{k^2}{2m}}\right\}+1\right],`$ (4)
$`v_k^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[ϵ_k^1\left\{n_0U_D(k)+{\displaystyle \frac{k^2}{2m}}\right\}1\right],`$ (5)
where

$$ϵ_k=\left\{\frac{n_0k^2U_D(k)}{m}+\left(\frac{k^2}{2m}\right)^2\right\}^{1/2}.$$
(6)
Here, $`ϵ_k`$ is the energy spectrum of the Bogoliubov transformation. Note that $`u_k^2v_k^2=1`$. $`U_D(k)`$ is the Fourier transformation of the $`D`$-dimensional coupling $`V(r)`$ and given by $`U_D(k)=DAQ^2/k^{D1}`$. Note that $`U_3(k)=4\pi Q^2/k^2`$, and $`U_2(k)=2\pi Q^2/k`$ as expected.

## III condensate fraction at finite temperatures and non-integer dimensions

The extension of the theory to finite temperatures is effected by the Bose distribution function, $`f_k=1/\left\{\mathrm{exp}(\beta ϵ_k)1\right\}`$. The condensate density at nonzero temperature in D-dimensions is given by

$`n_0`$ $`=`$ $`n<\psi ^{}(𝐫)\psi ^{}(𝐫)>`$ (7)
$`=`$ $`n{\displaystyle \underset{𝐤0}{}}\left[v_k^2+f_k(u_k^2+v_k^2)\right].`$ (8)
Substituting Eqs. (4), (5), and (6) into Eq. (8), we obtain the condensate fraction in D-dimensions as follows

$`1{\displaystyle \frac{n_0}{n}}`$ $`=`$ $`{\displaystyle \frac{1}{n}}{\displaystyle \underset{𝐤0}{}}\left[v_k^2+f_k(u_k^2+v_k^2)\right]`$ (9)
$`=`$ $`{\displaystyle \frac{1}{2n}}{\displaystyle _0^{\mathrm{}}}𝑑k{\displaystyle \frac{DAk^{D1}}{(2\pi )^D}}\left[{\displaystyle \frac{n_0U_D(k)+\frac{k^2}{2m}}{\sqrt{\frac{n_0k^2U_D(k)}{m}+\left(\frac{k^2}{2m}\right)^2}}}\left(1+{\displaystyle \frac{2}{e^{\beta ϵ_k}1}}\right)1\right].`$ (10)
Eq. (10) can be converted into a compact form after some mathematical steps,

$$\left(1\frac{n_0}{n}\right)\left(\frac{n_0}{n}\right)^{D/(D+1)}=F(D)I(T,D).$$
(11)
$`F(D)`$ is the coefficient which is defined by

$$F(D)=\frac{3}{D+1}\frac{2^{(D^2+3D+2)/(D+1)}}{D^{(D+2)/(D+1)}}\frac{r_s^{D/(D+1)}}{\mathrm{\Gamma }^2\left(D/2\right)}.$$
(12)
$`I(D,T)`$ is the integral which is defined by

$$I(D,T)=_0^{\mathrm{}}𝑑xx^{(D2)(D3)/D(D+1)}\left[\frac{p(x)}{q(x)}\left\{1+\frac{2}{e^{C(D,T)g(x)q(x)}1}\right\}2x^{6/D}\right],$$
(13)
where

$`p(x)`$ $`=`$ $`1+2x^{12/D},`$ (14)
$`q(x)`$ $`=`$ $`\sqrt{1+x^{12/D}},`$ (15)
$$C(D,T)=\left(\frac{Dn_0}{n}\right)^{2/(D+1)}\frac{2^{(3D)/(D+1)}}{r_s^{2D/(D+1)}}\frac{Q^2}{a_Bk_BT},$$
(16)
and

$$g(x)=x^{6(3D)/D(D+1)}.$$
(17)
The temperature dependence is enclosed in the function $`C(D,T)`$. The integral $`I(D,T)`$ is performed numerically.

The condensate fraction in D-dimensions is plotted as a function of temperature in FIG. 1. The unit of the temperature is $`Q^2/a_Bk_B`$. The three coupling strengths were chosen for the calculation: (a) $`r_s=0.8`$ for a weak coupling and (b) $`r_s=1.2`$ for an intermediate coupling, and (c) $`r_s=2.0`$ for a strong coupling, respectively. We see that the condensate fraction depends strongly on the interaction strength, $`r_s`$. We also observe that for all the different values of the coupling strength, the condensate fractions have significant nonzero values and are relatively flat near $`T0K`$ region.

## IV superfluid fraction at finite temperatures and non-integer dimensions

In order to compare the theory with available experimental data on superconducting films, it is necessary to obtain the superfluid fraction which is directly related to the penetration depth of superconducting films. The superfluid density of dilute Bose gas in $`D`$-dimensions has been studied by Fisher and Hohenberg . It is given by a Landau quasi-particle formula based on the $`D`$-dimensional Bogoliubov energy spectrum $`ϵ_k`$ in Eq. (6) as

$$\frac{n_s}{n}=1\frac{\beta }{2nm}\underset{𝐤0}{}k^2\frac{e^{\beta ϵ_k}}{(e^{ϵ_k}1)^2}.$$
(18)
Substitution of Eq. (6) into the Eq. (18) gives an explicit expression for the D-dimensional superfluid fraction,

$$\frac{n_s}{n}=1H(D,T)_0^{\mathrm{}}𝑑xx^{D+1}\frac{e^{B(D,T)\sqrt{x^{3D}+x^4}}}{\{e^{B(D,T)\sqrt{x^{3D}+x^4}}1\}^2}.$$
(19)
$`B(D,T)`$ and $`H(D,T)`$ are defined by the following expressions:

$$B(D,T)=2^{(3D)/(D+1)}\left(\frac{Dn_0r_s}{n}\right)^{2/(D+1)}\frac{Q^2}{a_Bk_BT},$$
(20)
and

$$H(D,T)=\frac{A^2D^{(2D+1)/(D+1)}}{(2\pi )^D}\left(\frac{4n_0r_s}{n}\right)^{D/(D+1)}B(D,T).$$
(21)
The condensate fraction $`n_0/n`$ of Eq. (11) is substituted into Eqs. (20) and (21) for the calculation of the superfluid fraction $`n_s/n`$.

The $`D`$-dimensional superfluid fraction is obtained from Eq. (19), and plotted as a function of temperature for three different values of the coupling strength in FIG. 2. The same $`r_s`$ values as in FIG. 1 were chosen for the comparison.

We find the basic structure between the condensate fraction in FIG. 1 and the superfluid fraction in FIG. 2 is similar for any coupling strength, $`r_s`$. Therefore, it is clear that condensate density gives a strong and useful hint for the superfluid density for the CBF system.

## V superfluid density of superconducting films

Measurements on the temperature dependent superfluid density of superconducting films $`La_{2x}Sr_xCuO_4`$ and $`Mo_{77}Ge_{23}`$ drew considerable interests, in connection with the puzzling splitting behaviors shown in FIG. 3. In FIG. 3 we note that $`1/\lambda ^2`$ is proportional to the superfluid density , and inverse of the sheet inductance, $`1/L(T)`$, is also proportional to the areal superfluid density, $`n_sd`$.

We observe that the two figures in FIG. 3 exhibit the same generic splitting behavior for samples with various thickness. Considering differences in the detailed physical properties between LSCO and MoGe films, such qualitative similarity is quite surprising and, thus, it strongly implies that the generic behavior stems from some hidden common properties.

There have been several theoretical efforts to explain the general features of the superfluid density of superconducting films. It has been shown that the simple BCS and the d-wave symmetry just lead into a quadratic dependence of $`1/\lambda ^2(T)`$ and, thus, can not explain the generic splitting behavior. Effect of thermal fluctuations is shown to explain the relative flat behavior of curves at low temperatures. But it could not also explain the above mentioned behavior.

Since the superconducting mechanisms of LSCO is believed to be basically different from that of MoGe, any theory which explains the behavior should not be based on the microscopic details of the superconducting mechanisms of the two materials. Instead, it should mainly reflect the common geometric nature of superconducting films.

It has been shown that condensation of bound pairs can be used to explain the superconductivity with proper scaling. In the previous section, we have shown that the superfluid density, which originates from the charged boson condensation in fractal dimensions, $`2<D<3`$, exactly duplicates the observed experimental data from the superconducting films. Therefore, we believe that the theoretical results in FIG. 2 are the natural explanation of the experimental observations in FIG. 3.

We note also that the general features of the superfluid density curves of the superconducting films are qualitatively similar to those of liquid helium-4 in films and porous media.

## VI Conclusions

In this paper, we have studied the condensate of charged boson fluid at finite temperatures in non-integer dimensions between 2 and 3. The condensate and superfluid fraction are obtained as functions of temperature and dimensions at various values of the coupling strength.

We have shown that the generic splitting behavior of the superfluid density universal for superconducting films may originate from the geometric nature of the films and detailed nature of the superconducting mechanism plays a much less prominent role than generally believed.

###### Acknowledgements.

This work was supported by the Korea Science and Engineering Foundation through a special program and the Center for Strongly Correlated Material Research(SNU).


text:
# Mesoscale simulations of surfactant dissolution and mesophase formation

(21 April 2002 — PREPRINT)

## Abstract

The evolution of the contact zone between pure surfactant and solvent has been studied by mesoscale simulation. It is found that mesophase formation becomes diffusion controlled and follows the equilibrium phase diagram adiabatically almost as soon as individual mesophases can be identified, corresponding to times in real systems of order $`10\mu \mathrm{s}`$.

When pure surfactant comes into contact with water, mesophases appear at the interface. This is an important process not only for the practical use of surfactants, but also from the point of view of fundamental surfactant phase science. Indeed, contact ‘flooding’ or penetration scan experiments can yield quantitative information on mesophases as a function of composition Laughlin (1992). The phenomenon is invariably *diffusion controlled* in the sense that the widths of the mesophases follow $`t^{1/2}`$ growth laws, and *adiabatic* in the sense that the local composition determines the mesophase boundaries according to the equilibrium phase diagram Warren and Buchanan (2001); Chen et al. (2000); Laughlin and Munyon (1987). But penetration scan experiments are restricted to observation times of minutes to hours: what happens on time scales shorter than this? How early does diffusion control set in, and how soon can one expect the mesophase boundaries to track the local composition adiabatically? Such questions are not just of scientific interest since time scales of seconds or less are important in modern processing and the everyday use of surfactants, for instance determining how rapidly washing powder dissolves. Experimentally, this regime is very difficult to access because of the short time scales and the relatively small amounts of mesophase involved. To probe these questions therefore, we have therefore undertaken novel mesoscale simulations of surfactant dissolution. We find adiabatic diffusion control is established remarkably rapidly in our model, on time scales in which only a few repeat units of the growing mesophases have appeared, corresponding to times in the real systems of order $`10\mu \mathrm{s}`$.

The model we have used is a ‘minimalist’ particle-based model of a binary surfactant / water mixture, based on the dissipative particle dynamics (DPD) method Groot and Warren (1997); Jury et al. (1999). In DPD, the particles are soft spheres, interacting with pairwise soft potentials of the form $`U=\frac{1}{2}A(1r/r_c)^2`$ ($`r<r_c`$) where $`r`$ is the particle separation, $`r_c`$ is the range of the interaction, and $`A`$ the amplitude. In the model we have three species of particles: A, B and C. The A and B particles are bound together in pairs as dimers with a fixed separation $`r_d`$, and represent the surfactants. The C particles are monomers representing the solvent (water). The different species are distinguished by their interaction amplitudes, and the trick is to find a set of amplitudes which recover suitable phase behaviour. Following earlier work Jury et al. (1999), we use $`A_{\mathrm{AA}}=A_{\mathrm{BB}}=A_{\mathrm{CC}}=25`$, $`A_{\mathrm{AB}}=30`$, $`A_{\mathrm{AC}}=0`$, $`A_{\mathrm{BC}}=50`$ and $`r_d=0.5`$ which gives phase diagram features lying in a convenient temperature range around $`k_\mathrm{B}T1`$ (we fix units by choosing $`m=r_c=1`$ where $`m`$ is the mass of the particles).

Fig. 1 shows the phase diagram as a function of dimer concentration and $`k_\mathrm{B}T`$, at an overall density $`\rho =(2N_{\mathrm{AB}}+N_\mathrm{C})/V=6`$, where dimer concentration is defined to be the mole fraction of particles in dimers: $`c=2N_{\mathrm{AB}}/(2N_{\mathrm{AB}}+N_\mathrm{C})`$. Similar to many real systems Laughlin (1994), micellar (L1), hexagonal (H1) and lamellar (Lα) phases are found in order of increasing dimer concentration, and the interactions are chosen so that on the pure dimer side there is an isotropic fluid (L2) phase. The fact that a certain amount of solvent is needed to induce the Lα phase reflects the real behaviour of many nonionic surfactants, and is crucial to the contact simulations below. We did not find any cubic phases, although these are not always present in real systems (an extensive study of lattice models by Larson sur  suggests that cubic phases could be engineered to appear particularly if the DPD model is elaborated beyond dimers). Also, we caution that $`k_\mathrm{B}T`$ in the model is not easily mapped onto a real temperature. Despite these deficiencies, it is nevertheless remarkable that such a simple model reproduces the main features of the phase diagram common to a large number of surfactants.

To simulate surfactant dissolution and mesophase formation with this model is now very simple. We take two simulation boxes, one containing an equilibrated fluid of pure dimers and the other containing equilibrated solvent particles. We place them next to each other and allow the dimers and solvent particles to interdiffuse. Just as in the real systems, mesophases start to appear at the interface over time. Here we report results for the $`k_\mathrm{B}T=1`$ isotherm, although the $`k_\mathrm{B}T=0.5`$ isotherm was also studied with similar results. We consider simulation boxes of size $`10^2\times 100`$, with the concentration gradient along the long axis. Fig. 1 shows some representative simulation snapshots.

Whilst this is straightforward, there are a couple of technical points to be considered. Firstly it is important to adjust the densities in the two boxes so that the pressures are equal. From separate simulations, we determined the equation of state for the pure components, and found that densities $`\rho =6.124`$ and 5.896 for the solvent and dimers respectively give a common pressure $`p100`$. These densities are within $`2`$% of the density $`\rho =6`$ used to construct the phase diagram in Fig. 1, which obviates the need to consider constant pressure simulations. The pressure matching is important though because it suppresses sound waves which would spoil the subsequent analysis. Secondly, if we were to use conventional periodic boundary conditions, there would be an unwanted second interface in the system. We eliminate this by bounding the simulation box in the long direction by hard reflecting walls, supplemented by short-range soft repulsive potentials of the form $`U=\frac{1}{2}A_{\mathrm{wall}}(1z/r_c)^2`$ ($`z<r_c`$), where $`z`$ is the distance from the hard wall. The soft repulsive force suppresses density oscillations which would otherwise arise due to the abrupt termination of the particle density. We find empirically $`A_{\mathrm{wall}}=25`$ gives a smoothly vanishing density. Periodic boundary conditions are retained in the other two directions.

Now we turn to an objective analysis of the simulations to determine the growth laws. This is done by introducing a local order parameter which can distinguish between mesophases as a function of distance normal to the original contact plane. Although we examined various possibilities such as the use of Minkowski functionals Michielsen and de Raedt (2001), an approach which works well in practice is motivated by the $`\rho _\mathrm{C}=1.3`$ isosurfaces shown in Fig. 1. The different mesophases are clearly distinguishable to the eye because the isosurface is predominantly spherical, cylindrical, planar, or fragmented, for the L1, H1, Lα and L2 phases respectively. This geometric insight can be made quantitative by constructing the second moment of the isosurface normal distribution $`p(𝐧)`$, namely the symmetric tensor $`𝐌=\mathrm{𝐧𝐧}p(𝐧)𝑑𝐧`$ VTK . The local geometric nature of the isosurface and hence the underlying mesophase is reflected in the eigenvalues $`\mu _i`$ of $`𝐌`$ as laid out in Table 1. If these eigenvalues are ranked in order of increasing size, we determined after some trials that a slice can be classified as H1 if $`\mu _1<0.05`$ and $`\mu _{2,3}>0.3`$, or Lα if $`\mu _1<0.05`$ and $`\mu _2<0.15`$ (note $`_{i=1}^3\mu _i=1`$).

For the present study, we divide the simulation box into slices of thickness $`0.25r_c`$ parallel to the original contact plane and determine the eigenvalues of $`𝐌`$ for each slice. Fig. 2(a) gives a representative example, showing that the various mesophases can be clearly distinguished. The transition between H1 and Lα is particularly sharp as can also be seen in Fig. 1(h). The above criteria are used to classify each slice and the positions of the boundaries between mesophases determined as a function of time. There is a short incubation period before the local order parameter can distinguish the different mesophases, but beyond this point, the boundaries can be reliably tracked as shown in Fig. 2(b). The key results are contained in Figs. 2(c) and (d), which show respectively the local composition at the mesophase boundary and the squared width of the middle two mesophases as a function of time.

Fig. 2(c) shows that the local compositions become established at constant values almost as soon as mesophases can be distinguished. Table 2 shows that the L1-H1 and H1-Lα boundaries in the kinetic simulation lie almost exactly at the point expected from the equilibrium phase diagram. The Lα-L2 boundary is displaced from that found in the equilibrium phase diagram, but it is clear from Fig. 2(a) that this boundary is not completely sharp. Apart from this boundary therefore, it appears that the mesophases grow adiabatically almost from the earliest moments that the local order parameter can distinguish the growing mesophases.

Fig. 2(d) shows that both mesophases follow a $`\mathrm{\Delta }z^2t`$ law, where $`\mathrm{\Delta }z`$ is the mesophase width. This is a classic verification of diffusion control. To further examine this, we fit the slope to what would be expected if the local composition followed a simple diffusion law:

$$\mathrm{\Delta }z^2=4D[\mathrm{erf}^1(2c_21)\mathrm{erf}^1(2c_11)]^2t,$$
(1)
where $`D`$ is an *effective* diffusion coefficient, and $`c_{1,2}`$ are the fixed compositions at the edges of the mesophase of interest. The results are shown in Table 3, where for comparison we have included the diffusion coefficients for the two pure components, ie dimers in a pure dimer fluid and solvent particles in pure solvent, determined from separate simulations. The effective diffusion coefficients are slightly different in the two phases, similar to the finding our previous experiments Chen et al. (2000), although the difference here is within the statistical errors. Thus we conclude that, whilst the diffusion coefficient may be weakly composition- and mesophase-dependent, the whole process is diffusion controlled from the earliest times for which mesophases can be distinguished. Intriguingly the effective diffusion coefficients in the mesophases are considerably greater than the self diffusion coefficients in the pure components. This suggests the existence of a sizeable free energy driving force enhancing diffusion and is reminiscent of a detailed study on the H1 phase of the $`\mathrm{C}_{12}\mathrm{E}_6`$ / water system Sallen et al. (1997). Note that the orientation of the H1 and Lα phases in the simulation is such that mutual diffusion occurs in the ‘easy’ direction: along the hexagonal rods, or parallel to the lamellar layers (we have also extracted orientation information from the eigenvectors of $`𝐌`$ which will be reported elsewhere).

We can use the effective diffusion coefficients to map elapsed times in our mesoscale simulation onto real times. The corresponding diffusion coefficients in nonionic surfactant systems are typically $`2\times 10^{10}\mathrm{m}^2\mathrm{s}^2`$ Chen et al. (2000). The lamellar repeat spacing $`L`$ can be used as a common measure of length: $`L2.5`$ DPD units in the simulation, and typically $`L5\mathrm{nm}`$ in real systems. By using $`L^2/D`$ to scale time, we conclude one DPD time step is approximately equivalent to $`50\mathrm{ns}`$. Thus the duration of the whole simulation (2000 DPD time units) represents more than $`0.1\mathrm{ms}`$ of real time, equivalent to $`10^8`$ molecular dynamics time steps. Fig. 2 shows that adiabatic diffusion controlled dissolution sets in after about $`200`$ DPD time units, equivalent to about $`10\mu \mathrm{s}`$, which is the origin of the time scale quoted in the introduction.

Our simulations, whilst addressing very directly intriguing questions concerning surfactant dissolution, are also interesting because they are simulations of phase formation kinetics starting from a highly inhomogeneous state rather than the usual approach which is to quench a homogeneous system. To our knowledge, the only comparable study in the past has been the disappearance of an interface in a vapour-liquid or binary liquid mixture after an ‘antiquench’ (sudden increase in temperature) Ma et al. (1993). It is indeed striking that we do *not* find any evidence of the delayed appearance of mesophases due to metastability, unlike that seen experimentally in temperature-jump experiments Clerc et al. (1995), nor do we find any systematic displacement of the mesophase boundaries by the strong concentration gradients present at such early times.

Our mesoscale model is also well suited to other fundamental (meso)phase studies such as temperature quenches Jury et al. (1999); Clerc et al. (1995), the effect of shear on phase boundaries, the exploration of epitaxial relationships such as between the H1 and Lα phases in the present simulations, amongst many other possibilities.


text:
# Surface modes and vortex formation in dilute Bose-Einstein condensates at finite temperatures

## Abstract

The surface mode spectrum is computed self-consistently for dilute Bose-Einstein condensates, providing the temperature dependence of the surface mode induced vortex nucleation frequency. Both the thermodynamic critical frequency for vortex stability and the nucleation frequency implied by the surface excitations increase as the critical condensation temperature is approached from below. The multipolarity of the destabilizing surface excitation decreases with increasing temperature. The computed finite-temperature critical frequencies support the experimental observations and the zero-temperature calculations for vortex nucleation.

Since it became feasible to create and observe singly quantized vortices Matthews1999b ; Madison2000a  and vortex lattices Madison2000b ; Abo2001a  in dilute Bose-Einstein condensates, a desire for quantitative understanding how vortices penetrate these unique quantum systems emerged. When the system is externally rotated above a certain critical nucleation frequency $`\mathrm{\Omega }_\nu `$, vortices enter the condensate from its boundary. However, the exact dependence of $`\mathrm{\Omega }_\nu `$ on the external variables and the experimental conditions remain the topic of continuing discussions. The reported minimum ratios between the external rotation and the harmonic trapping frequencies for which vortices have been observed thus far include approximately 0.4 Haljan2001a  and 0.7 Hodby2001a  in oblate, and 0.1 Raman2001a , 0.3 Abo2001a , and 0.7 Chevy2000a ; Madison2001a  in prolate trap geometries. Also, a number of theoretical interpretations have been proposed for explaining these frequencies Isoshima1999a ; Feder1999a ; Dalfovo2000a ; Feder2001a ; Garcia-Ripoll2001a ; Mitsubishi ; Garcia-Ripoll2001b ; Sinha2001a ; Anglin2001a ; Williams2001a ; Garcia-Ripollpreprint ; Anglinpreprint ; Muryshevpreprint ; Kramerpreprint .

First, there exists an angular frequency above which the free energy of an axially symmetric vortex state becomes lower than that of the vortex-free state. In most experiments, this thermodynamic critical frequency $`\mathrm{\Omega }_c`$, has been found to be substantially below the lowest frequency for which vortices have been observed. This may be explained in terms of an entry barrier at the surface of the condensate, preventing the transition to the lower free energy state. This, in turn, implies hysteresis effects when the rotating drive is ramped above the nucleation threshold and back again. Hence, vortices should occur in condensates at much lower rotation frequencies than those used to create them. In fact, according to the finite-temperature Popov, G1 and G2 approximations, axisymmetric singly quantized vortex states are predicted to be locally energetically (meta)stable even below $`\mathrm{\Omega }_c`$ Isoshima1999b ; Virtanen2001a  (see, however, the discussion in Ref. Virtanen2002a ).

Secondly, there exist nodeless surface excitations, or surfons, in the spectrum of the Bose-Einstein condensate that are localized near the surface of the condensate. These ‘quantum tidal waves’ Khawaja1999a  are often interpreted as being responsible for the instability of the condensate, leading to the onset of vortex nucleation as they represent a density perturbation on the condensate surface that may evolve into a quantized vortex line, thus changing the topology of the condensate. For rotationally invariant systems, an analog of the Landau criterion Landau1941  may be expressed for the surfons as

$$\mathrm{\Omega }_s=\underset{l}{\mathrm{min}}\left\{\frac{\omega _l}{l}\right\},$$
(1)
where $`\mathrm{}\omega _l`$ is the excitation energy of a surfon carrying the angular momentum $`\mathrm{}l`$. This angular frequency of the rotating perturbation corresponds to the instability of the vortex-free condensate since beyond it there emerges anomalous negative-energy excitations with high multipolarities in the spectrum of the condensate Isoshima1999a . However, the value for the nucleation frequency predicted by Eq. (1) is still too low for serving to explain the observed vortex nucleation thresholds in Refs. Chevy2000a ; Madison2001a ; Hodby2001a . Accordingly, it has been suggested that in those cases the external perturbation used to rotate the system has mainly driven the ($`l=2`$) quadrupole surface mode, for which $`\omega _l/l`$ yields a much higher critical angular frequency. Also, the bending of vortex lines and anomalous vortex core modes have been argued to explain certain aspects of the vortex nucleation process Garcia-Ripoll2001b ; Feder2001a .

In this paper, we study finite-temperature effects on vortex nucleation and the thermodynamic stability of vortices in dilute Bose-Einstein condensates. The temperature dependence of the surfon-induced vortex nucleation frequencies, as well as the thermodynamic critical frequency, are computed. The numerical results obtained are consistent with the experimental observations and support the present theoretical understanding of the vortex nucleation process based on zero-temperature field theories. Although our results are in qualitative agreement with those of Ref. Mitsubishi , there exist quantitative differences besides the system parameters, since we do not discard any low-energy quasiparticle excitations from the self-consistent sums, and we do employ the local density approximation taking properly into account also the high-lying quasiparticle states. Furthermore, due to the thermal gas, there is a slight difference between the prediction of Eq. (1) and the frequency at which there emerges anomalous surface mode(s) in the spectrum calculated in the rotating frame of reference.

The stationary generalized Gross-Pitaevskii equation

$$\mu \varphi (𝐫)=\left[_0+g|\varphi (𝐫)|^2+2g\stackrel{~}{n}(𝐫)\mathrm{\Omega }L_z\right]\varphi (𝐫)$$
(2)
for the macroscopic wavefunction $`\varphi (𝐫)`$ models the Bose-Einstein condensate at finite temperatures. Above, $`\mu `$ stands for the chemical potential, $`_0(\mathrm{}^2/2M)^2+V(𝐫)`$, and the interaction strength $`g=4\pi \mathrm{}^2a/M`$, where $`a`$ is the scattering length for atoms of mass $`M`$. The wavefunction is normalized according to $`[|\varphi (𝐫)|^2+\stackrel{~}{n}(𝐫)]d𝐫=N`$ and the system is rotated at the angular velocity $`𝛀`$, providing the $`N`$ trapped atoms with the angular momentum $`𝑳`$. Furthermore, the density of the noncondensate atoms $`\stackrel{~}{n}(𝐫)`$ is determined self-consistently via the relation

$$\stackrel{~}{n}(𝐫)=\underset{q}{}[(|u_q(𝐫)|^2+|v_q(𝐫)|^2)n_q+|v_q(𝐫)|^2],$$
(3)
where $`n_q=(e^{E_q/k_\mathrm{B}T}1)^1`$ is the Bose-Einstein distribution function and the quasiparticle amplitudes $`u_q(𝐫),v_q(𝐫)`$, and eigenenergies $`E_q`$ are obtained from the Bogoliubov equations

$`u_q(𝐫)+g\varphi ^2(𝐫)v_q(𝐫)`$ $`=`$ $`E_qu_q(𝐫),`$ (4a)
$`^{}v_q(𝐫)+g\varphi ^2(𝐫)u_q(𝐫)`$ $`=`$ $`E_qv_q(𝐫).`$ (4b)
Within the Popov approximation employed here, $`_0\mu +2g[|\varphi (𝐫)|^2+\stackrel{~}{n}(𝐫)]\mathrm{\Omega }L_z`$, and $`q`$ labels the quasiparticle states whose amplitudes must obey the normalization condition $`[|u_q(𝐫)|^2|v_q^{}(𝐫)|^2]d𝐫=\delta _{qq^{}}`$. Negative-energy quasiparticle excitations with positive norm are referred to as anomalous modes. The free energy of the system, $`F=H_{\mathrm{eff}}TS`$, is evaluated using the equation Freenote

$`F=\mu N`$ $``$ $`{\displaystyle \frac{g}{2}}{\displaystyle |\varphi (𝐫)|^4d𝐫}2g{\displaystyle \stackrel{~}{n}(𝐫)|\varphi (𝐫)|^2d𝐫}`$
$`+`$ $`{\displaystyle \underset{q}{}}E_qn_q{\displaystyle E_q|v_q(𝐫)|^2d𝐫}`$
$``$ $`k_\text{b}T{\displaystyle \underset{q}{}}[(1+n_q)\mathrm{ln}(1+n_q)n_q\mathrm{ln}n_q].`$
We shall consider cylindrically symmetric states, for which $`\varphi (𝐫)=\varphi (r)e^{im\theta }`$, the integer winding number $`m`$ determining the number of circulation quanta in the vortex. Consequently, the equations may be separated into the radial, angular and axial parts, and the surfons are uniquely labeled by their angular momentum quantum number $`l`$ (the radial and axial quantum numbers for these states equal zero). The condensate is radially confined by a harmonic trapping potential $`V(𝐫)=\frac{1}{2}M\omega _r^2r^2`$. In the axial direction we impose periodic boundary conditions. The relevant physics of this system is expected to closest resemble spherical/oblate condensate geometries, as the comparison of the critical frequencies with those in Ref. Dalfovo2000a  also suggests. In order to guarantee the accuracy of our results Numerics , we compute the discrete quasiparticle states using Eqs. (4) up to $`E_q=50\mathrm{}\omega _r`$ and employ the local density approximation Giorgini1997a ; Reidl1999a  for the remaining states. In the results presented, $`\omega _r=2\pi \times 10`$ Hz and the effective gas parameter $`Na/a_{\mathrm{ho}}90`$, where $`a_{\mathrm{ho}}=\sqrt{\mathrm{}/m\omega _r}`$ denotes the harmonic-oscillator length.

Energies per angular momentum of the surface modes with different multipolarities are plotted in Fig. 1 for a range of temperatures. The solid line is the hydrodynamic prediction $`1/\sqrt{l}`$ for the surfon dispersion relation Stringari1996a . At finite temperatures, Kohn’s theorem for parabolic confinement, according to which the dipole mode ($`l=1`$) should oscillate exactly at the trap frequency $`\omega _r`$, is not satisfied due to the neglect of the noncondensate dynamics. For all temperatures, the nucleation frequency $`\omega _l/l`$ implied by the quadrupole mode is found to be close to the value 0.7 observed in the experiments for which the quadrupolar drive may be expected to be the primary source for nucleation Chevy2000a ; Madison2001a ; Hodby2001a . This confirms the negative result found in Ref. Hodby2001a  for the significant temperature dependence of the vortex nucleation induced by a quadrupolar drive. That particular nucleation threshold is also expected to be fairly insensitive to the particle number and the geometry of the trap Dalfovo2000a .

However, the energies of the surface modes with higher multipolarities progressively increase at higher temperatures, and the lowest critical frequency $`\mathrm{\Omega }_s`$ due to the surfon excitation is found to increase approximately from 0.4 $`\omega _r`$ to 0.5 $`\omega _r`$ with increasing temperature. Also, the corresponding surface mode shifts to a lower value of angular momentum, see the inset in Fig. 1, which is a direct consequence of the shrinking of the radius of the condensate.

The computed radial distributions for the densities of the condensate, noncondensate and the destabilizing surface mode are depicted for $`\mathrm{\Omega }=0`$ in Fig. 2 for two different temperatures. These surfons are localized near the edge of the condensate, of which only a part is seen within the area of the figure. The shift in the multipolarity of the destabilizing surfon as a function of temperature is illustrated in the figure. As the temperature increases the radial extent of the condensate cloud shrinks. Hence, the angular momentum of the surfon localized in the vicinity of the condensate boundary has the larger angular momentum the lower the temperature is. Notice the noncondensate bulge, due to the mutual repulsion between the condensate and noncondensate particles, on the outer edge of the condensate.

In Fig. 3 we have plotted the lowest surfon-related ($``$) nucleation frequencies $`\mathrm{\Omega }_s`$ and two different thermodynamic critical frequencies. Also shown is the critical condensation temperature $`T_c(\mathrm{\Omega })=T_c^0(1\mathrm{\Omega }^2/\omega _r^2)^{1/3}`$ Stringari1999a  and the nucleation frequency suggested by the anomalous vortex core mode ($``$), calculated within the nonselfconsistent Bogoliubov approximation. Both the thermodynamic and the surfon critical frequencies remain nearly constant below $`T_c/2`$. At higher temperatures, especially the thermodynamic equilibrium critical frequency $`\mathrm{\Omega }_c`$ $`()`$ strongly deviates from its zero-temperature value, in agreement with the calculation of Ref. Stringari1999a . Furthermore, the critical frequency, given by the energy difference between the vortex and nonvortex states in the absence of rotation ($`\mathrm{\Omega }=0`$), is shown ($`\mathrm{}`$). Physically, this corresponds to the situation in which the noncondensate stays at rest, which has probably been the case in most of the available experiments, excluding those of Ref. Haljan2001a . The decrease of the critical frequency with increasing temperature is probably due to the neglect of the noncondensate dynamics by the Popov approximation whose predictions also for the excitation energies are known to deviate from the experimental values for temperatures $`TT_c/2`$.

We have also computed the surfon spectrum for an axisymmetric vortex state and found that the corresponding nucleation frequencies are only slightly higher than those presented in Fig. 3. This is in accordance with the fact that the angular rotation frequency needed for nucleating vortices grows in proportion to the number of vortices present in the system. Moreover, a presence of multiply quantized vortices in the harmonic trap Simula2002a  should not affect the qualitative process of vortex formation at the condensate boundary.

In conclusion, we have computed the surface mode spectrum and the thermodynamic phase diagram for vortex stability in the $`T\mathrm{\Omega }`$ plane for dilute Bose-Einstein condensates at finite temperatures. According to the current understanding, the highest observed nucleation thresholds with $`\mathrm{\Omega }/\omega _r0.7`$ can be explained in terms of the excitation of the quadrupole surface modes, whereas the lower values $`0.30.4`$ lie close to the minima of $`\omega _l/l`$ for the corresponding systems. This scenario for vortex nucleation in dilute Bose-Einstein condensates is suggested to be valid by our computations also at finite temperatures. However, the nucleation value $`0.1`$ Raman2001a  coinciding with the thermodynamic critical frequency, cannot solely be explained by the present analysis. A remedy of this observation in the favor of the surfon instability mechanism could be obtained by noting that the local fluid velocity around the small stirrer used in the experiment may well exceed the velocity of the stirrer and thus the effective rotation frequency could be locally much higher Anglin2001a .

###### Acknowledgements.

We thank the Center for Scientific Computing for computer resources, and the Academy of Finland and the Graduate School in Technical Physics for support. One of us (TPS) is grateful to the Research Council of Helsinki University of Technology for a postgraduate scholarship.