text:
# Globular Cluster Systems in the Hydra I Galaxy Cluster

## 1. Introduction

The study of globular clusters (GCs) in the center of galaxy clusters always has been of special interest due to their extraordinary rich abundance around the central (often cD) galaxy. Whether all central GCs belong to the bulge of the central galaxy or to the cluster as a whole is under discussion. Kinematic studies of GCs around NGC 1399 in the Fornax cluster (Richtler et al., this volume; Kissler-Patig et al. 1999) and M87 in the Virgo cluster (Bridges et al., this volume; Cohen & Ryzhov 1997) seem to indicate that most GCs outside a certain radius follow the gravitational potential of the galaxy cluster rather than that of the galaxy. Also the detection of GCs far outside the center of NGC 1399 (Dirsch et al., this volume) seems to confirm the existence of intra-cluster globular clusters, as predicted by some authors (e.g. West et al. 1995, and references therein). However, the process that is responsible for the formation of a very rich central GCS as well as a cD halo has not been understood yet. The main scenarios are i) a multi-phase in situ formation within a common host halo (Forbes et al. 1997), ii) stripping of GCs from nearby galaxies (Kissler-Patig et al. 1999), or iii) formation in small sub-halos and the subsequent assembly into a giant elliptical, as for example the accretion of a large number of dwarf galaxies (e.g. Côté et al. 1998; Hilker et al. 1999). In order to decide which one of the processes is dominant, further observational constraints are needed.

The availability of 8m-class telescopes has made it possible for the first time to investigate globular cluster systems in cluster of galaxies in great detail beyond the nearby Virgo and Fornax cluster. At the distance of the Hydra I cluster ($`44`$ Mpc or $`(mM)_V33.2`$ mag; Struble & Rood 1991), only a few CCD fields are needed to cover a wide range in radius from the cluster core to the outside. With good seeing conditions and moderate exposure times, globular clusters can be easily detected down to the peak of their luminosity function ($`V26`$). Moreover, dwarf galaxies as faint as the Local Group dwarf spheroidals ($`\mu _V25`$ mag/arcsec2) can be resolved and morphologically classified, and their properties can be related to those of the globular clusters. The same is true for the extended stellar halo of the central cD galaxies whose light can be followed out to faint surface brightnesses. Such a dataset enables us to study not only the GCSs of individual galaxies, but the GCS of the galaxy cluster as a whole.

In this contribution, I present first results on globular clusters in the Hydra I galaxy cluster. This cluster has a compact regular core shape and a pair of bright galaxies in the center, whose brighter component, NGC 3311, possesses a cD halo and a extraordinary rich globular cluster system (e.g. McLaughlin et al. 1995).

## 2. Observation, reduction, and candidate selection

Images of seven fields in the Hydra I cluster (see L-shape configuration in Fig.1) and one background field have been taken through $`V`$ and $`I`$ filters. The observations have been performed in a service mode run in April 2000 with UT1/FORS1 at ESO/Paranal. All exposures were taken during dark time and with a seeing between $`0.^{\prime \prime }5`$ and $`0.^{\prime \prime }7`$. The exposure times for $`V`$ and $`I`$ were $`3\times 8`$ min and $`9\times 5.5`$ min, respectively.

On the combined images, the light of the main elliptical galaxies has been subtracted by ellipse fitting before applying SExtractor for the finding and photometry of the objects. Point sources have been measured through an aperture of 2′′diameter, and then corrected for the total magnitude and photometric zero points from the ESO standard calibration plan.

All objects that have a FMHM$``$2′′and a color $`0.4<(VI)<1.6`$ mag have been defined as globular cluster candidates. The adopted reddening throughout the analysis was $`E_{VI}=0.10`$ mag, the average of the values given by Schlegel et al. (1998) and Burstein & Heiles (1984).

## 3. Globular Cluster color and spatial distribution

Some years ago it was thought that there is a lack of metal-poor, blue GCs in central galaxy NGC 3311. Using Washington photometry, only red GCs with a mean metallicity of \[Fe/H\]=$`0.31`$ dex were found (Secker et al. 1995). Only recently, this surprising result has been overruled by HST photometry, which shows a normal color distribution of the GCS, comparable to that of other giant ellipticals (Brodie el al. 2000). This has also been confirmed by our analysis (see Fig.2). For the color distribution in Fig.2, only point sources with $`V<25.5`$ mag and a photometric error less than 0.1 mag have been selected. Applying the same selection criteria for the control field, the color distribution of background sources has been subtracted statistically. NGC 3311 as well as NGC 3309 possess a blue globular cluster sub-population with a well defined peak at about $`(VI)=0.93`$ mag (\[Fe/H\]$`1.5`$ dex). The red sub-population has a broader distribution peaking around $`(VI)=1.08`$ mag. Most probably, this reflects a superposition of several red sub-populations of different metallicities. The existence of an extended blue tail in the color distributions points either to a population of very metal-poor globular clusters (\[Fe/H\]¡$`2.0`$ dex) or a population of young clusters of a rather recent star formation event (that might be related to the central dust lane in NGC 3311).

In the outer halo of NGC 3311 (beyond 16 kpc from the galaxy center) blue globular clusters start to dominate the distribution. The lower right panel of Fig.2 shows that (mainly blue) globular clusters exist out to very large distances from the cluster center ($`<250`$ kpc). In this plot all globular cluster systems of the major galaxies in Hydra have been excluded.

The dominance of blue globular clusters in the halo and the large extent of the intra-cluster GCS is also reflected in the radial distribution of the different sub-populations. Figure 3 shows the high concentration of red clusters ($`1.0<(VI)<1.2`$) within a radius of 10. Interestingly, the distribution of dwarf galaxy candidates (see Fig.1) closely follows the distribution of blue globular clusters ($`0.7<(VI)<0.9`$), except a possible deficiency in the very center of the cluster.

## 4. Future Analysis

The study of the color and spatial distribution has shown that the investigation of globular clusters in galaxy clusters should not only be restricted to the analysis of GCSs around individual member galaxies, but all the intra-cluster GCs has to be taken into account. The next steps for further analyses of the presented data are:
$``$ Constructing a density and color map for the sub-populations of the intra-cluster GCS in order to define its center and compare it to that of the central galaxy and the X-ray gas halo.
$``$ Modelling the cD halo light and study the specific frequency as a function of centro-cluster distance.
$``$ Studying the individual GCSs for all member galaxies down to the dwarf galaxy regime (colors, density profiles, specific frequency, etc.).
$``$ Confirming the membership of dwarf galaxies and bright compact objects by follow-up spectroscopy (Magellan I data in hand).

### Acknowledgments.

I wish to thank L. infante and F. Barrientos for their contributions in this project and FONDECyT for its support through ‘Proyecto FONDECYT 3980032’.

## References

Brodie, J.P., Larsen, S.S., & Kissler-Patig, M. 2000, ApJ, 543, L19

Burstein, D., & Heiles, C., 1984 ApJS, 54, 33

Cohen, J.G., & Ryzhov, A. 1997, ApJ, 486, 230

Côté, P., Marzke, R.O., & West, M.J. 1998, ApJ, 501, 554

Forbes, D.A., Brodie, J.P., & Grillmair, C.J. 1997, AJ, 113, 1652

Hilker, M., Infante, L., & Richtler, T. 1999, A&AS, 138, 55

Kissler-Patig, M., Grillmair, C.J., Meylan, G., Brodie, J.P., Minniti, D., & Goudfrooij, P. 1999, AJ, 117, 1206

McLaughlin, D.E., Secker, J., Harris, W.E., & Geisler, D. 1995, AJ, 109, 1033

Schlegel, D.A., Finkbeiner, D.P., & Davis, M. 1998, ApJ, 500, 525

Secker, J., Geisler D., McLaughlin D.E., & Harris W.E. 1995, AJ, 109, 1019

Struble, M.F., & Rood, H.J. 1991, ApJS, 77, 363

West, M.J., Côté, P., Jones, C., Forman, W., & Marzke, R.O. 1995, ApJ, 453, L77

## Discussion

Grillmair:  Has anyone ever found an offset between centers of globular cluster systems and the centers of their host galaxies/galaxy clusters?

Hilker:  No, but our data set enables us to examine this for the first time. It would ber of special interest to investigate this issue in clusters that show an X-ray halo whose center is displaced with respect to the position of the central galaxy.

Gallagher:  Have you compared the radial distributions of blue GCs in Hydra with that of the dEs?

Hilker:  First analysis seems to indicate that the spatial distribution of both populations are closely related (see Fig.3). However, the membership of the individual dwarf galaxy candidates has to be confirmed (by our follow-up spectroscopy for example) and their globular cluster systems have to be identified, before further interpretations can be made.


format: text

text:
# Dilute Bose gas in a quasi two-dimensional trap

## I Introduction

Over the last decade or so, there have been a number of theoretical studies of Bose-Einstein condensation (BEC) in inhomogeneous low-dimensional Bose systems. In principle, such low-dimensional systems can be created in the laboratory through a suitable manipulation of the potential used to trap the ultra-cold atoms. In this paper, we will be primarily interested in the situation for which the confinement in the $`z`$-direction is much steeper than the in-plane confinement, so that the system can be viewed as being (quasi) two-dimensional (2D) in the sense that the motion of the atoms in the $`z`$-direction is “frozen out”.

It is well-known that a homogeneous 2D Bose gas has no long range order (at $`T0`$), and so a BEC transition at any (non-zero) finite temperature cannot occur. However, in a 2D harmonically confined system with oscillator frequency $`\omega _0`$, it is easily shown that an ideal Bose gas has a critical BEC temperature

$$T_c^{(0)}=(6/\pi ^2)^{1/2}N^{1/2}\mathrm{}\omega _0,$$
(1)
which, for $`N1`$, is much larger than the oscillator gap $`\mathrm{}\omega _0`$. It is interesting to consider the limit $`N\mathrm{}`$, $`\omega _00`$, such that $`N^{1/2}\omega _0=`$ constant. In this situation, the critical BEC transition temperature $`T_c^{(0)}`$ remains the same. Note that this is not the usual thermodynamic limit, which in 2D would demand that $`N\omega _0=`$ constant, resulting in no BEC. We have previously shown that a self-consistent Thomas-Fermi (SCTF) model (i.e., with no condensate density), has a solution at all temperatures provided the mean-field contribution of a repulsive zero range interaction is included, no matter how weak. This result contradicts a previous claim in the literature that no self-consistent solution of the TF equation in a trap is possible below a certain critical temperature $`T_c`$. (It is worth pointing out that in Ref. , the absence of a solution below $`T_c`$ was interpereted as a phase transition in the sytem, but not a phase transition to the BEC state.) The absence of a phase transition for the repulsive case was also confirmed analytically by examining the branch-points of the density in the complex fugacity plane. Note, however, that the SCTF method cannot address the important issue of phase fluctuations in the condensate. It was also found that the such a system could be described by noninteracting haldons obeying the generalized fractional exclusion statistics (FES). This has been recently confirmed by Hansson et al. .

For a dilute Bose gas in a three dimensional trap, where the average distance between the atoms is much larger than the scattering length, the Gross-Pitaevskii (GP) approach of using a zero-range pseudo-potential with a strength given by $`4\pi \mathrm{}^2a/m`$, (where $`a`$ is the s-wave scattering length, and m is the mass of the atom) is very successful. For a strictly 2D system however, the strength has a logarithmic energy dependence. Recently, Petrov et al.  have deduced the strength of the quasi-2D system by considering free two-particle scattering in the plane of the gas, and strong harmonic confinement in the $`z`$direction (i.e., $`\omega _0\omega _z`$). The effective two dimensional pseudo-potential may be written as a zero-range momentum-dependent interaction with a strength given by

$$\frac{2\pi \mathrm{}^2}{m}g=\frac{2\sqrt{2\pi }\mathrm{}^2}{m}[\mathrm{}_z/a\frac{1}{\sqrt{2\pi }}\mathrm{ln}(\pi p_{12}^2\mathrm{}_z^2)]^1,$$
(2)
where $`𝐩_{\mathrm{𝟏𝟐}}=(𝐩_\mathrm{𝟏}𝐩_\mathrm{𝟐})/\mathrm{𝟐}`$ is the relative momentum between two bosons. In the above equation, $`m`$ is the mass of the atom, $`a`$ the s-wave scattering length, and $`\mathrm{}_z=\sqrt{\mathrm{}/m\omega _z}`$. Note that the kinetic energy of relative motion, $`E=p_{12}^2/m`$, reduces to twice the single-particle energy in the frame where the CM momentum $`𝐏=(𝐩_1+𝐩_2)/2=0`$. To simplify actual calculations in the many-body problem with this coupling strength, the momentum-dependence in Eq.(2) is replaced  by a density dependence by putting $`p_{12}^2=mE=2m|\mu |`$, where $`\mu `$ is the chemical potential. For the bosons in a harmonic oscillator trap, the spatial density is inhomogeneous, and the local chemical potential $`\mu (r)=\mu \frac{1}{2}m\omega _0^2r^2`$ will be used. For the situation when $`\mathrm{}_za`$, the first term on the right-hand side (RHS) of Eq. (2) dominates over the second, and the logarithmic term may be neglected. This yields a “constant” interaction strength,

$$\frac{2\pi \mathrm{}^2}{m}\stackrel{~}{g}=\frac{2\sqrt{2\pi }\mathrm{}^2}{m}a/\mathrm{}_z,$$
(3)
which agrees with our definition for the coupling constant in an earlier work.

In this paper, we perform finite-temperature SCTF calculations with the above density dependent strength factor in the pseudo-potential. Both the repulsive, and the attractive cases are considered for realistic values of the physical parameters. For the repulsive case ($`g>0`$), the earlier conclusion that there is no BEC at any temperature remains unaltered. For the attractive case ($`g<0`$) on the other hand, we find compelling numerical and analytical evidence suggesting the existence of a phase transition at a $`T^{}`$ which is greater than the critical temperature for BEC given by Eq. (1). The outline for the rest of the paper is as follows. In Sec. II, we present numerical results of the SCTF calculations for both $`\stackrel{~}{g}`$ and $`g`$ as defined above, as well as some analytical arguments explaining the origin of the critical temperature $`T^{}`$. In Sec. III we provide additional support for our numerical results by analytically investigating the thermodynamic properties of the system at zero and finite-temperatures. Finally, in Sec. IV, we present our concluding remarks.

## II Thomas-Fermi Calculation

### A Repulsive Interaction

We consider bosonic atoms in a two-dimensional harmonic oscillator trap above the critical temperature, if any. In this subsection, a repulsive interaction is used, while in the next subsection the interaction is taken to be attractive. In the SCTF method, the inhomogeneous number density of the interacting bosons is given by

$$n(r)=\frac{d^2p/(2\pi \mathrm{}^2)}{[\mathrm{exp}[(\frac{p^2}{2m}+\frac{1}{2}m\omega _0^2r^2+U(n(r))\mu )\beta ]1]},$$
(4)
where $`U(n(r))=\frac{2\pi \mathrm{}^2}{m}\stackrel{~}{g}n(r)`$ is the mean-field potential generated by the zero-range two-body pseudo-potential with the constant strength $`\stackrel{~}{g}`$, and $`\mu `$ is the chemical potential that serves to fix the total number of atoms in the trap,

$$N=n(𝐫)d^2𝐫.$$
(5)
To obtain $`n(r)`$ for a fixed $`N`$, Eqs. (4) and (5) have to be solved self-consistently. (In our numerical work, we have scaled all lengths by the characteristic in-plane oscillator length $`\mathrm{}_0=\sqrt{\mathrm{}/m\omega _0}`$, and all energies by $`\mathrm{}\omega _0`$.) The case for $`\stackrel{~}{g}>0`$, has already been addressed in an earlier paper, where it was shown that the self-consistent solution for $`n(r)`$ could be obtained for any $`T0`$, indicating that there is no BEC. We first want to check if the same conclusion holds when $`\stackrel{~}{g}`$ is replaced by the more realistic momentum dependent $`g`$, as given by Eq. (2). As noted before, we replace the $`p_{12}`$dependence in $`g`$ by an effective density dependence. Under this assumption, the $`p`$integration in (4) can be done analytically. At zero-temperature, the result is given by

$$\frac{2\pi \mathrm{}^2}{m}gn(r)=(\mu \frac{1}{2}m\omega _0^2r^2).$$
(6)
The RHS of Eq. (6) is precisely the local chemical potential $`\mu (r)`$. Hence we replace the $`p_{12}^2`$ on the RHS of Eq. (2) by $`(2\pi \mathrm{})^2gn(r)`$ and $`g`$ itself is determined self-consistently by Eqs. (5) and (6). At finite-temperatures, we continue to follow the same procedure, but in this case $`g(r)`$ and $`n(r)`$ must be found self-consistently from Eqs. (4) and (5). Note that Petrov et al.  in one of their footnotes recommend an identical prescription. Following this method, we still find that for a fixed boson number $`N`$, Eqs. (4) and (5) may be solved self-consistently for all $`T0`$, indicating that there is no phase transition for a repulsive interaction. The behavior of the chemical potential $`\mu `$ as a function of $`T/T_c^{(0)}`$ is shown in Fig. 1. The lowest solid curve, for noninteracting bosons, shows the usual discontinuity in $`\mu `$ at $`T/T_c^{(0)}=1`$. The next higher pair of closely spaced curves show the results for $`N=10^4`$ Rb87 atoms in a trap with with $`\mathrm{}\omega _0=0.6`$ nK, $`\omega _z=10^3\omega _0`$, and $`a=5.8`$ nm. The solid curve shows the result for $`\stackrel{~}{g}=0.05`$, obtained from the above physical parameters, while the slightly lower dashed curve is for the corresponding $`g`$ given by Eq. (2). We note that the difference between the two is negligible for the parameters chosen above, which suggests that FES is still obeyed to a good approximation for relatively small values of $`g`$. In the same figure, the top three curves are for an artificially boosted interaction, with the continuous curve given by $`\stackrel{~}{g}=0.25`$. Note that this could in principle be achieved in a variety of ways by altering both $`\mathrm{}_z`$ and $`a`$. Two limiting cases are given by keeping $`\mathrm{}_z`$ constant and decreasing $`a`$ by a factor of five, or keeping $`a`$ fixed and increasing $`\mathrm{}_z`$ by a factor of five. Although both procedures yield the same $`\stackrel{~}{g}`$, we see from Eq. (2) that $`g`$ gets affected more when $`\mathrm{}_z`$ alone is altered. In Fig. 1, the lowest (dashed) curve is for this case, while the next higher one is for the case where only $`a`$ is decreased by a factor of five. We see that in all cases with a repulsive interaction, the main conclusion that there is no strict BEC continues to hold.

Even though there is no BEC in the sense that the SCTF solution is obtained without paying any special attention to the lowest single-particle quantum state, it is found that there is a significant enhancement in the peak of the density distribution as the temperature is reduced below $`T_c^{(0)}`$, the critical temperature for the noninteracting case. This is shown in Fig. 2 for $`\stackrel{~}{g}=0.25`$ (continuous curves) at $`T/T_c^{(0)}=0.8`$ and $`1.1`$. When the density-dependent version of the coupling $`g`$ is used, the enhancement is even more pronounced. The reason for this behavior is clear from the inset, where the effective $`g(r)`$, as obtained from the self-consistent solution, is shown as a function of $`r`$. We note that for the density-dependent case, the interaction $`g(r)`$ is considerably weaker than the constant $`\stackrel{~}{g}`$, thus causing less suppression of the central density.

### B Attractive Interaction

We first take the example of Li7 atoms for which the s-wave scattering length is $`a=1.45`$ nm. For $`\mathrm{}\omega _0=7`$ nK (the corresponding oscillator angular frequency $`\omega _0`$ in the plane is $`2\pi \times 145`$ rads/s ) and $`\omega _z=10^3\omega _0`$, the coupling constant $`\stackrel{~}{g}=0.012`$. On the other hand, the interaction strength could be considerably larger for Cs135 atoms, which have $`a=60`$ nm. In Fig. 3, we display for $`N=10^4`$ atoms of Li7 the behavior of the chemical potential $`\mu `$ as a function of $`T`$ obtained self-consistently. No solution is found for $`TT^{}`$, irrespective of whether one uses $`\stackrel{~}{g}`$ or $`g`$. The highest lying solid curve is for $`\stackrel{~}{g}=0.012`$, and difference between this and the density-dependent form of $`g`$ for this weak interaction is too small to be shown in the same scale. In the same figure, we also display the effect of boosting the attractive interaction strength nearly ten-fold to $`\stackrel{~}{g}=0.1`$. Not only does the temperature $`T^{}`$ increases significantly, but it is also possible now to differentiate the curves for the cases where the density-dependent form of $`g`$ is used rather than $`\stackrel{~}{g}`$. Furthermore, when the boosting of the coupling $`g`$ is done by altering $`\mathrm{}_z`$ rather than $`a`$, $`T^{}`$ is larger. Nevertheless, the behavior of the chemical potential with temperature, for $`T>T^{}`$ does not alter significantly with the strength of the attractive interaction. Comparing Fig. 3 with the corresponding repulsive cases (see Fig. 1), we note that $`\mu `$ showed a much larger variation, but only in the range $`TT_c^{(0)}`$, which is in any case inaccessible for the attractive interactions.

Some insight as to why no self-consistent solution for attractive interaction strength $`\stackrel{~}{g}`$ may be found by performing the $`p`$integration in Eq. (4) analytically (valid for either sign of $`\stackrel{~}{g}`$ ). The result is given by

$$(\mu \frac{1}{2}m\omega _0^2r^2)=\frac{2\pi \mathrm{}^2}{m}\stackrel{~}{g}n(r)+\frac{1}{\beta }\mathrm{ln}[1\mathrm{exp}(2\pi \mathrm{}^2n(r)\beta /m)].$$
(7)
For a negative $`\stackrel{~}{g}`$, both terms on the RHS are negative, and therefore $`\mu `$, which is independent of $`r`$, must be negative at all temperatures. For this case, if we take the above expression literally at $`T=0`$, the density is given by

$$n(r)=\frac{m}{2\pi |\stackrel{~}{g}|\mathrm{}^2}(|\mu |+\frac{1}{2}m\omega _0^2r^2).$$
(8)
This solution for $`n(r)`$ is rejected since it keeps increasing with $`r`$ ! We demand a normalized density $`n(r)`$ that is a monotonically decreasing function of $`r`$, which is only found for $`T>T^{}`$. By differentiating Eq. (7) with respect to $`r`$, we find that for $`dn/dr0`$ for all values of $`r`$, we must have

$$[\mathrm{exp}(\frac{2\pi \mathrm{}^2}{m}n(r)\beta )1]^1|\stackrel{~}{g}|0.$$
(9)
This may be satisfied only for $`TT^{}`$, which is related to the central density at this temperature $`n(0)`$ by the relation

$$T^{}=\frac{\frac{2\pi \mathrm{}^2}{m}n(0)}{\mathrm{ln}(1+1/|\stackrel{~}{g}|)}.$$
(10)
It is interesting that the above equation may also be derived using the the generalized exclusion statistics of Haldane. It has been shown  that when the interaction strength of the zero-range pseudo-potential is taken to be a constant (which we denote by $`\stackrel{~}{g}`$ in this paper), the interacting boson system may be mapped on to a noninteracting system of haldons, whose statistical occupancy factor of a single particle state $`i`$ with energy $`ϵ_i`$ is given by

$$\eta _i=(w_i+\stackrel{~}{g})^1,$$
(11)
which must be $`0`$. Note that $`w_i`$ obeys the relation

$$w_i^{\stackrel{~}{g}}(1+w_i)^{1\stackrel{~}{g}}=\mathrm{exp}[(ϵ_i\mu )\beta ].$$

In what follows, we go to the limit $`\omega _00`$, $`N\mathrm{}`$, such that $`N^{1/2}\omega _0=`$ constant. We then identify the central density $`n(0)`$ of the harmonic oscillator potential as the uniform density in this limit. The number density $`n(0)`$ of these haldons at the center of the harmonic oscillator potential is given by

$$n(0)=\frac{m}{2\pi \mathrm{}^2\beta }\mathrm{ln}\frac{1+w_0}{w_0},$$
(12)
where $`w_0`$ defines the occupancy of the lowest quantum state through Eq. (11). We now use these relations for the attractive interaction strength $`\stackrel{~}{g}=|\stackrel{~}{g}|`$. Note from Eq. (11) that for the occupancy $`\eta _0`$ to be positive definite, $`w_0|\stackrel{~}{g}|`$, and that the occupancy tends to infinity when $`w_0=|\stackrel{~}{g}|`$. Substituting this limiting value of $`w_0`$ in Eq. (12) at $`T^{}`$, we immediately obtain Eq. (10). We should regard $`T^{}`$ as the limiting temperature below which the system cannot be described as a collection of ideal haldons.

In Fig. 4, the self-consistent spatial density $`n(r)`$ is shown for $`N=10^4`$ Li7 atoms in the quasi 2D trap, with $`\stackrel{~}{g}=0.012`$. For this very weak interaction, the corresponding curve for the density-dependent coupling strength $`g`$ cannot be differentiated on this scale. We see that not only is there there a marked enhancement in the central density when $`T`$ is decreased from $`1.1`$ $`T^{}`$ to $`T^{}`$, but the shape of the density distribution is also different. We cannot, of course, get any solution for $`T<T^{}`$. It may be instructive to find a relation between $`T^{}`$ and $`N`$ analogous to Eq. (1), which will of course also depend on the interaction strength $`\stackrel{~}{g}`$. Empirically, we find that this is given by

$$T^{}=\gamma (|\stackrel{~}{g}|)\left(\frac{N}{\zeta (2)}\right)^{1/2}\mathrm{}\omega _0,$$
(13)
where $`\gamma (|\stackrel{~}{g}|)=[1+\frac{6}{5}(|\stackrel{~}{g}|)^{2/3}]`$. This relation is quite accurate even for small $`N`$ and large $`\stackrel{~}{g}`$, as depicted in Fig. 5. In this figure, we show the numerically calculated points by circles, and the continuous curves are derived from Eq. (13). Although the solid curves fitting the data have been extended down to $`N=0`$ at $`T=0`$, we were not able to numerically access stable values of $`N`$ for $`T^{}8`$, and so Eq. (13) is only used to fit the data down to a finite temperature $`T^{}8`$. Also note from Eqs. (1) and (13) that $`T^{}/T_c^{(0)}`$ is independent of $`N`$, and therefore we can draw a “universal” curve of $`\stackrel{~}{g}`$ versus $`T^{}/T_c^{(0)}`$, which is shown by the inset in Fig. 5. The normal bosonic phase is above the continuous curve shown in this inset, but the SCTF method used cannot access the region below.

## III Analytical considerations

In this section, we wish to address the question that naturally arises from the numerical work presented in Sec. IIB; namely, are our numerical calculations suggesting that there is phase transition in the quasi-2D trap for attractive interactions ? Since we are unable to answer this question within the framework of the SCTF model, we now turn our attention to various analytical approaches which may aid in our assessment of the numerical results presented above.

### A Zero-Temperature Variational Calculation

Let us consider the 2D analogue of the GP energy functional, namely

$$E=d^2r\left[\frac{\mathrm{}^2}{2m}|\psi |^2+\frac{m}{2}r^2\omega _0^2|\psi |^2+\frac{1}{2}\frac{2\pi \mathrm{}^2}{m}\stackrel{~}{g}|\psi |^4\right],$$
(14)
where $`\stackrel{~}{g}`$ is given by Eq. (3), and may be of either sign. In Eq. (14), the function $`\psi (r)`$ represents the wave function of a system for which there is macroscopic occupation of the lowest quantum state (i.e., all $`N`$ bosons are in the ground state):

$$\psi (r)=\frac{1}{b}\sqrt{\frac{N}{\pi }}\mathrm{exp}(\frac{r^2}{b^2}),$$
(15)
where $`b`$ is taken to be a variational parameter. We find that the energy is a minimum for

$$b=\sqrt{\frac{\mathrm{}}{m\omega _0}}(1+N\stackrel{~}{g})^{1/4}.$$
(16)
For an attractive interaction, $`\stackrel{~}{g}=|\stackrel{~}{g}|`$, Eq. (16) implies that the condensate collapses for

$$N|\stackrel{~}{g}|1.$$

This rough estimate is in good agreement with the numerical estimate of $`N|\stackrel{~}{g}|0.94`$ by Adhikari  which was carefully obtained through a self-consistent solution of the 2D GP equations. For a weakly interacting gas like Li7 with $`\stackrel{~}{g}=0.012`$, we find that the condensate is unstable for $`N>83`$.

### B Finite-Temperature Thermodynamics

This subsection, deals with the study of the behavior of the branch points of the quantity $`u=2\pi \mathrm{}^2\beta n(r)/m`$ in the complex fugacity plane $`Z=\mathrm{exp}(\beta \mu )`$ as a function of the coupling constant $`\stackrel{~}{g}`$. This is analogous to the analysis that was made by Sutherland while studying the thermodynamics of a one-dimensional gas interacting with an inverse square potential. The important point to be taken from Sutherland’s work is that a branch point of $`u`$ on the positive real axis of $`Z`$ indicates a phase transition in the system.

We begin by noting that Eq. (7) can be rewritten as

$$Ze^{\beta m\omega _0^2r^2/2}=2e^{(\stackrel{~}{g}1/2)}\mathrm{sinh}\frac{u}{2}.$$
(17)
As before, going to the limit $`N\mathrm{}`$, $`\omega _00`$ with $`N^{1/2}\omega _0=`$ constant, Eq. (17) reduces to the $`r`$-independent relation

$$Z=e^{\stackrel{~}{g}u}e^{(\stackrel{~}{g}1)u}.$$
(18)
We note that for $`\stackrel{~}{g}=0`$, the above equation gives $`u=\mathrm{ln}(1Z)`$, implying that $`u`$ has a branch point at $`Z=1`$, i.e. for $`\mu =0`$. This is precisely where BEC takes place for a noninteracting gas in a harmonic trap. For $`\stackrel{~}{g}=1/2`$, the branch points of $`u`$ are at $`Z=\pm 2i`$, and for $`\stackrel{~}{g}=1`$, the branch point is at $`Z=1`$. For arbitrary values of $`\stackrel{~}{g}`$, the branch points of $`u`$ can be obtained by examining the analytical structure of the solutions to the equation ($`x\mathrm{exp}(u)`$)

$$x^{\stackrel{~}{g}}x^{\stackrel{~}{g}1}Z=0.$$
(19)
Following Sutherland’s prescription, we obtain the branch points of $`u`$, denoted by $`Z_0`$, for all $`\stackrel{~}{g}0`$:
(i) $`\underset{¯}{\stackrel{~}{g}1}`$

$$Z_0=\frac{1}{\sqrt{\stackrel{~}{g}(\stackrel{~}{g}1)}}\left(\frac{\stackrel{~}{g}1}{\stackrel{~}{g}}\right)^{(\stackrel{~}{g}1/2)},$$
(20)
(ii) $`\underset{¯}{1>\stackrel{~}{g}0}`$

$$Z_0=e^{\pm \pi i\stackrel{~}{g}}\frac{1}{\sqrt{\stackrel{~}{g}(1\stackrel{~}{g})}}\left(\frac{1\stackrel{~}{g}}{\stackrel{~}{g}}\right)^{(\stackrel{~}{g}1/2)}.$$
(21)
In Fig. 6, the loop and the line along the negative real axis $`1Z0`$ in the complex fugacity plane shows the branch points of $`u`$ for positive values of $`\stackrel{~}{g}`$ given by Eqs. (20) and (21). Since there is no branch point on the real fugacity axis for $`\stackrel{~}{g}>0`$, there is no phase transition, in agreement with our numerical calculations.

The situation for $`\stackrel{~}{g}<0`$ can be addressed by either solving Eq. (19) directly, or by noticing that the transformation

$$\stackrel{~}{g}\stackrel{~}{g}+1,uu,ZZ,$$
(22)
leaves Eq. (18) unchanged. More specifically, under Eq. (22), the branch points of $`u`$ on the negative $`Z`$axis between $`Z_0[1,0]`$ for $`1\stackrel{~}{g}<\mathrm{}`$ are now shifted to the positive side of the axis, namely, $`Z_0[0,1]`$ for all $`\mathrm{}<\stackrel{~}{g}0`$. The branch points for $`0<\stackrel{~}{g}<1`$ are all mapped back onto Eq. (21), thereby leaving the loop structure in the complex $`Z`$-plane invariant under (22) (see Fig. 6). In other words, for every $`\stackrel{~}{g}<0`$ we have an associated branch point on the positive real $`Z`$-axis, suggesting that there is a phase transition in the system for attractive interactions. It is important to note, however, that the above analysis can offer us no insight into the precise nature of the phase transition; is it BEC ? Is it the Kosterlitz-Thouless transition  ? Is it the collapse of the gas ?

### C Random-Phase Approximation

In this subsection, we will try to address the nature of the phase transition discussed above by applying the well-known random-phase approximation (RPA) to investigate the stability of the trapped 2D Bose gas with attractive interactions. Our approach parallels the recent investigation by Meuller and Baym, with the main difference being that we confine ourselves to strictly two-dimensions. The basic idea behind their method is to identify the collapse of the gas with a divergence in the density-response function $`\chi (k,\omega ;\beta )`$, which is a measure of the response of the system to some external probe with wave vector $`k`$ and frequency $`\omega `$. As in Ref. , $`\chi `$ is calculated within the local-density approximation (LDA), which amounts to replacing the response function of the inhomogeneous system by that of a system of uniform density $`n`$. In order to make contact with the trapped gas, the uniform density $`n`$ is then taken to represent the central density of the trap $`n(r=0)`$.

For a confined Bose gas, the collapse is associated with an instability of the lowest energy “breathing” mode of the system. In 2D, the breathing mode is given by

$$\delta \rho (r)\left(\frac{r^2}{\mathrm{}_0^2}1\right)\mathrm{exp}(r^2/\mathrm{}_0^2),$$

where $`r`$ is the 2D radial coordinate, and $`\mathrm{}_0=\sqrt{\mathrm{}/m\omega _0}`$ is the characteristic size of the cloud. The Fourier transform of $`\delta \rho (r)`$ has a maximum at $`k=2/\mathrm{}_0`$, and so we look for an instability at this wave vector. Owing to the fact the breathing mode has a vanishing frequency, the line of collapse is characterized by

$$\chi (k=2/\mathrm{}_0,\omega =0;T^0)^1=0,$$
(23)
where $`T^0`$ denotes the temperature for which the gas undergoes collapse. Note that $`T^0`$ need not coincide with the BEC transition temperature $`T_c`$ of the interacting gas. Since our SCTF calculations deal only with the normal phase of the Bose gas, we will restrict our attention to the noncondensate contribution to $`\chi (k,\omega ;T)`$ above, and present a more detailed analysis of the RPA in 2D attractive Bose gases elsewhere. Within the RPA, the response function for the interacting Bose gas has the structure

$$\chi ^{\mathrm{RPA}}(k,\omega ;T)=\frac{1}{(1\stackrel{~}{g}\chi _0^n)},$$
(24)
where $`\chi _0^n`$ is the “bare” noncondensate response of the noninteracting system. In 2D, we have

$$\chi _0^n(k,\omega ;T)=\frac{1}{(2\pi )^2}d^2𝐪\frac{f(qk/2)f(q+k/2)}{\mathrm{}\omega (\epsilon (q+k/2)\epsilon (qk/2))}.$$
(25)
In Eq. (25), $`f(q)`$ is the Bose distribution function, and $`\epsilon (k)=\mathrm{}^2k^2/2m`$ is the free particle kinetic energy. After integrating out the angular dependence, we have

$$\chi _0^n(k,0;T)=\frac{8\pi m}{(2\pi \mathrm{})^2k}𝑑qf(q)\frac{q}{\sqrt{k^24q^2}}.$$
(26)
Using the above expression for $`\chi _0^n`$ we can obtain the line of collapse by setting $`k=2/\mathrm{}_0`$ and solving for the roots of the denominator in Eq. (24). This procedure yields the line of collapse in terms of $`\mu `$ and $`T_{\mathrm{RPA}}^0`$. Following the arguments made in Ref. , we can then relate $`\mu `$ to the total number of normal phase atoms, $`N`$, in the harmonically confined gas through the relation

$`N`$ $`=`$ $`\left({\displaystyle \frac{k_BT}{\mathrm{}\omega _0}}\right)^2g_2(e^{\beta \mu }),`$ (27)
where $`g_\nu (z)=_pz^p/p^\nu `$ is the polylogarithm function. Note that evaluating Eq. (27) at $`\mu =0`$ defines the line of condensation for a noninteracting 2D Bose gas. (This result can easily be derived by integrating Eq. (4) with $`\stackrel{~}{g}=0`$ over all space.)

The results of our calculations for $`\stackrel{~}{g}=0.012,0.1,0.2`$ are shown in Fig. 7. For comparison, we have also included as dashed lines in the figure the $`\stackrel{~}{g}=0.1,0.2`$ SCTF results from Sec. IIB (see also Fig. 5). This figure suggests that the line of collapse derived from the RPA has the same functional relationship as the empirical formula given by Eq. (13) in Sec. IIB. In order to confirm this observation, we recall that the relationship used to link the uniform gas results to the harmonically confined gas, viz., Eq. (27), is strictly valid only for a noninteracting system. As we saw in Sec. IIB, the attractive interactions have the effect of “renormalizing” the ideal $`N`$ vs. $`T`$ dependence (see Eq. (13)), which will result in different curves for the RPA and SCTF data. With this in mind, we rescale Eq. (27) by $`1/\gamma (|\stackrel{~}{g}|)`$, and comapre the results with the SCTF data. After this procedure, we find that the SCTF and RPA curves are almost indistinguishable on the scale of Fig. 7. In this sense, we can identify $`T_{\mathrm{RPA}}^0`$ with the $`T^{}`$ of Sec. IIB. This analysis leads us to believe that the phase-boundary generated from the SCTF in Sec. IIB is in fact the line of collapse of the Bose gas as one approaches the instability from above the BEC transition temperature $`T_c`$.

To close this section, we mention that we have also performed RPA calculations for the 3D interacting Li7 Bose gas with attractive interactions (exactly as in Ref.  but without exchange) and found that the 3D gas is more stable than its 2D counterpart; that is, for the same temperature and particle number, the 3D system is always in the stable uncollapsed phase relative to the 2D gas. To illustrate this, we have included as an inset to Fig. 7, the lines of collapse (solid curve) and condensation (dashed curve) for the 3D system, focusing only on the $`N`$ and $`T`$ regimes that are relevant for a comparison with the corresponding 2D system.

## IV Summary

One of the first results of this paper was a confirmation of the fact that, at least within the SCTF method, there is no strict BEC in a quasi-2D trap for atoms with repulsive interactions, even if the more “realistic” momentum-dependent pseudo-potential recently derived by Petrov et al. , and given by Eq. (2) is used. Indeed, our results indicate that the much simpler density-independent coupling constant $`\stackrel{~}{g}`$ (Eq. (3)) is quite adequate for numerical work on 2D Bose systems. This observation provides further support to the conclusion obtained earlier  using a momentum-independent coupling constant $`\stackrel{~}{g}`$.

In the case of attractive interactions, we have found strong numerical and analytical evidence for the existence of a phase transition in the system. Since the qualitative features of this phase transition are not strongly influenced by using a more complicated momentum-dependent pseudo-potential, it is again sufficient to use the simpler $`\stackrel{~}{g}`$ in numerical work. Our numerical investigations reveal that the phase transition is signaled by a “critical” temperature $`T^{}>T_c^{(0)}`$ below which there is no self-consistent solution for the normal ground state of the system. From the point of view of FES, the $`T^{}`$ was interpreted as the limiting temperature below which the system cannot be described as a collection of noninteracting haldons.

In order to further our understanding of the numerical results, we considered several analytical approaches. This began in Sec. IIIA with a zero-temperature 2D variational calculation of the GP energy functional. Surprisingly, this simple analysis revealed that if there is a BEC for attractive interactions in 2D, then the region of stability for the condensate (i.e., prior to collapse) is severely restricted by the condition $`N|\stackrel{~}{g}|1`$. For Li7, this implies an instability for $`N>83`$. The $`T=0`$ calculation was followed in Sec. IIIB by an investigation of the branch points of $`u`$ in the complex fugacity plane. We found that for every $`\stackrel{~}{g}<0`$, there is an associated branch point of $`u`$ on the positive real $`Z`$-axis. This result strongly supports our numerical calculations which indicate that there is a phase transition taking place in the system for attractive interactions.

Finally, in Sec. IIIC, an RPA-type calculation was used to investigate the line of collapse for a 2D Bose gas with attractive interactions. The purpose of this study was to try and determine the nature of the phase transition suggested by the analysis in Sec. IIIB. We found that the line of instability in the RPA calculation had the same $`N`$ vs. $`T`$ dependence as the SCTF curves in Sec. IIIB. The correspondence between the two formulations is highly suggestive that the phase transition is in fact associated with the collapse of the 2D gas as one approaches the instability from above $`T_c`$. We also performed analgous RPA calculations for the 3D Bose gas with attractive interactions, and found that for the same $`N`$ and $`T`$, the 3D system is always in a stable, uncollapsed phase relative to the 2D gas. We will present a more detailed comparison of the 2D and 3D RPA calculations elsewhere.

## V Acknowledgments

One of us (BVZ) would like to thank Dr. J. P. Carbotte for financial support as well as stimulating discussions. RKB would like to thank Dr. Niall Whelan for relation (22). This research was supported by grants from the Natural Sciences and Engineering Research Council of Canada.


text:
# Time-delayed Spatial Patterns in a Two-dimensional Array of Coupled Oscillators

## Abstract

We investigated the effect of time delays on phase configurations in a set of two-dimensional coupled phase oscillators. Each oscillator is allowed to interact with its neighbors located within a finite radius, which serves as a control parameter in this study. It is found that distance-dependent time-delays induce various patterns including traveling rolls, square-like and rhombus-like patterns, spirals, and targets. We analyzed the stability boundaries of the emerging patterns and briefly pointed out the possible empirical implications of such time-delayed patterns.
PACS numbers: 47.54.+r, 05.45.Ra, 82.40.Ck, 89.75.Kd

Spatiotemporal patterns arise in numerous physical, chemical, and biological systems . The brain, one of the most complex systems, is now also known to generate spatiotemporal patterns such as plane waves and spirals . Early studies, in this context, have concentrated on exploring the dynamics arising in a set of coupled oscillators. It has been, however, generally assumed that interactions among the individual oscillators are instantaneous. Recently, noting that an inclusion of time delay is more natural in realistic systems, several authors have investigated the effects of time delay, and have found the consequent dynamic phenomena such as multistability, desynchronization, clustering, amplitude death, anticipated synchronization and slow switching . It has been also found recently that a time delay, proportional to a distance between elements, could produce a propagating structure in a one-dimensional coupled-oscillators.

In this paper, we discuss further the effect of time delays on the dynamics in a complex system in two dimensions, employing a two-dimensional coupled oscillator model. It is found that time delays among constitutive elements alone can induce definite spatiotemporal patterns in two-dimensional complex system.

We start with the following coupled oscillator model equations for phase variables $`\theta _{ij}(t)`$ ($`i`$ and $`j`$ are integers):

$$\dot{\theta }_{ij}(t)=w+\frac{K}{N(r_0)}\underset{k,l}{\overset{0<r_{kl,ij}r_0}{}}W(r_{kl,ij})\mathrm{sin}\left[\theta _{kl}(t\frac{r_{kl,ij}}{v})\theta _{ij}(t)\right],$$
(1)
with $`N(r_0)=\underset{k,l}{\overset{0<r_{kl,ij}r_0}{}}W(r_{kl,ij})`$. Here $`r_{kl,ij}`$ denotes the distance between element $`(i,j)`$ and element $`(k,l)`$. The metric properties of the ensemble are fixed by setting $`r_{kl,ij}=\sqrt{(ik)^2+(jl)^2}`$. In this model an interaction is non-local and time delayed, characterized by coupling length $`r_0`$ and signal propagation speed $`v`$, and time delays are expressed by $`r_{kl,ij}/v`$. Assuming that an individual element interacts with its neighbors located within a finite radius, we introduce a coupling radius denoted by $`r_0`$; We also discussed the case of $`r_0\mathrm{}`$ in this study. Further we introduce a weighting function $`W(r_{kl,ij})`$ which mimics the interaction strength depending on a distance between interacting elements. Since there is no law about $`W(r_{kl,ij})`$, as far as we know, we choose $`W(r_{kl,ij})=1/r_{kl,ij}`$ assuming that the Gauss’ law of interactions in two dimensional cases is also valid here.

The above equation without time delays is a model of coupled limit cycle oscillators in the so-called phase approximation . It is known that, in the absence of delays $`r_{kl,ij}/v`$ , for any positive values of the coupling coefficients $`K`$, all the oscillators finally form a planar solution where all the oscillators are synchronized, oscillating with frequency equal to that of the individual oscillators . The present study of the model (1), on the other hand, shows that the final steady states correspond to planar solutions with frequencies different from $`w`$ or to frequency-synchronized nontrivial phase configurations.

To find out possible phase configurations, let us assume steady solutions in which all the oscillators have the same frequency $`\mathrm{\Omega }`$. Then the solutions of the model can be written as $`\theta _{ij}(t)=\mathrm{\Omega }t+\varphi _{ij}`$. Substitution of these solutions into Eq. (1) yields

$$\mathrm{\Omega }=w+\frac{K}{N(r_0)}\underset{k,l}{\overset{0<r_{kl,ij}r_0}{}}\frac{1}{r_{kl,ij}}\mathrm{sin}(\varphi _{kl}\varphi _{ij}\mathrm{\Omega }r_{kl,ij}/v).$$
(2)
This self-consistency condition requires that the sums $`\frac{1}{N(r_0)}\underset{k,l}{\overset{0<r_{kl,ij}r_0}{}}\frac{1}{r_{kl,ij}}\mathrm{sin}(\varphi _{kl}\varphi _{ij}\mathrm{\Omega }r_{kl,ij}/v)`$ be same for all $`(i,j)`$ . Let us consider the situation in the continuum limit. Differentiating the term with respect to $`\stackrel{}{r}[=(i,j)]`$, we get the following.

$$\frac{\varphi (\stackrel{}{r^{}})}{\stackrel{}{r^{}}}=\frac{\varphi (\stackrel{}{r})}{\stackrel{}{r}}\text{ , for all }\stackrel{}{r}\text{ and }\stackrel{}{r^{}}.$$
(3)
This implies that, in the continuum limit, $`\varphi (\stackrel{}{r})=\gamma \stackrel{}{r}\stackrel{}{e}+\varphi _0`$ is a unique solution, up to the arbitrariness of $`\varphi _0`$ and $`\stackrel{}{e}`$, as long as $`\stackrel{}{e}`$ is a unit vector embedded in the lattice plane. The case of $`\gamma =0`$ corresponds to synchronized planar solutions, while nonzero $`\gamma `$’s make roll solutions traveling in the direction of $`\stackrel{}{e}`$ with the phase velocity $`\mathrm{\Omega }/\gamma `$.

However, the space discreteness of the system admits topological defects, so other phase configurations are also possible . In this study, we set $`w=\pi /5`$, and solved the equation on a 64 by 64 square plane ($`1i64`$, $`1j64`$) with periodic boundary conditions. We checked that different size of lattice does not alter the result. We tried several kinds of initial conditions. As the initial conditions in phase of the form $`\theta _{ij}(0)=\theta _0+\eta _{ij}`$ in this study, where $`\eta _{ij}`$ denotes a random number, we considered two different cases: (i) $`0\eta _{ij}<2\pi `$, and (ii) $`0\eta _{ij}<10^6`$. In addition, for the dynamic history of individual oscillators for $`t<0`$, we considered two situations (i) where the oscillators were stationary, and (ii) where the oscillators have evolved independently from one another with their natural frequency $`w`$.

Figure 1, which is made with the parameters $`K=0.6`$ and $`r_0=10`$, shows effects of the delay $`1/v`$. When the delay $`1/v`$ is small ($`<0.6`$), only planar solutions($`\gamma =0`$) occur. When the delay is large enough($`>1.1`$), planar solution is not available but patterns with nonzero $`\gamma `$’s occur. In the intermediate region ($`0.61/v1.1`$), planar solutions as well as patterns are possible. The lines in Fig.1 are the numerically obtained solutions of Eq. (2) with $`\varphi _{ij}=\gamma i`$, $`\gamma =0,\frac{5\pi }{32},\frac{6\pi }{32}`$, and $`\frac{7\pi }{32}`$. We see that the wavelengths of patterns become shorter generally as delay increases.

Figure 2 provides the phase diagrams of the model. If the coupling strength $`K`$ is above $`K_{c2}`$ (denoted by open squares), only a synchronized planar solution appears, but if $`K`$ is below the critical value $`K_{c1}`$ (denoted by filled squares), planar solution is not available and the system evolves into patterns such as rolls \[as in Fig. 3(a)\], or other phase distributions \[as in Figs. 3(b) - 3(d)\] depending on initial configurations. In the regions of $`K_{c1}<K<K_{c2}`$, we see both patterns and planar solutions depending on initial configurations. Figure 4 shows the patterns that appear only in the limited regions in the phase diagram. When $`r_0`$ is relatively small ($`5`$), in other words, when the effect of space discreteness becomes important, spirals occur more frequently \[See Fig. 4(a)\] than other configurations. Furthermore, the model sometimes organizes target patterns from slightly perturbed homogeneous phase configurations for parameter $`KK_{c1}`$ \[see Fig. 4(b)\]. With $`KK_{c2}`$, there emerge patterns embedded in a planar solution as shown in Fig. 4(c).

Note that the arrangement of the phase singularities plays an essential role in forming the patterns of Figs. 3(b) - (d). The singularities are possible due to the discreteness of the space as in spiral patterns, but are dynamically arranged in a rhombus \[Fig. 3(b)\] or a square lattice \[Fig. 3(c)\]; the phases of oscillators make rotating waves around those singularities . However, rolls are more prevalent than squares or rhombus patterns in the entire region of the diagram.

There seem to exist three kinds of multistabilities in our model. First, as stated above, we can see planar solutions, rolls or other phase distributions for the same parameters depending on initial configurations \[see Figs. 1-3\]. Second, there are multiple states with different wavelengths and frequencies even for the same parameters in some cases, which is why we see multiple $`\mathrm{\Omega }`$’s of patterns in Fig. 1. This kind of multistability is also reported in one-dimensional lattices . Third, different initial configurations make the different angles between two primitive translation vectors of patterns like Figs. 3(b) and 3(c). It seems that multistabilities are the major characteristics of time delayed systems .

In the systems of two-coupled oscillators , or globally coupled oscillators with uniform time delays , the dynamics are determined by the value of $`\mathrm{\Omega }\tau `$, where $`\mathrm{\Omega }`$ is the frequency of the synchronized state and $`\tau `$ is a time delay between oscillators. The physical meaning of the value is clear. It denotes the virtual phase difference between oscillators due to delay. If it is near $`2n\pi `$, where $`n`$ is an integer, each oscillator can keep pace with others to make the whole system synchronized, otherwise it cannot . Likewise, we suggest that the dynamics of our model is determined by the quantity $`\mathrm{\Omega }r_0/v`$, which we call $`\mathrm{\Theta }`$ here, where $`\mathrm{\Omega }`$ is a frequency of a synchronized planar solution. Our conjecture is that there exists some value $`\mathrm{\Theta }_c`$ so that planar solutions are possible only for $`\mathrm{\Theta }<\mathrm{\Theta }_c`$ .

We measured $`\mathrm{\Theta }_c`$ directly from numerical simulations using Eqs. (1) and (2) in Fig. 5. We see that $`\mathrm{\Theta }_c`$ is a nearly constant value ($`2.26`$), from which we can calculate $`K_{c1}`$ in Fig. 2. Substituting the roll solution $`\varphi =\gamma x`$ to the continuum limit of Eq.(2) makes the following equation.

$$\mathrm{\Omega }=w\frac{K}{r_0}_0^{r_0}\mathrm{sin}(r\mathrm{\Omega }/v)J_0(\gamma r)𝑑r,$$
(4)
where $`J_0`$ is a Bessel function. In the case of planar solutions ($`\gamma =0`$), the above equation can be integrated to give

$$K=\frac{\mathrm{\Theta }w}{2\mathrm{sin}^2\frac{\mathrm{\Theta }}{2}}(1\frac{\mathrm{\Theta }}{wr_0\frac{1}{v}}).$$
(5)
The solid lines in Figs. 2 (a) and (b) are obtained by substituting $`\mathrm{\Theta }=2.26`$ to Eq.(5), and show good agreements with the simulation results.

It is useful to look at what the existence of constant $`\mathrm{\Theta }_c`$ implies. It means the following. Coupled oscillators tend to be synchronized, but when too large delays make the $`\mathrm{\Theta }`$ value go beyond the tolerance limit $`\mathrm{\Theta }_c`$, they fail to become a synchronized planar solution and take other configurations. The occurrence of incoherence in the coupled oscillators with uniform time delay  can now be explained in the same manner.

It is also noteworthy that, in the limit of $`r_0\mathrm{}`$, Eq. (5) converts to $`K_c=\frac{\mathrm{\Theta }_cw}{2\mathrm{sin}^2\frac{\mathrm{\Theta }_c}{2}}`$, which is independent of $`1/v`$. This is surprising, for it means that we must consider the effects of time delay however small when the interaction length is infinite as in most physical systems.

To summarize, we have investigated the effects of time delayed interactions in an ensemble of coupled oscillators in two dimensions, and found that distance-dependent time delay alone can induce various spatial patterns, while without delay or with uniform delay alone only synchronized planar solutions might appear. We have analyzed also that the stability of planar solutions can be determined by the quantity $`\mathrm{\Omega }r_0/v`$. Our study may find its relevance in recent experimental observations that non-locally interacting oscillators, which are especially probable for neurobiological systems, exhibit traveling rolls and spirals . Finally, our results seem to suggest that time delay may play a significant role in the study of memory storage and information processing based on the spatiotemporal activities of neurons .

We thank Hwa-Kyun Park for useful discussions. This work was supported in part by the Interdisciplinary Research Program (Grant No. R01-1999-00019) of the Korea Science and Engineering Foundation, and in part by the Brain Korea 21 Project in 2000. We also appreciate the support of the Korea Research Foundation (Grant No. KRF-2000-015-DP0097).

Figure Captions

Synchronization frequency $`\mathrm{\Omega }`$ as a function of a delay ($`K=0.6,r_0=10`$); squares indicate planar solutions; circles indicate patterns. Lines are the solutions of Eq.(2). See the text for details.

The phase diagrams of the model with $`r_0=10`$ in (a), and with $`1/v=1`$ in (b). Above open squares, there is no pattern available, and only planar solutions appear. Between open and filled squares, the system evolves into planar solutions or patterns according to initial conditions. Below the filled squares, planar solution is not possible here, and only patterns appear. Solid lines are from Eq.(5). See the text for details.

Typical patterns generated in the model : (a) a traveling roll, (b) a rhombus-like pattern, (c) a square-like pattern, (d) coexistence of a rhombus-like pattern and a roll. $`K=1.25`$, $`r_0=25`$, and $`1/v=1`$, for (a), (b), and (c). Note that planar solutions are also found for this set of parameter values.; $`K=0.6`$, $`r_0=10`$, and $`1/v=1.1`$, for (d). See Fig. 4 (d) for the color coding.

Patterns obtained in simulations of the model : (a) spirals ($`K=1.0`$, $`r_0=4`$, and $`1/v=1`$), (b) a target pattern ($`K=0.4`$, $`r_0=7`$, and $`1/v=1`$), (c) a spiral embedded in a planar solution ($`K=0.8`$, $`r_0=8`$, and $`1/v=1`$). In (c), the planar solution oscillates much more slowly than the spiral. Figure (d) shows the color code used in this paper.

The quantity $`\mathrm{\Omega }r_0/v`$ determines the stability of planar solutions. We cannot see planar solutions above the filled squares.


text:
# Comparative Studies of Line and Continuum Positron Annihilation Radiation

## I Introduction

Nine years of observations made with the Oriented Scintillation Spectrometer Experiment (OSSE) on-board NASA’s COMPTON observatory (1991-2000)john93 , eight years of observations made with the Gamma-Ray Spectrometer on-board the Solar Maximum Mission (SMM) (1980-1989)shar88 , and two years of observations made with the Transient Gamma-Ray Spectrometer (TGRS) on-board the WIND mission (1995-1997) harr98  have been utilized to study the galactic distribution of positron annihilation radiation. The OSSE instrument featured a 3.8 x 11.4 FWHM FoV, a $``$3 x 10-5 photons cm-2 s-1 line sensitivity (106 s on-source time), and a 45 keV energy resolution at 511 keV. These detector attributes have permitted the first detailed studies of the distribution of annihilation radiation in the inner radian of the Galaxy. The annihilation of positrons with electrons gives rise to two spectral features, a line emission at 511 keV and a positronium continuum emission (which increases in intensity with energy roughly as a power law up to 511 keV and falls abruptly to zero above 511 keV)ore49 . The TGRS instrument, which featured a germanium detectors with excellent energy resolution, has demonstrated that the integrated flux from the inner radian is best described as a narrow 511 keV line (FWHM $``$ 1.8 keV) and a positronium continuum to 511 keV line ratio of $``$ 3.6 (which corresponds to a positronium fraction of fPs=0.94)harr98 .

Purcell et al. (1997) (hereafter PURC97)purc97  reported results from OSSE/SMM/TGRS studies of the 511 keV line component of annihilation radiation. They found the 511 keV emission to be comprised of three components; 1) an intense bulge emission, 2) a fainter disk emission, and 3) an enhancement of emission at positive latitudes (hereafter called a PLE). The PLE was also reported by Cheng et al. (1997)chen97 , and has been interpreted to be an “annihilation fountain” by Dermer & Skiboderm97 . PURC97 characterized the emission via mapping, employing the SVD matrix inversion algorithm, and via model fitting, testing the combination of a spheroidal Gaussian bulge, a disk that is flat in longitude to $`\pm `$40 and Gaussian in latitude (FWHM = 9), and a spheroidal PLE. The two characterizations differ in the thickness of the Gaussian disk (SVD being narrower) and the extension of the PLE. The enhancement of the PLE differed between the two characterizations, varying from 1.5 x 10-4 photons cm-2 s-1 for the SVD map to 9 x 10-4 photons cm-2 s-1 for the broad 2D Gaussian PLE (FWHM = 16.4). A parallel study of both line and continuum annihilation radiation along the galactic plane by Kinzer et al. (1996,2001)kinz96 ; kinz01  reported that positronium continuum emission is similarly distributed as 511 keV line emission. The Kinzer studies did not investigate the PLE.

We report here updates from our continuing analysis which extends the study of PURC97 (see also Milne et al. (1998,1999) miln98 ; miln99 ). The primary differences between current studies and PURC97 are; 1) the inclusion of more observations, both archival and data collected after PURC97, and 2) reporting maps of the positronium continuum emission in addition to the 511 keV line. To extract the positronium continuum component from the underlying galactic continuum emission, we widened the spectral modeling to include thermal bremsstrahlung and exponentially-truncated power-law models. We also removed high-energy diffuse continuum emission following a prescription from Kinzer et al. (1999), distributing the emission spectrally according to an $`\alpha `$ = -1.65 power-law and spatially according to a 90 x 5 2D Gaussiankinz99 . Two maps of the 511 keV line emission and two maps of the positronium continuum emission are shown in Figure 1.

## II Discussion

Although not identical, the two 511 keV maps share certain fundamental features. Both exhibit an intense bulge emission and a fainter planar emission. The regions that appear anomalous in the RL map (relative to symmetrical bulge and disk emissions) are also anomalous in the SVD map. The positronium continuum maps are also dominated by an intense bulge and a fainter disk component. Pairings of bulge and disk components suggest the same families of solutions. Both suggest that the bulge-to-disk ratio can vary from 0.2 -3 depending upon whether the bulge component features a halo (which leads to a large B/D). The most noticeable difference between the four maps is how the PLE is characterized.

Both mapping algorithms suggest enhancement of 511 keV line emission from the region reported in PURC97, although at lower flux levels ($``$ 1.0 x 10-4 photons cm-2 s-1). The broad 2D Gaussian PLE described in PURC97 is not an acceptable solution; the quality of fit is much worse than solutions without a PLE. By contrast, there is no suggestion of a PLE at the location reported in PURC97 ($`l,b`$ = -2, 7) in either of the positronium continuum maps. Shown in Figure 2 are 1 wide cuts cuts through the four maps, taken at 70 relative to the negative-longitude galactic plane. It is apparent that the four curves are consistent with the same distribution towards negative latitudes, but differ at positive latitudes. The lower panel of Figure 2 shows the positive latitude cut after subtraction of the mirror negative latitude cut. The enhancement at 511 keV has a corresponding deficit for positronium continuum for both mapping algorithms.

An enhancement of 511 keV emission that is due to an additional source of positrons (relative to the Galaxy-wide bulge and disk components) would be expected to feature an enhancement of both annihilation components. For the lowest fPs value of zero, we would expect the positronium continuum emission to be symmetric and the 511 keV line emission to be enhanced. A positronium continuum deficit at positive latitudes would not be expected. If, alternatively, the PLE is a region with no excess of positrons but is where the local fPs varies from the integrated fPs=0.94 (suggested both by OSSE analysiskinz01 , and by wide FoV germanium detectorsharr98 ), then the 511 keV line emission could be enhanced and the positronium continuum emission could be deficient. However, variations of the fPs it do not conserve photon flux. The 1:1 enhancement-to-deficit ratio suggested in Figure 2b would not be expected. A third possibility is that the enhancement is instead due to the influence of gamma-ray sources that corrupt the spectral fitting. Two possible mechanisms for this biasing could be from observations made while a source is exhibiting a hard X-ray flare, and/or if the source exhibits a previously undetected hard tail. The fact that both the line and positronium continuum components peak at 511 keV combined with the 45 keV FWHM energy resolution of the OSSE instrument mean that as many as 30% of the counts at that energy are due to positronium continuum photons. It is unclear whether flaring gamma-ray sources can bias the the spectral fitting to a large enough extent to entirely account for the apparent asymmetry in the annihilation radiation.

A few compact sources which produce this type of biasing have been identified. At the present time, it has not been established which of the three explanations is correct. The CGRO/BATSE instrument made observations of these sources which were nearly simultaneous with the OSSE observationsharm92 . It is an objective of the current analysis effort to determine whether joint analysis of OSSE & BATSE observations of these sources will permit the unambiguous extraction of annihilation radiation from this complex environment. Fortunately, although other regions of the inner galactic radian may be similarly biased, the majority of the region is not expected to have been affected. The wealth of information available both from the complete OSSE data-set as well as from the expanding data-set of monitoring of compact sources with the BATSE instrument may permit definitive statements as to the galactic distribution of positron annihilation radiation, particularly of the existence of a PLE.


text:
# Galaxy clustering in the Sloan Digital Sky Survey (SDSS): A first comparison with the APM Galaxy Survey

## 1 Introduction

The Sloan Digital Sky Survey (SDSS) uses a dedicated 2.5 m telescope and a large format CCD camera to obtain images of over 10,000 $`deg^2`$ of high Galactic latitude sky in five broad bands. The SDSS represents a new paradigm of scientific project, where digital data is made available on-line to the community via ”virtual observatories”. On June 2001 the SDSS collaboration made an early data release (EDR) publicly available. In this paper, we take this opportunity to do a study of some large scale structure aspects of the EDR. We compare this new data with the APM Galaxy Survey (Maddox etal 1990), which is based on 185 UK IIIA-J Schmidt photographic plates each corresponding to $`6\times 6`$ square degrees on the sky limited to $`b_J20.5`$ and covering $`b<40`$ and $`\delta <20`$ degrees. These fields where scanned by the APM machine and carefully matched using the plate overlaps.

The APM Survey has produced one of the best estimates of the angular galaxy 2-point correlation function $`w_2(\theta )`$ to date. Its shape on large scales led to the discovery of “extra” large-scale power, and gave early indications for the paradigm shift out of the standard (Einstein-de Sitter) CDM model (Maddox etal 1990). Frieman & Gaztañaga (1999, FG99 hereafter) estimated the 3-point galaxy correlation function in the APM Galaxy Survey and its comparison with theoretical expectations (see also Peebles 1980, Fry 1984, Juszkiewicz, Bouchet & Colombi 1993, Fry & Gaztañaga 1993, Gaztañaga 1994, Bernardeau 1994, Fosalba & Gaztañaga 1998, Buchalter, Jaffe & Kamionkowski 2000, Scoccimarro etal. 2001, and references therein). For the first time, the APM measurements extended to sufficiently large scales to probe the weakly non-linear regime with a reliable Survey. The results are in good agreement with gravitational growth for a model with initial Gaussian fluctuations. They also indicate that the APM galaxies are relatively unbiased tracers of the mass on large scales and provide stringent constraints upon models with non-Gaussian initial conditions.

Linear CCD devices provide more reliable magnitudes than the traditional non-linear system of visual magnitudes. Scanned plates should do better, but there has been extended discussions regarding variable sensitivity inside individual Schmidt plates and large-scale gradients in the APM survey calibration. A CCD survey, such as the SDSS, reduces these systematic effects and could provide more accurate photometry and luminosity function estimators (eg Blanton etal 2001, see also Gaztañaga & Dalton 2000). Thus, a natural question arrises: is the large scale clustering trace by the new SDSS digital data compatible with the one measured in the old photographic APM scans? As there is no spatial overlap between these two data sets, we will compare the statistical properties, such as correlation functions and moments of counts in cells. The answer to this question will bring new evidence to our understanding of structure formation in the lines mentioned above.

## 2 SDSS sample and pixel maps

We download data from the SDSS public archives using the SDSS Science Archive Query Tool (sdssQT, http://archive.stsci.edu/sdss/software/). We select objects from an equatorial SGC (South Galactic Cap) strip 2.5 wide ($`1.25<DEC<1.25`$ degrees.) and 66 degrees long ($`351<RA<56`$ degrees.) which seems to be homogeneous. This strip (SDSS numbers 82N/82S) corresponds to some of the first runs of the early commissioning data (runs 94 and 125) with variable seeing conditions (a few tenths of arc-second within scales of a few degrees. 111See http://www-sdss.fnal.gov:8000/skent/seeingStatus.html). To avoid confusion with the EDR equatorial strip in the North (studied by the SDSS collaboration, eg in Scranton et al 2001) we will refer to this strip as EDR/SGC.

We first select all galaxies brighter than $`u^{}=22.3,g^{}=23.3,r^{}=23.1,i^{}=22.3,z^{}=20.8`$, which corresponds to the SDSS limiting magnitudes for 5 sigma detection in point sources (York etal. 2000). Galaxies are found from either the $`1\times 1`$, $`2\times 2`$ or $`4\times 4`$ binned CCD pixels and they are deblended by the SDSS pipeline (Lupton etal 2001). There are about 375000 objects classified as galaxies in the EDR/SGC. Figure 1 shows the number counts (surface density) for these 375000 galaxies as a function of the magnitude in each band, measured by the SDSS modified Petrosian magnitudes $`m_u^{}`$, $`m_g^{}`$, $`m_r^{}`$, $`m_i^{}`$ and $`m_z^{}`$ (see Yasuda et al 2001 for a discussion of the SDSS counts). Continuous diagonal lines show the $`10^{0.6m}`$ expected for a low redshift homogeneous distribution with no k-correction, no evolution and no-extinction. The limiting magnitudes for detecting galaxies in this strip seems to be closer to $`m_u^{}<21.5`$, $`m_g^{}<21`$, $`m_r^{}<20.5`$, $`m_i^{}<20`$ and $`m_z^{}<19.5`$.

We next select galaxies with SDSS modified Petrosian magnitudes $`16.8<m_g^{}<19.8`$. We choose $`g^{}`$ because this is the closest SDSS band to the APM blue $`b_J`$. The range $`16.8<m_g^{}<19.8`$ results from requiring that the SDSS surface density in a 3 magnitude slice equals that of the APM in the $`17<b_J<20`$ slice (ie $`300`$ gal per $`deg^2`$). We have checked that small departures from this range produce similar clustering results. There are about 54000 galaxies that match this criteria in the EDR/SGC. Figure 1 shows that the number counts in this range of magnitudes (limit by the two dashed vertical lines in the figure) follow well the $`10^{0.6m}`$ relation (continuous line), indicating that the sample selected is homogeneous and completed to this magnitude. We do not apply any mask to the EDR/SGC sample or try to make any correction for variable seeing or variable extinction across the strip. Visually, there are no obvious holes or inhomogeneities in the pixel maps at this depth. We will assume that variations in extinction and seeing do not affect clustering of the relatively bright galaxies we are considering (see a detailed study by Scranton etal 2001).

Finally, we produce equal area projection pixel maps of various resolutions 222Available on www.inaoep.mx/ $`\stackrel{~}{}`$ gazta/Mapg18.9.cut.pgm.gz, similar to those made for the APM (see Plate 1 in Maddox etal 1990). Figure 2 shows the EDR/SGC pixel maps with a resolution of about $`7`$ arc-min. It is apparent that there are important surface density variations along the strip, even on the largest scales. One can also see filamentary structure and a overall texture similar to that in the APM pixel maps (see Plate 1 in Maddox etal. 1990) or the Lick maps (see cover of Peebles 1993).

## 3 Clustering comparison

Our goal is to compare the APM to the new SDSS data. Out of the APM Survey we considered a $`17<b_J<20`$ magnitude slice in an equal-area projection pixel map with a resolution of $`3.5`$ arc-min, that covers over $`4300`$ $`deg^2`$ around the SGC. The APM can fit over 25 strips of similar area, shape and depth as the SDSS strip shown in Figure 1. To study sampling and biasing effects on the SDSS clustering estimators (due to the limited survey area), we have cut 10 strips out of the APM map. Each of these 10 sub-samples are identical in size and shape to the EDR/SGC SDSS strip, and are separated from each other by few degrees, so that they can be considered as independent. This will allow us to set $`90\%`$ confidence intervals against the null hypothesis that the EDR/SGC strip is compatible with the APM map. In all cases we correct the clustering in the APM maps for a $`5\%`$ contamination of randomly merged stars (see Maddox etal 1990), ie we scale fluctuations up by $`5\%`$ (see also Gaztañaga 1994).

### 3.1 Smoothed 1-point Moments

We first compare the lower order moments of counts in cells of variable size $`\theta `$ (larger than the pixel map resolution). We follow closely the corresponding APM analysis in Gaztañaga (1994, see also Szapudi etal 1995) and use the same software and estimators here for the SDSS. The top panel in Fig. 3 shows the variance of fluctuations in density counts $`\delta \rho /\overline{\rho }1`$ smoothed over a scale $`\theta `$: $`\overline{w_2}<\delta ^2(\theta )>`$, which is plotted as a function of the smoothing radius $`\theta `$. The shaded region shows the $`90\%`$ confidence interval for the SDSS results to be compatible to the APM map. This shaded region brackets the minimum and maximum results at each scale in the 10 APM sub-samples that mimic the EDR/SGC. The individual results in each subsample are strongly correlated so that the whole curve for each subsample scales up an down inside the shaded region, ie there is a strong covariance at all separations due to large scale density fluctuations. Triangles correspond to the mean of the 10 APM sub-samples, which are slightly bias down as compared to the whole APM (continuous line). The SDSS results (open squares) are all within the $`90\%`$ confidence region set by the APM. Because of the strong covariance mentioned above, we do not expect the SDSS values to scatter around the APM values, but rather we expect the whole curved to be shifted, as shown in the Figure. This is also the case for the other statistical measures presented below.

The extend of the shaded area would narrow if we increase the total solid angle of the subsamples. Numerical simulations show that the corresponding error-bar is about a factor of 4-6 smaller for a survey as large as the APM.

The bottom panel in Fig. 3 shows the corresponding comparison for the normalized angular skewness:

$$s_3(\theta )\frac{<\delta ^3(\theta )>}{<\delta ^2(\theta )>^2}\frac{\overline{w_2}(\theta )}{\overline{w_3}(\theta )}$$
(1)
Again the mean of the individual APM zones (triangles) does not equal the whole APM estimation (continuous line). This is due to estimation (or ratio) bias (see Hui & Gaztañaga 1999). The SDSS EDR/SGC data on $`s_3`$ shows an excellent agreement with the APM at the smaller scales (in contrast with the EDSGC results, see Szapudi & Gaztañaga 1998). On larger scales the SDSS values are smaller, but the discrepancy is not significant within the $`90\%`$ region and given the strong covariance of individual subsamples. Similar results are found for higher order moments. As we approach the scale of $`2`$ degrees, the width of our strip, it becomes impossible to do counts for larger cells. It is therefore interesting to study the correlation functions, which should be less affected by boundary effects (although sampling effects will still be important).

### 3.2 N-point Correlation functions

We next study the angular 2-point, $`w_2<\delta _1\delta _2>`$, and 3-point, $`w_3<\delta _1\delta _2\delta _3>`$, correlation functions. Here we follow closely the notation, estimators and software used in FG99. The 3-point function is normalized as:

$$q_3\frac{w_3(\theta _{12},\theta _{13},\theta _{23})}{w_2(\theta _{12})w_2(\theta _{13})+w_2(\theta _{12})w_2(\theta _{23})+w_2(\theta _{13})w_2(\theta _{23})}$$
(2)
where $`\theta _{12}`$, $`\theta _{13}`$ and $`\theta _{23}`$ correspond to the sides of the triangle form by the 3 angular positions of $`\delta _1\delta _2\delta _3`$. Here we will consider isosceles triangles, ie $`\theta _{12}=\theta _{13}`$, so that $`q_3=q_3(\alpha )`$ is given as a function of the interior angle $`\alpha `$ which determines the other side of the triangle $`\theta _{23}`$ (FG99). We also consider the particular case of the collapsed configuration $`\theta _{23}=0`$, which corresponds to $`<\delta _1\delta _2^2>`$ and is normalized in slightly different way (see also Szapudi & Szalay 1999):

$$c_{12}\frac{<\delta _1\delta _2^2>}{<\delta _1\delta _2><\delta _1^2>}2q_3(\alpha =0).$$
(3)
Top panel in Fig.4 shows the 2-point correlation as a function of the separation $`\theta `$ between points (here small pixel cells). Results from both surveys agree remarkably well up to 1 degree. In a similar way as to what happened with the variance (ie Fig.3) the SDSS EDR/SGC sample shows a flatting of the $`w_2`$ slope at large scales. At scales bigger than $`3`$ degrees this discrepancy becomes significant and it is incompatible with the APM subsamples with a confidence level higher that $`90\%`$. This deviation could be due to systematic gradients in the SDSS raw data that we are using (eg due to large seeing variation across the strip, see §2 above). But it could also be due to an underestimation of the sampling errors in the APM (eg due to smoothing of the APM calibration on the scale of the subsamples). We need to consider larger and better calibrated areas to confirm or refute with the SDSS the features in APM galaxy power spectrum shape discussed in Gaztañaga & Baugh (1998, and references therein) and Gaztañaga & Juszkiewicz (2001). At smaller scales, the agreement is excellent and the measured shape of the SDSS 2-point function confirms the idea that galaxies are ’anti-bias’ with respect to the $`\mathrm{\Lambda }`$CDM model (eg Gaztañaga 1995, Jenkins etal 1998, and references therein).

Bottom panel of Fig.4 and Fig.5 shows the reduced 3-point function $`q_3`$ for different configurations. Fig.5 compares $`q_3(\alpha )`$ for isosceles triangles of side $`\theta _{12}=\theta _{13}=0.5`$ degrees. Bottom panel of Fig.4 shows the collapsed case $`c122q_3(\alpha =0)`$ as a function of the other triangle side $`theta=\theta _{12}=\theta _{13}`$. The reduced 3-point function in Fig.5 is close to zero and the signal seems stronger on smaller scales. In all cases, the SDSS results are within the $`90\%`$ confidence of the APM values. Note again the bias between the mean of the 10 APM strips (triangles) and the whole APM (continuous line), and also the strong covariance across the curves. As mentioned above both of these results are expected and due to (very) large scale density fluctuations.

## 4 Discussion and conclusions

We have presented a rigorous $`90\%`$ confidence test for the hypothesis that the clustering in the SDSS EDR/SGC strip is compatible with the one in the APM data. The test is passed for clustering at scales smaller than 2 degrees (which corresponds to $`15Mpc/h`$), and seems to fail (at $`90\%`$ confidence) at larger scales. This could be due to large scale variations in the seeing (on scales of a few degrees) present in the early SDSS commissioning data (EDR/SGC) that we are using (see §2), but they could also reflect the fact that the APM sampling errors are not accurate at this level. In principle, it is possible to correct for these seeing gradients by using the appropriate mask (see Scranton etal 2001), but we have not attempted to do this here.

The agreement between these two catalogues, on scales smaller than 2 degrees, is not trivial: correlation functions are quite sensitive to systematic effects and clustering for galaxies of fainter magnitudes or different colors is quite different. Moreover, we have not rescale or match the amplitudes of the two catalogues: the agreement comes naturally when we just select galaxies with the same colors, magnitude range and similar surface density in both samples. It is hard to belive that this is just a coincidence. These results rather illustrates that EDR data and software pipelines from SDSS work well and are suited to do analysis of large scale clustering.

We find it remarkable that two totally independent surveys, separated in space (ie angular position and spectral filtering) and time (over 10 years in technological and software development) provide very similar results for the large scale angular clusterings trace by galaxies. This is true both visually and by accurate statistical comparisons. It illustrates that large scale clustering analysis is becoming ”repeatable” and therefore that our conclusions could stand on solid scientific grounds. In particular, one conclusion that we have been able to confirm in this work, is that the higher order correlations in galaxy samples indicate that gravitational growth from Gaussian initial conditions is responsible for the hierarchical structures we see in the sky (eg FG99, Scoccimarro etal 2001, and references therein). When more of the SDSS data is analyzed we can hope to make much higher precession testing of our models for structure formation (eg Tegmark etal 1998; Colombi etal 2000).

## Acknowledgments

I acknowledge support by grants from IEEC/CSIC and DGI/MCT BFM2000-0810. Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/.


text:
# Phase transitions in the mesoscopic superconducting square

\[

## Abstract

We solve the Ginzburg-Landau equation (GLE) for the mesoscopic thin film of the square shape in the magnetic field. In the limit of Ginzburg-Landau parameter $`\kappa \mathrm{}`$ we find a series of first and second order phase transitions as temperature and/or magnetic field changes. First order phase transitions between giant flux states can be described with a simple variational procedure. We discuss the similarity with rotating liquid $`He^4`$ and derive a simple formula for $`H_{c1}`$. We identify order parameters based on symmetry arguments and we propose Landau functional describing the second order phase transition.

\]

Advances in nanotechnology and constantly shrinking semiconductor devices have motivated researches to study properties of mesoscopic superconducting samples. One line of research in this field has focused on the problem of the phase transitions in the mesoscopic superconducting sample under the influence of the external magnetic field . There are two characteristic limits in which phase transitions have different properties. If the size of the sample $`a\xi `$, where $`\xi `$ is the superconducting coherence length, and applied field are large enough, there are many vortices in the sample. In this case long-range interaction between vortices and image vortices is screened by spontaneous creation of vortex loops near the sample boundary. It leads to the decrease of the surface barrier for the vortex to penetrate into the sample. In the opposite case, when $`a\xi \lambda `$, with $`\lambda `$ being the London penetration depth, there are only few vortices in the sample. The standard Abrikosov approach  must be modified because of the strong influence of the sample boundaries. In this case magnetic field is almost uniform throughout the sample and long-range forces between vortices are not important. Thermodynamics of this system is determined by the short-range repulsion of vortices and Bean-Livingston barrier forces.

Different approaches have been applied for the investigation of phase transitions in the latter limit. Most of them consider disk geometry. Buzdin and Brison applied electrostatic formalism to consider influence of the barrier on the vortex structure of the thin superconducting disk . Within this approach vortices are replaced by the hard-core particles interacting through Coulomb forces and the giant vortex has never been discussed. Numerical solution of GLE for the same geometry reveals a series of the first and second order phase transitions in the superconducting disk. Such transitions take place between giant vortex states with different vorticity as well as between a giant vortex state and a multi-vortex state as the external field changes . We emphasize an important difference between the disk and the square geometry. Solution of the linearized GLE, describing the nucleation of superconducting order parameter near $`H_{c2}`$ line for the disk, always corresponds to the giant vortex state. On the other hand, as it was demonstrated by Chibotaru et al. , there are many well separated zeros of the order parameter in the case of the square sample. Consequently, the behavior of the square sample near $`H_{c2}`$ line should be qualitatively different from the disks. On the basis of the solution of the linearized GLE the appearance of the antivortex in the center of the sample has been predicted.

In the this paper we investigate phase transitions in the superconducting film of a square shape as a function of temperature $`T`$ and external magnetic field $`H`$ in the limit $`a\xi \lambda `$. We solve GLE for the thin superconducting film with the thickness $`d\xi `$. We also assume that superconducting film is an extreme type two superconductor, where the Ginzburg-Landau parameter $`\kappa =\lambda /\xi 1`$. We show that configuration with one anti-vortex in the center and four vortices on the diagonals of the square is unstable when we move away from the $`H_{c2}`$ line and nonlinear term in the GLE is considered. On the contrary, at higher magnetic field, the configuration with four vortices on diagonals of the square remains stable. We find a sequence of phase transitions of the first order between giant vortex states as well as between multi-vortex states with different vorticity. Second order phase transition takes place when a giant vortex state splits into a multi-vortex state with simultaneously breaking the $`C_4`$ symmetry. Such transitions are discussed in terms of the phenomenological theory of Landau.

GLE for dimensionless complex order parameter $`\psi `$ has the following form:

$$\xi ^2(i+\frac{2\pi 𝐀}{\mathrm{\Phi }_0})^2\psi \psi +\psi |\psi |^2=0$$
(1)
here $`\xi =\frac{\mathrm{}^2}{4m|\alpha |}`$, $`\alpha `$ is temperature dependent parameter of Ginzburg-Landau expansion for the free energy, $`\mathrm{\Phi }_0`$ is the flux quantum, $`𝐀`$ is the vector potential $`𝐇=\times 𝐀`$. The second GLE equation for the vector potential can be written as:

$$\times \times 𝐀=i\frac{\mathrm{\Phi }_0}{4\pi \lambda ^2}(\psi ^{}\psi \psi \psi ^{})\frac{|\psi |^2𝐀}{\lambda ^2}.$$
(2)
Since we consider the case of a small mesoscopic square where $`a\xi \lambda `$, the magnetic field is uniform in the film. The correction to the external field is of the order of $`1/\kappa ^2`$ and may be found by solving Eq.(1) while assuming uniform magnetic field and substituting the solution of Eq.(1) to Eq.(2). Such a solution is equivalent to the expansion of the free energy in $`1/\kappa ^2`$ series. In addition to Eq.(1) we have to supply the boundary condition for the superconductor-insulator junction:

$$(i+\frac{2\pi 𝐀}{\mathrm{\Phi }_0})𝐧\psi =0,$$
(3)
where $`𝐧`$ is normal vector to the surface of the sample.

Introducing $`N\times N`$ discrete points in the square we rewrite Eq.(1) in the form of nonlinear discrete Schrödinger equation:

$$\underset{𝐥}{}t_{𝐢+𝐥,𝐢}\psi _{𝐢+𝐥}4t_{𝐢,𝐢}\psi _𝐢\psi _𝐢+\psi _𝐢|\psi _𝐢|^2=0$$
(4)
where summation index $`𝐥=(\pm 1,0)`$, $`(0,\pm 1)`$ points toward nearest neighbours and $`t_{𝐢_\mathrm{𝟏},𝐢}=(\xi N/a)^2exp(\frac{2\pi i}{\mathrm{\Phi }_0}_𝐢^{𝐢_\mathrm{𝟏}}𝐀(𝐫)𝑑𝐫)`$. Equivalent discretization of the boundary conditions, Eq.(3), provides an additional equation which can be directly solved and substituted into Eq.(4). As a result, equation close to the boundary are slightly different from the ’bulk’:

$$\underset{𝐥}{}t_{𝐢+𝐥,𝐢}\psi _{𝐢+𝐥}ϵ(𝐢)t_{𝐢,𝐢}\psi _𝐢\psi _𝐢+\psi _𝐢|\psi _𝐢|^2=0,$$
(5)
where $`\psi _𝐢=0`$ if $`𝐢`$ is outside of the sample, $`ϵ(𝐢)=4\delta _{i_x,1}\delta _{i_x,N}\delta _{i_y,1}\delta _{i_y,N}`$ and $`𝐢=(i_x=1,\mathrm{},N,i_y=1,\mathrm{},N)`$. There is one important advantage of such a treatment of the boundary condition. When neglecting the nonlinear term in the Eq.(4), the system of linear equations reduces to the problem of eigenvalues and eigenfunction of the $`hermitean`$ matrix. On the other hand, the solution of nonlinear equations requires iterations and inversion of the hermitean matrix.

Let us first discuss the solution of the linearized GLE and compare our results with previous studies . The lowest eigenvalue of the linear GLE determines the upper critical field of the sample. We have calculated eigenvalues of the linear problem expressed in units $`(a/\xi (T))^2`$ as a function of the dimensionless external magnetic field $`h=\mathrm{\Phi }/\mathrm{\Phi }_0`$ where $`\mathrm{\Phi }`$ is the total flux through the sample. Our results for a few lowest eigenvalues agree within the linewidth with the results of Ref.. Spatial pattern of the order parameter is also similar. For the field $`h5.5`$ we have observed 5 zeros of the order parameter near the center of the square. The solution corresponds to 4 vortices on the diagonals and one anti-vortex in the center of the square with total vorticity $`m=3`$. The distance between vortices is of the order of $`\delta 0.12\xi <<\xi `$. All these zeros are pinned in the region of the square where $`|\psi (x,y)|`$ is small (four orders less then the value of the order parameter near the sample boundary). This indicates that vortex-anti-vortex structure becomes unstable when we move away from the $`H_{c2}`$ line and nonlinear term and $`1/\kappa `$ corrections are considered. Moreover, the value of the order parameter and the screening current are small near zeros. In that case correction to the external field is determined by the current flowing around all 5 zeros of the order parameter rather than by the current between them. Therefor, we do not expect suppression of the field in the core of the anti-vortex.

In the following we consider the changes in the vortex structure when we move away from the $`H_{c2}`$ line. Taking into account the nonlinear term, we were unable to detect more than one zero of the order parameter near the origin of the square for the value of the field $`h=5.5`$. Consequently, the solution with more than one zero survives only very close to $`H_{c2}`$ line. We do not expect any phase transition at the point where all zeros are joined together since total flux and the symmetry of the solution do not change. This situation is different from the case of higher field $`h=7`$, where solution with four zeros of the order parameter survives far from the $`H_{c2}`$ line. In fact, the giant flux solution was not detected in that case.

In the Fig. (1a) we plot the difference in the free energy between solutions with different integral vorticity $`\mathrm{\Delta }F=F(m=3)F(m=2)`$ as a function of $`(a/\xi )^2`$ for the fixed magnetic field $`h=5.5`$. As it is clearly seen from the figure, near $`(a/\xi )^243`$ a first order phase transition takes place. At that point the high-temperature phase corresponding to the giant vortex with $`m=3`$ becomes metastable, and the phase corresponding to giant vortex, shown in Fig. (2 top), with $`m=2`$ becomes a ground state. At this point, the slope of the first derivative of the free energy as a function of $`T`$ is discontinuous, which corresponds to the latent heat of the transition. With further decrease of the temperature and $`\xi `$, the second transition takes place. At $`(a/\xi )^266`$ the giant vortex located in the center splits along one of the diagonals, Fig. (2 bottom). This transition is the second order phase transition, where $`|\psi |^2`$ is no longer invariant under four-fold axis of the square. The phase transition is clearly observed by computing the order parameter $`\eta =xy|\psi (x,y)|^2𝑑x𝑑y`$, presented in Fig. (1b). A nonvanishing $`\eta `$ is followed by a jump in the second derivative of the free energy $`\mathrm{\Delta }=d^2F/d((a/\xi )^2)^2`$, presented in Fig. (1c). We estimated the magnitude of the effect in terms of the specific heat jump: $`\mathrm{\Delta }C/a^2d=\frac{510^3\mathrm{\Phi }_0^2T_c^{^{}}}{(4\pi )^3\lambda ^2\xi ^2T_c^2}`$, where $`T_c^{^{}}`$ is the critical temperature of the transition determined by the condition $`(a/\xi )^266`$, and $`T_c`$ is the critical temperature of the sample at $`h=0`$. Landau functional, describing this phase transition, is defined as $`\mathrm{F}=\alpha _1\eta ^2+\beta _1\eta ^4`$, where $`\alpha _1`$ and $`\beta _1`$ are Landau coefficients with $`\alpha _1(TT_c^{^{}})`$. This phase transition corresponds to the one dimensional corepresentation $`B`$ of the nonunitary $`C_{4v}(C_4)`$ group.

Since we restrict our calculation to the limit of $`\kappa \mathrm{}`$ and $`\lambda a`$, long range forces ($`r\xi `$) between vortices are irrelevant. At short distances $`r\xi `$ there is short range repulsion between them. There is another interaction of vortices with the boundaries of the sample, known as Bean-Livingston barrier . In the vicinity of $`H_{c2}`$ line the interaction with the boundaries is larger than the repulsion between vortices. This is the case for the magnetic field $`h=5.5`$, where the giant vortex with $`m=3`$ is located in the center of the sample. With the decrease of $`\xi `$ the energy difference between two different solutions with different vorticity decreases and the first order phase transition to a giant flux state with vorticity $`m=2`$ takes place at $`(a/\xi )^243`$. A giant flux state (with lower vorticity) remains stable, since the interaction with the sample boundaries still dominates over the repulsion between vortices. Further decrease of $`\xi `$ leads to the decrease of the interaction of the vortices with the boundaries. Bean-Livingston force decreases and repulsion of the vortices in the giant flux state starts to dominate. As a result of the interplay of the vortex-vortex repulsion and the repulsion of vortex from the boundaries, the second order phase transition (at $`(a/\xi )^266`$) between the giant vortex state and the multivortex state takes place, preserving the integral vorticity. Separation between vortices is determined by the vortex-vortex repulsion that tends to separate them as far as possible, while in the contrary, repulsion from the boundaries prevents vortices from approaching the boundaries.

Situation is different when external magnetic field is increased to $`h7`$. In that case near $`H_{c2}`$ line ground state corresponds to a multivortex state with the total vorticity $`m=4`$, Fig. (3 top). When temperature decreases, the first order phase transition takes place at $`(a/\xi )^2=110`$, see Fig. (1a). At this point multivortex state with $`m=4`$ becomes unstable while the multivortex state with $`m=3`$, presented in Fig. (3 bottom), represents the solution with the lowest free energy. Apart from the change of the vorticity, the symmetry is also reduced at the transition point. Phase transition takes place in accordance with joint corepresentation $`E^\pm `$ of the nonunitary group $`C_{4v}(C_4)`$. Consequently, four

orientations of the pseudodipolar moment of the vortices are possible. Two component order parameter corresponding to given change of symmetry, presented in Fig. (1b), can be determined as follows: $`\eta _x=x|\psi (x,y)|^2𝑑x𝑑y`$, $`\eta _y=y|\psi (x,y)|^2𝑑x𝑑y`$. Free energy in that case depends on the vorticity $`m`$ and order parameter $`\eta _x,\eta _y`$. $`\mathrm{F}(m=4,\eta _x,\eta _y)`$ always has a minimum at $`\eta _x=\eta _y=0`$. For $`m=3`$, $`\mathrm{F}(m=3,\eta _x,\eta _y)=\alpha _2(\eta _x^2+\eta _y^2)+\beta _2(\eta _x^4+\eta _y^4)+\gamma _2\eta _x^2\eta _y^2`$, where $`\alpha _2,\beta _2,\gamma _2`$ are Landau coefficients, $`\alpha _2(TT_c^{^{}})`$ where $`T_c^{^{}}`$ is the temperature of the transition between the giant vortex state and the multivortex state for the case of $`m=3`$. $`T_c^{^{}}`$ is determined by the condition $`(a/\xi )^256`$. This transition is unobservable because $`T_c^{^{}}`$ is lower than the transition temperature for the first order phase transition where the vorticity $`m`$ changes form $`m=3`$ to $`m=2`$ (compare Figs. (1a) and (1b)).

The first order phase transition between giant vortex states with different vorticity can be described qualitatively on the basis of a simple variational function for the order parameter. Spatial dependence of the order parameter in the giant vortex state with vorticity $`m`$ can be approximated by the function:

$$\psi (r,\varphi )=\{\begin{array}{c}\hfill (r/\xi )^m\mathrm{exp}(im\varphi );r<\xi ,\\ \hfill \mathrm{exp}(im\varphi );r>\xi .\end{array}$$
(6)
Substituting this function into the Ginzburg-Landau functional and keeping only leading terms in $`(a/\xi )^2`$, we obtain a simple expression for the free energy which should be minimized as a function of vorticity $`m`$:

$$F\frac{\mathrm{\Phi }_0^2d}{(8\pi )^2\lambda ^2}(m^2\mathrm{ln}(a^2/\pi \xi ^2)/2\mathrm{\Phi }m/\mathrm{\Phi }_0).$$
(7)
Minimization of the Eq.(7) in the large $`m1`$ limit provides the expression for the vorticity:

$$m\frac{\mathrm{\Phi }}{\mathrm{\Phi }_0}/\mathrm{ln}(a^2/\pi \xi ^2).$$
(8)
Minimization of the Eq.(7) for $`h=5.5`$ provides that the phase transition from $`m=3`$ to $`m=2`$ takes place at $`(a/\xi )^2=28`$ which is lower than calculated value $`(a/\xi )^2=43`$. Nevertheless, this simple variational formula provides simple explanation of the transition with the change of vorticity in the giant vortex. It is interesting to note that due to logarithmic term, Eq.(7) for the free energy is similar to the free energy of the rotating superfluid liquid . This similarity appears because $`\lambda >a`$ and all integrals are cut at $`a`$, rather than at $`\lambda `$. Eq.(7) yields also the estimate for $`H_{c1}`$, the field at which the first vortex appears in the sample. Substituting $`m=1`$ to the Eq.(1) and solving equation $`F(m=1)=0`$ we obtain:

$$H_{c1}=\frac{\mathrm{\Phi }_0}{2a^2}\mathrm{ln}(a^2/\pi \xi ^2).$$
(9)
This expression is similar to the bulk $`H_{c1}`$ where $`\lambda ^2`$ is substituted to $`a^2/\pi `$.

In conclusion we have solved GLE in the limit of $`\kappa \mathrm{}`$ for the thin square film. We have predicted a series of the first and second order phase transitions with the change of temperature and magnetic field. On the basis of the symmetry we constructed Landau functional for the second order transitions. First order transitions are described on the basis of variational estimates.

We wish to thank D. Mihailović, V.V. Moshchalkov and A.S.Alexandrov for many useful discussions and suggestions. We are also grateful to I. Sega for critically reading the manuscript.