text:
# Dark Matter and Cosmic Strings in Particle Models11footnote 1 Invited Plenary Talk Given at COSMO’99, ICTP, Trieste, Sept.27-Oct.2, 1999

## 1  Non-thermal Production of Neutralino Cold Dark Matter from Cosmic String Decays

To begin with, we consider a general case and calculate the relic mass density of the lightest supersymmetric particle(LSP), then we will move on to a discussion of some implications. Consider a phase transition which is induced by the vacuum expectation value (vev) of some Higgs field $`\mathrm{\Phi }`$, $`|\mathrm{\Phi }|=\eta `$, and takes place at a temperature $`T_c`$ with $`T_c\eta `$. The strings are formed by the Higgs field $`\mathrm{\Phi }`$ and gauge field $`A`$. The mass per unit length of the strings is given by $`\mu =\eta ^2`$. During the phase transition, a network of strings forms, consisting of both infinite strings and cosmic string loops. After the transition, the infinite string network coarsens and more loops form from the intercommuting of infinite strings. Cosmic string loops loose their energy by emitting gravitational radiation. When the radius of a loop becomes of the order of the string width, the loop releases its final energy into a burst of $`\mathrm{\Phi }`$ and $`A`$ particles. Those particles subsequently decay into LSP, which we denote by $`\chi `$, with branching ratios $`ϵ`$ and $`ϵ^{}`$. For simplicity we now assume that all the final string energy goes into $`\mathrm{\Phi }`$ particles. A single decaying cosmic string loop thus releases $`N2\pi \lambda ^1ϵ`$ LSPs which we take to have a monochromatic distribution with energy $`E\frac{T^c}{2}`$ and $`\lambda `$ is the Higgs self coupling constant.

In such scenarios, we thus have two sources of cold dark matter(CDM) which will contribute to the matter density of the universe. We have CDM which comes from the standard scenario of thermal production; it gives a contribution to the matter density $`\mathrm{\Omega }_{therm}`$. And we also have non-thermal production of CDM which comes from the decay of cosmic string loops and gives a contribution $`\mathrm{\Omega }_{nonth}`$. The total CDM density is $`\mathrm{\Omega }_{CDM}=\mathrm{\Omega }_{therm}+\mathrm{\Omega }_{nonth}`$. During the temperature interval between $`T_c`$ and the LSP freezeout temperature $`T_\chi `$, LSPs released by decaying comic string loops will thermalise very quickly with the surrounding plasma, and hence will contribute to $`\mathrm{\Omega }_{therm}`$, which should not sensitively deviate from the value calculated by the standard method. However, below the LSP freezeout temperature, since the annihilation of the LSP is by definition negligible, each CDM particle released by cosmic string decays will contribute to $`\mathrm{\Omega }_{nonth}`$. We obviously must have

$$\mathrm{\Omega }_{nonth}<1.$$
(1)
This will lead us to a constraint (a lower bound) on the cosmic string forming scale. We now calculate $`\mathrm{\Omega }_{nonth}`$.

We assume that the strings evolve in the friction dominated regime so that the very small scale structure on the strings has not formed yet. The network of strings can then be described by a single length scale $`\xi (t)`$. In the friction dominated period, the length scale $`\xi (t)`$ has been shown to scale as :

$$\xi (t)=\xi (t_c)\left(\frac{t}{t_c}\right)^{\frac{3}{2}}$$
(2)
where $`\xi (t_c)(\lambda \eta )^1`$. The number density of cosmic string loops created per unit of time is given by :

$$\frac{dn}{dt}=\nu \xi ^4\frac{d\xi }{dt}$$
(3)
where $`\nu `$ is a constant of order 1. We are interested in loops decaying below $`T_\chi `$. The number density of LSP released from $`t_{lsp}`$ till today is given by:

$$n_{lsp}^{nonth}(t_0)=N\nu _{\xi _F}^{\xi _0}\left(\frac{t}{t_0}\right)^{\frac{3}{2}}\xi ^4𝑑\xi $$
(4)
where the subscript $`0`$ refers to parameters which are evaluated today. $`\xi _F=\xi (t_F)`$ where $`t_F`$ is the time at which cosmic string loops which are decaying at time $`t_\chi `$ (associated with the LSP freezeout temperature $`T_\chi `$) have formed. Now the loop’s average radius shrinks at a rate  $`\frac{dR}{dt}=\mathrm{\Gamma }_{loops}G\mu `$, where $`\mathrm{\Gamma }_{loops}`$ is a numerical factor $`1020`$. Since loops form at time $`t_F`$ with an average radius $`R(t_F)\lambda ^1G\mu M_{pl}^{\frac{1}{2}}t_F^{\frac{3}{2}}`$, they have shrunk to a point at the time $`t\lambda ^1\mathrm{\Gamma }_{loops}^1M_{pl}^{\frac{1}{2}}t_F^{\frac{3}{2}}`$. Thus $`t_F(\lambda \mathrm{\Gamma })_{loops}^{\frac{2}{3}}M_{pl}^{\frac{1}{3}}t_\chi ^{\frac{2}{3}}`$. Now the entropy density is $`s=\frac{2\pi ^2}{45}g_{}T^3`$ where $`g_{}`$ counts the number of massless degrees of freedom in the corresponding phase. The time $`t`$ and temperature $`T`$ are related by $`t=0.3g_{}^{\frac{1}{2}}(T)\frac{M_{pl}}{T^2}`$ where $`M_{pl}`$ is the Planck mass. Thus using Eqs.(2) and (4), we find that the LSP number density today released by decaying cosmic string loops is given by:

$$Y_{LSP}^{nonth}=\frac{n_{lsp}^{nonth}}{s}=\frac{6.75}{\pi }ϵ\nu \lambda ^2\mathrm{\Gamma }_{loops}^2g_{_{T_c}}^{\frac{9}{4}}g_{_{T_\chi }}^{\frac{3}{4}}M_{pl}^2\frac{T_\chi ^4}{T_c^6},$$
(5)
where the subscript on $`g^{}`$ refers to the time when $`g^{}`$ is evaluated.

The LSP relic abundance is related to $`Y_\chi `$ by:

$`\mathrm{\Omega }_\chi h^2`$ $``$ $`M_\chi Y_\chi s(t_0)\rho _c(t_0)^1h^2`$ (6)
$``$ $`2.82\times 10^8Y_\chi ^{tot}(M_\chi /\mathrm{GeV})`$
where $`h`$ is the Hubble parameter and $`M_\chi `$ is the LSP mass. Now $`Y_{LSP}^{tot}=Y_\chi ^{therm}+Y_\chi ^{nonth}`$; hence by setting $`h=0.70`$, Eqs. (6) and (1) lead to the following constraint:

$$5.75\times 10^8Y_\chi ^{nonth}(M_\chi /\mathrm{GeV})<1.$$
(7)
We thus see that Eqs. (5) and (7) lead to a lower bound on the cosmic string forming temperature $`T_c`$.

Our results have important implications for supersymmetric extensions of the standard model with extra $`U(1)`$’s (or grand unified models with an intermediate $`SU(3)_c\times SU(2)_L\times U(1)_Y\times U(1)^{}`$ gauge symmetry). Most importantly, the requirement $`\mathrm{\Omega }_{nonth}<1`$ imposes a new constraint on supersymmetric model building and rules out many models with a low scale of a new symmetry breaking which produces defects such as cosmic strings.

Consider, for example, the model with an extra $`U_{BL}(1)`$ gauge symmetry. In this model, the strings will release not only right-handed neutrinos $`N_i`$, but also their superpartners $`\stackrel{~}{N}_i`$. The heavy neutrinos $`N_i`$ and their scalar partners $`\stackrel{~}{N}_i`$ can decay into various final states including the LSP. The superpotential relevant to the decays is

$$W=H_1ϵLy_lE^c+H_2ϵLy_\nu N^c,$$

where $`H_1,H_2,L,E^c`$ and $`N^c`$ are the chiral superfields and $`y_l,y_\nu `$ are Yukawa couplings for the lepton and neutrino Dirac masses, $`m_l=y_lv_1,m_D=y_\nu v_2`$, with $`v_{1,2}`$ being the vacuum expectation values of the Higgs fields. At tree level, the decay rates of $`N_i`$ into s-lepton plus Higgsino and lepton plus Higgs are the same and they are smaller than the rate of $`\stackrel{~}{N}_i`$ decaying into s-lepton plus Higgs and Higssino plus lepton by a factor of 2. If the neutralino is higgsino-like, the LSP arise directly from the decays of the $`N_i`$ and $`\stackrel{~}{N}_i`$. If the neutralino is bino- or photino-like, subsequent decays of s-lepton into binos or photinos plus leptons will produce the LSP. For reasonable values of the parameters, we estimate the branching ratio $`ϵ`$ of the heavy particle decay into LSP to be between 0.1 and 0.5. From Eq. (5) it follows that string decays can easily produce the required amount of LSP. However too many LSPs will be generated unless the $`BL`$ breaking scale, $`\mathrm{\Lambda }_{BL}`$ is higher than about $`10^8`$ GeV . In turn, this will set a lower limit on the neutrino masses generated by the see-saw mechanism, $`m_\nu m_D^2/\mathrm{\Lambda }_{BL}`$. Inserting numbers and taking $`m_Dm_\tau 1.8`$ GeV, one obtains that $`m_\nu 30`$ eV.

Our lower limit on the $`BL`$ symmetry breaking scale in gauged $`BL`$ models and in general models with an extra $`U(1)`$  pushes the mass of the new gauge boson far above the Fermi scale, rendering it impossible to test the new physics signals from the extra $`Z^{}`$ in accelerators.

As one more implication of our results, we consider the limit on the lifetime of the Z-string in the MSSM. Since the Z-string is produced during the electroweak phase transition, the stable Z-string would produce too much LSP and overclose the universe. However if the Z-string decays before the temperature of the LSP freezeout, then the LSP produced from the decay of the Z-string would be thermalized immediately, which would result in a negligible $`\mathrm{\Omega }_{nonth}`$.

## 2 $`\pi `$ and $`\eta ^{}`$ strings in QCD

We consider an idealization of QCD with two species of massless quarks $`u`$ and $`d`$. The lagrangian of strong interaction physics is invariant under $`SU_L(2)\times SU_R(2)`$ (we will come back to the discussion of the $`U_A(1)`$ at end of this section) chiral transformations

$$\mathrm{\Psi }_{L,R}\mathrm{exp}(i\stackrel{}{\theta }_{L,R}\stackrel{}{\tau })\mathrm{\Psi }_{L,R},$$
(8)
where $`\mathrm{\Psi }_{L,R}^T=(u,d)_{L,R}`$. However this chiral symmetry does not appear in the low energy particle spectrum since it is spontaneously broken due to the quark condensate. Consequently, three Goldstone bosons, the pions, appear and the (constituent) quarks become massive. At low energy, the spontaneous breaking of chiral symmetry can be described by an effective theory, the linear sigma model, which involves the massless pions $`\stackrel{}{\pi }`$ and a massive $`\sigma `$ particle.

As usual, we introduce the field

$$\mathrm{\Phi }=\sigma \frac{\tau ^0}{2}+i\stackrel{}{\pi }\frac{\stackrel{}{\tau }}{2},$$
(9)
where $`\tau ^0`$ is unity matrix and $`\stackrel{}{\tau }`$ are the Pauli matrices with the normalization condition $`Tr(\tau ^a\tau ^b)=2\delta ^{ab}`$. Under $`SU_L(2)\times SU_R(2)`$ chiral transformations, $`\mathrm{\Phi }`$ transforms as

$$\mathrm{\Phi }L^+\mathrm{\Phi }R.$$
(10)
The renormalizable effective lagrangian of the linear sigma model is given by

$$=_\mathrm{\Phi }+_q,$$
(11)
where

$$_\mathrm{\Phi }=Tr[(_\mu \mathrm{\Phi })^+^\mu \mathrm{\Phi }]\lambda [Tr(\mathrm{\Phi }^+\mathrm{\Phi })\frac{f_\pi ^2}{2}]^2,$$
(12)
and

$$_q=\overline{\mathrm{\Psi }}_Li\gamma ^\mu _\mu \mathrm{\Psi }_L+\overline{\mathrm{\Psi }}_Ri\gamma ^\mu _\mu \mathrm{\Psi }_R2g\overline{\mathrm{\Psi }}_L\mathrm{\Phi }\mathrm{\Psi }_R+h.c..$$
(13)
During chiral symmetry breaking, the field $`\sigma `$ takes on a nonvanishing vacuum expectation value, which breaks $`SU_L(2)\times SU_R(2)`$ down to $`SU_{L+R}(2)`$. This results in a massive sigma $`\sigma `$ and three massless Goldstone bosons $`\stackrel{}{\pi }`$, as well as giving a mass $`m_q=gf_\pi `$ to the constituent quarks. Numerically, $`f_\pi 94`$ MeV and $`m_q300`$ MeV.

We studied the classical solutions of this model and discovered a class of vortex-like configurations which we refer to as pion string. Like the Z string of the standard electroweak model, the pion string is not topologically stable. Nevertheless, as demonstrated recently in numerical simulations in the case of semilocal strings (which are also not topologically stable) pion strings are expected to be produced during the QCD phase transition in the early Universe (and also in heavy-ion collisions). The strings will subsequently decay.

The pion string is a static configuration of the lagrangian $`_\mathrm{\Phi }`$ of Eq. (12). To construct these solutions, we define new fields

$`\varphi `$ $`=`$ $`{\displaystyle \frac{\sigma +i\pi ^0}{\sqrt{2}}},`$ (14)
$`\pi ^\pm `$ $`=`$ $`{\displaystyle \frac{\pi ^1\pm i\pi ^2}{\sqrt{2}}}.`$ (15)
The lagrangian $`_\mathrm{\Phi }`$ now can be rewritten as

$$=(_\mu \varphi )^{}^\mu \varphi +_\mu \pi ^+^\mu \pi ^{}\lambda (\pi ^+\pi ^{}+\varphi ^{}\varphi \frac{f_\pi ^2}{2})^2.$$
(16)
For static configurations, the energy functional corresponding to the above lagrangian is given by

$$E=d^3x\left[\stackrel{}{}\varphi ^{}\stackrel{}{}\varphi +\stackrel{}{}\pi ^+\stackrel{}{}\pi ^{}+\lambda (\pi ^+\pi ^{}+\varphi ^{}\varphi \frac{f_\pi ^2}{2})^2\right].$$
(17)
The time independent equations of motion are:

$`^2\varphi `$ $`=`$ $`2\lambda (\pi ^+\pi ^{}+\varphi ^{}\varphi {\displaystyle \frac{f_\pi ^2}{2}})\varphi ,`$ (18)
$`^2\pi ^+`$ $`=`$ $`2\lambda (\pi ^+\pi ^{}+\varphi ^{}\varphi {\displaystyle \frac{f_\pi ^2}{2}})\pi ^+.`$ (19)
The pion string solution extremising the energy functional of Eq. (17) is given by

$`\varphi `$ $`=`$ $`{\displaystyle \frac{f_\pi }{\sqrt{2}}}\rho (r)e^{i\theta },`$ (20)
$`\pi ^\pm `$ $`=`$ $`0,`$ (21)
where the coordinates $`r`$ and $`\theta `$ are polar coordinates in the $`xy`$ plane (the string is assumed to lie along the z axis), and $`\rho (r)`$ satisfies the following boundary conditions

$`r0,`$ $`\rho (r)0;`$ (22)
$`r\mathrm{},`$ $`\rho (r)1.`$ (23)
The pion-string is not topologically stable, since any field configuration can be continuously deformed to the vacuum.

In the discussions above we have neglected the electromagnetic interaction. When turning on $`U_{em}(1)`$, the derivatives are replaced by the covariant derivatives. Recently Nagasawa and Brandenberger have pointed out that the interactions with a finite temperature plasma will lead to corrections to the effective potential and may stabilize the $`\pi `$ strings.

Since a pion string is made of $`\sigma `$ and $`\pi ^0`$ fields, it is neutral under the $`U_{em}(1)`$ symmetry. However, the $`\pi ^0`$ will couple to photons via the Wess Zumino type interaction.

$$_{low}=\frac{f_\pi ^2}{4}Tr(_\mu \mathrm{\Sigma }^+^\mu \mathrm{\Sigma })\frac{1}{4}F_{\mu \nu }F^{\mu \nu }\frac{N_c\alpha }{24\pi }\frac{\pi ^0}{f_\pi }ϵ^{\mu \nu \alpha \beta }F_{\mu \nu }F_{\alpha \beta },$$
(24)
where $`N_c=3`$, $`\mathrm{\Sigma }=exp(i\stackrel{}{\tau }\stackrel{}{\pi }/f_\pi )`$, and $`\alpha `$ is the electromagnetic fine structure constant.

Recently Brandenberger and I propose a new mechanism for the generation of primordial magnetic fields. It is based on the realization that anomalous global strings couple to electricity and magnetism via the anomalous interactions shown above and induce magnetic fields . The major advantage of this mechanism is that the coherence scale of the magnetic fields induced by these global vortex lines is set by the length scale $`\xi (t)`$ of the strings (the typical curvature radius of the strings).

Before concluding this section, I will speculate about the existence of an unstable $`\eta ^{}`$ string. In QCD, in the limit of massless quarks, there is an additional $`U_A(1)`$ chiral symmetry. This chiral symmetry, when broken by the quark condensate, predicts the existence of a goldstone boson. There is no such a light meson, however. This is resolved by the Adler-Bell-Jackiw $`U_A(1)`$ anomaly together with the properties of non-trivial vacuum structure of non-abelian gauge theory, QCD. The $`U_A(1)`$ symmetry is badly broken by instanton effects at zero temperature.

As the density of matter and/or the temperature increases, it is expected that the instanton effects will rapidly disappear, and one thus has an additional $`U_A(1)`$ symmetry (besides $`SU_L(2)\times SU_R(2)`$) at the transition temperature of the QCD chiral symmetry. When the $`U_A(1)`$ symmetry is broken spontaneously by the quark condensate, a topological string, the $`\eta ^{}`$-string, results. Differing from the pion-string, the $`\eta ^{}`$-string is topologically stable at high temperatures, but will decay as the temperature decreases. The $`\eta ^{}`$-string can form during the chiral phase transition of QCD. In the setting of cosmology, it will exist during a specific epoch below the QCD chiral symmetry breaking temperature during the evolution of the universe. In the context of heavy-ion collisions, it will exist in the plasma created by the collision. The strings then become unstable as the temperature decreases and when the instanton effects become substantial.

## Acknowledgments

This work was supported in part by the NSF of China.


format: text

text:
# Galactic Halo substructure in the Sloan Digital Sky Survey: the ancient tidal stream from the Sagittarius dwarf galaxy

## 1 Introduction

In the standard hierarchical picture of galaxy formation, galaxies like the Milky Way grew by repeated mergers with other galaxies, and in particular with many smaller galaxy “fragments”. In the outer regions of galactic halos, disruption timescales can be long, and we can hope to observe galaxy remnants as well as to catch galaxies in the process of being accreted. That this accretion is still an on-going process is clear, as evidenced by the many studies of tidal features in the dwarf satellite galaxies and globular clusters.

However, what is of prime interest is to compare the observed structures to the results of galaxy formation theories. These now make strong predictions about the structure and substructure of galaxy halos. Indeed the standard cold dark matter models appear to be inconsistent with observational constraints, as they over-predict halo clumpiness (see, e.g., Moore et. al 1999).

In an attempt to constrain such theories, several groups have undertaken large scale surveys to uncover substructure in the halo. These include, among others, the “spaghetti” survey (Morrison et al., 2000), the APM carbon star survey (Totten & Irwin, 1998; Totten, Irwin & Whitelock, 2000; Ibata et al., 2000), and the SDSS survey (Yanny et al., 2000; Ivezić et al., 2000).

The APM carbon star survey, an almost all high-latitude sky survey, revealed the existence of a gigantic band of these intermediate age ($`6\mathrm{Gyr}`$) stars around the sky (Ibata et al. 2000; hereafter paper 1). 111Extant optical and NIR photometry and detailed optical spectroscopy of the APM colour-selected cool Halo carbon stars show that these stars are identical in nature to those found in AGB populations in other Local Group galaxies, particuarly the Magellanic Clouds and many of the dwarf spheroidal Galactic satellites (Totten & Irwin, 1998; Totten, Irwin & Whitelock, 2000). Furthermore the central part of Sgr has a known AGB carbon star population (Whitelock, Irwin & Catchpole, 1996) that also has NIR and optical properties indistinguishable from the APM Halo sample. Normal AGB carbon stars are carbon-enriched during dredge-ups driven by helium burning thermal pulses during the post main sequence phase of 1-2 solar mass stars. As such they are thought to be $`>12`$ Gyrs old but less than 10 Gyrs old (see e.g. the reviews by Da Costa 1997 and also Azzopardi & Lequeux 1992). Approximately half of the “Halo” carbon stars detected in this survey are located in this ten degree wide band, which is centered on the expected orbit of the Sagittarius dwarf galaxy in a spherical Halo potential (paper 1). Sagittarius, the closest dwarf satellite galaxy to the Milky Way (Ibata, Gilmore & Irwin, 1994), is clearly in the process of being tidally torn apart. As shown in paper 1, those carbon stars trace a gigantic tidal stream of debris disrupted off the Sagittarius dwarf galaxy which has now wrapped itself around the Halo. Interestingly, this stream is observed to follow a great circle on the sky, which means that it is probably a coplanar rosette-shaped structure that we observe from within (since it should trace relatively closely the orbit of the center of mass of the dwarf galaxy). Most importantly, however, the fact that the stream is confined to a plane, implies that it does not suffer significant precession, which is only possible if the Galactic halo potential is not significantly oblate.

Two SDSS teams have recently presented results from their first year’s commissioning data (Yanny et al., 2000; Ivezić et al., 2000). Due to initial technical constraints, the drift-scanning camera took data only on the celestial equator during the commissioning period. Their data consists of two wedges: a NGC sample $`87^{}`$ long and a SGC sample $`60^{}`$ long, with each covering a declination range of only $`1.25^{}\stackrel{<}{}\delta \stackrel{<}{}1.25^{}`$. Despite the fact that their data provide essentially a two dimensional slice through the Galaxy, they find astonishing sub-structure in the distribution of A-stars. Sharply-defined concentric rings of stars are detected. The number of RR-lyrae variables also cuts off sharply at $`50\mathrm{kpc}`$. The nature of these structures is important, as they are the dominant structure seen in old Halo stars, comprising approximately half of the total number of Halo stars seen in the survey. In this contribution, we show that these A-star Halo sub-structures are morphologically very similar to the expected structure of the tidal debris of the Sagittarius dwarf galaxy, are at similar distances and are highly likely to be part of the Sagittarius stream.

At first sight it may seem surprising that such a limited survey can detect Halo sub-structure. However, since most long-lived Galactic satellite systems appear to be on close to polar orbits the chances of a stellar debris stream intersecting the current SDSS zero-declination survey strips is $``$50%. The Sagittarius stellar stream is currently the only known large relatively intact debris system and is within $`13^{}`$ of a polar orbit.

## 2 The expected debris of the Sagittarius dwarf galaxy

In paper 1, we undertook a large suite of numerical simulations of the interaction between the Sagittarius dwarf galaxy and the Milky Way, in order to determine how the APM carbon star sample constrains the shape of the Galactic halo. Figure 1 shows a typical high resolution simulation with $`10^5`$ particles, at the end of the integration, at $`T=12\mathrm{Gyr}`$ (details of these simulations are given in paper 1). For this particular simulation, the initial structure of the dwarf galaxy is a fluffy King model with mass $`10^8\mathrm{M}_{}`$, half-mass radius of $`r_{1/2}=0.5\mathrm{kpc}`$ and initial concentration $`c=0.5`$. The Milky Way model (termed model H1 in paper 1) has a bulge, a thin disk, a thick disk and an interstellar medium, with parameters identical to the Dehnen & Binney (1998) model 2; the Galactic halo has a flat rotation curve, with a constant density core of radius $`r_0=3\mathrm{kpc}`$, and an oblate spheroidal density structure, with density flattening of $`q_m=0.9`$. The total mass of the Galaxy model is such that the circular velocity at $`50\mathrm{kpc}`$ is $`v_c=220\mathrm{km}\mathrm{s}^1`$.

After accounting for a $`50`$% “background” from other accretion events, such simulations can give a reasonable representation of the observed carbon star distribution, both in spatial position and in kinematics (cf. paper 1). The Aitoff projection of this simulation onto the sky (Figure 1), shows that the majority of the particles in this model are clustered close to a very obvious great circle stream. However, the potential is not quite spherical, so orbits suffer some precession. The material first stripped from the dwarf galaxy has had longest to drift away from its progenitor, and so also experiences the largest difference in precession. The particles not on the great stream are those that were stripped from the dwarf galaxy at earlier times; indeed, the precession away from the great circle stream effectively provides a chronometer of the time since particles were disrupted from their parent galaxy.

To compare the simulation of Figure 1 to the SDSS commissioning run data, we take a $`2.5^{}`$ equatorial cut though the simulation. As discussed by Yanny et al. (2000) and Ivezić et al. (2000), there may be some mixture of blue star populations in their dataset, however, for simplicity, we assume that they are all blue horizontal branch stars with absolute magnitude $`M_g^{}=1.0`$. The resulting distribution of these equatorial stars as a function of right ascension and apparent $`g^{}`$ magnitude is displayed in Figure 2. The resemblance of this distribution with that found by Yanny et al. (2000) (their Figures 3 and 11) is striking. This simulation (which preceded the SDSS result) predicts the azimuthal rings seen by the Sloan Survey. Figure 3 is identical to Figure 2, but the simulation it is derived from was undertaken in a more flattened halo, with $`q_m=0.5`$. The sharp ring features seen in Figure 2 and in the SDSS dataset are smeared-out in Figure 3. A spherical Halo ($`q_m=1.0`$) gives rise to a slower precession rate (there is still some precession due to the modelled disk component), so a smaller fraction of particles wander off the great circle stream. Unfortunately, we cannot be quantitative in our comparisons to the models without a published list of SDSS A-stars.

The color-coding of Figures 2 and 3 indicates the time elapsed since the particles left the tidal radius of their progenitor. In Figure 2, the tight clumps seen near 2 hours and 14 hours are due to recently disrupted material, which can be composed of young as well as old stars. This is consistent with the detection of a clump of carbon stars in these regions (paper 1). The particles that were disrupted long ago, displayed here in yellow and red, have had many Gyrs to precess away from the great circle stream. The only remaining visible tracers of this material will be old stars, such as RRLyrae variables. This model suggests that the ancient disruption and star-formation history of the Sagittarius dwarf may be traced and analysed with the population of RRLyrae variables first detected by SDSS. In the main body of the Sagittarius dwarf galaxy, blue horizontal branch stars make up a much higher fraction of the total stellar population in its globular clusters than in its field stars (Layden & Sarajedini, 2000). The fact that the older parts of its tidal stream display apparently high concentrations of RRLyrae stars suggests a common origin with the Sagittarius dwarf’s globular clusters (though this may just reflect continual star-formation in the dwarf galaxy).

Our simulations make strong predictions. First, the radial velocity of the A-stars should vary smoothly as a function of right ascension, as displayed on the top panel of Figure 1. The radial velocity dispersion of these stars in a small range in right ascension should be small also, $`20\mathrm{km}\mathrm{s}^1`$, except in regions where multiple streams cross each other, at $`\alpha 2`$ hours and $`\alpha 14`$ hours. Furthermore, the majority of the A-stars should have similar radial velocities to those of the carbon stars in the 2 hour and 14 hour regions of paper 1. On the equatorial plane of the sky (where the extant SDSS data lie) at $`\alpha 2`$ hours, the velocity distribution should be bimodal, with two clumps at Heliocentric radial velocity of $`v+200\mathrm{km}\mathrm{s}^1`$ and $`v150\mathrm{km}\mathrm{s}^1`$. Further along the equatorial plane at $`\alpha 14`$ hours the mean Heliocentric velocity of the stream stars should be $`50\mathrm{km}\mathrm{s}^1`$, with significant dispersion of $`\sigma 60\mathrm{km}\mathrm{s}^1`$. Secondly, the distances of the A-stars, and RR-lyraes, should also closely follow the range of distances found from the C-star studies. In particular, an upper cutoff of $``$50 kpc in the RR-lyrae distances is exactly what is seen in the C-star distribution.

## 3 Conclusions

We have shown that numerical simulations of the dynamical evolution of the Sagittarius dwarf galaxy, previously reported in paper 1, predict the presence of sharp shell-like features in the Galactic halo. This shell-like structure is due to the effect of orbital precession acting on stars which were disrupted from the Sagittarius dwarf galaxy many billion years ago. Making a cut through our simulated dataset along the equatorial plane of the sky, we find a distribution very similar to that found in A-stars in the SDSS data (Yanny et al., 2000; Ivezić et al., 2000).

The A-stars furthest away from the dominant stream detected in paper 1, belong to those populations that were first disrupted from the dwarf galaxy. By comparing the relative densities of young to old stars in these tidal structures, it may be possible to reconstruct the dynamical evolution of the Sagittarius dwarf galaxy stream in great detail (Johnston et al., 1999).

The SDSS data that will soon become available along other scan-lines on the sky offer the very useful opportunity to distinguish between the different dynamical models of the Sagittarius dwarf galaxy. Such data will also allow us to distinguish competing models for the mass distribution of the dark halo. For instance, a comparison of the particle distributions of Figures 2 and 3, shows that the enhanced stream precession in the flattened halo of Figure 3 gives rise to a substantial density between $`270^{}<\alpha <0^{}`$, which is not seen in Figure 2. Data in that right ascension range will provide a means to further constrain the shape of the Galactic halo.

With the almost all high-latitude sky coverage of the APM carbon star survey, only two Halo streams were found (the other probably being due to the Large Magellanic Cloud). Given that SDSS has also intersected the tidal stream from the Sagittarius dwarf galaxy, this stream appears to be by far the dominant contributor to Halo substructure at Galactocentric distances less than $`50\mathrm{kpc}`$. It seems unlikely, therefore, that any more accreted systems of such magnitude will be found by SDSS or other surveys. However, there is hope to discover the much smaller accreted galaxy fragments, if they exist.

The simulation result of Figure 1, has a cautionary tale to tell about using starcounts to determine the shape of the “Halo”. That a substantial fraction of the outer Galactic halo stars appear to be clumped on an inclined band on the sky clearly does not mean that the Halo is a thin inclined flattened band! Such structures just tell a complicated story about the process of stellar accretion into the Halo. Since stars only contribute a miniscule amount to the total halo mass, the only way to reveal the dark matter distribution is through dynamical analyses. For this purpose, the dynamically cold streams discussed in this paper are ideal.


text:
# Glassy dynamics near zero temperature

## Abstract

We numerically study finite-dimensional spin glasses at low and zero temperature, finding evidences for (i) strong time/space heterogeneities, (ii) spontaneous time scale separation and (iii) power law distributions of flipping times. Using zero temperature dynamics we study blocking, clustering and persistence phenomena.

The study of non-equilibrium dynamics in quenched spin systems provides a rigorous and numerical test ground for modeling glassy behavior. Among the open issues which arise in this context there is the connection between spatial and temporal heterogeneities and the global off-equilibrium dynamical features characteristic of glassy structures such as aging . In the last decade, much progress has been made in the study of domain evolution and aging in both ordered and disordered systems . In particular, the mean-field picture of the glassy transition  plays a central role in interpreting results of numerical simulations on spin systems, but it can not account for heterogeneities. On the other hand a set of interesting numerical and analytical results on the off-equilibrium dynamics in non-frustrated systems have been obtained , which can be now extended to frustrated ones.

The scope of this study is to provide a further link between the different aspects of the glassy transition unfolded by the various approaches. More specifically, we shall focus on models with discrete ($`\pm J`$) couplings in 2 and 3 dimensions. Of particular interest is the onset of time scale separation in the low temperature (low $`T`$) regime and its connections with the blocking and persistence phenomena in the zero-temperature ($`T=0`$) limit.

The choice of discrete couplings stems from the fact that they naturally arise as interactions incorporating geometrical frustration and constraints. Moreover, models with discrete couplings are of central relevance in other areas such as combinatorial optimization and are known to provide very rich dynamical features even at $`T=0`$ .

In what follows, we report results from extensive numerical simulations concerning the following issues.

I) The onset and the dependence on waiting time of the spontaneous time scale separation and space clustering at sufficiently low $`T`$ in 2D and 3D models.

II) The related emergence at $`T=0`$ of the “blocking” phenomenon together with the clustering of a finite fraction of spins that flip infinitely often (“fast” spins), along with its interpretation in terms of the functional mean-field order parameter.

III) The $`T=0`$ distributions of flipping times and fast spins clusters sizes.

IV) The persistence phenomenon.

Our studies focus on the $`\pm J`$ Edwards-Anderson (EA) model in 2 and 3 spatial dimensions (square and cubic lattices). The Hamiltonian of these models reads

$$=\underset{(i,j)E}{}s_iJ_{ij}s_j,$$
(1)
where $`E`$ is the set of lattice edges, the spins are of Ising type and the couplings take the values $`\pm 1`$ with equal probability.

The study of 3D EA model in the very low temperature region requires huge thermalization times and therefore it has been severely limited by available computer resources. Only very recently these studies have been pushed to very low temperatures, thanks to the use of parallel computers . The data presented in this letter referring to the 3D case has been obtained with the help of the APE100 parallel computer . The 2D case is indeed much simpler in virtue of the existence of polynomial algorithms for ground-states calculations.

It is known that in both 2D and 3D EA models  lowering temperature produces a surprising increase in the number of high frequency flipping spins. Here we push the study of the flipping times distribution to the very low temperature region in the hard case of the 3D EA model. Moreover we investigate the dependence of this distribution on the waiting time, which is a relevant feature of the aging regime in glassy phases.

According to the typical scheme used in off-equilibrium dynamical studies , we simulate large systems (of at least $`32^3`$ spins) and we start the experiment with an instantaneous quench from infinite temperature to one in the glassy phase ($`T<T_c1.1`$). Next we let the system evolve for a waiting time $`t_w`$ (times are expressed in terms of Monte Carlo sweeps, MCS). Finally we calculate the flipping rates probability distribution function (pdf) measuring the number of flips done by every spin within time windows extending from $`t_w+t`$ to $`t_w+2t`$, where $`t=2^k`$ (with $`k26`$) is also the time window size. We simply define the mean flipping time $`\tau `$ as the time window size divided by the number of flips and we construct the pdf of its logarithm, $`P_{t,t_w}(\mathrm{ln}\tau )`$, by taking the histogram over all spins (with the label $`t`$ identifying the time window). This distribution is expected to be self-averaging (like correlation functions) and so we prefer to simulate few very large samples. The choice of working with the $`\mathrm{ln}(\tau )`$ pdf instead of that for $`\tau `$ is dictated by the broadness of the latter. Note that a pure exponential tail in the former \[$`P(\mathrm{ln}\tau )e^{\lambda \mathrm{ln}\tau }`$\] corresponds to a power law tail in the latter \[$`P(\tau )\tau ^{\lambda 1}`$\].

In the $`t_w=0`$ case we find, as expected, that $`P_{t,0}(\mathrm{ln}\tau )`$ does not depend on $`t`$ and so we consider only the pdf measured during the last time window, which has the largest support (more than $`10^6`$ MCS) and the highest statistics. In Fig. 1 we show such a distribution for different temperatures. While at very high temperatures the pdf is strongly peaked around the mean, it acquires very large tails approaching the glass transition temperature  and finally develops a clear bimodal shape at very low temperatures, as can be seen in Fig. 1 (note the logarithmic scale on the $`y`$ axis). The emergence of a spontaneous time scales separation, allows us to naturally divide the spins into two very general classes: the “fast” spins belonging to the left peak and the “slow” ones to the right peak. While the shape of the left part of the distribution does not change with temperature, the position of the “slow” peak is clearly moving towards higher values when the temperature is decreased. As shown in the inset of Fig. 1, this process verify the simple scaling $`T\mathrm{ln}(\tau )`$, commonly found in every activated process in spin glasses. It follows that slow spins are the ones that have to overcome a barrier in order to flip, while the fast ones will eventually have a zero local field at some time (which could even happen very rarely). In the zero temperature limit we expect that the slow peak moves towards unreachable time-scales, whereas fast spins are the only responsible for the dynamics. The adjective “fast” could be somehow misleading, in turn their actual flipping times follow a very broad distribution.

It is well known that in the aging regime two-times quantities depend on both times and not only on their difference. In the case of our $`P_{t,t_w}(\mathrm{ln}\tau )`$ this dependence on the waiting time is shown in Fig. 2, where we present data for a low temperature ($`T=0.35`$), a very large waiting time ($`t_w=2^{24}`$) and many values of $`t`$. During the aging process two different regimes can be identified : the quasi-equilibrium regime ($`tt_w`$) where the system relaxes inside a quasi-state and the aging regime ($`tt_w`$) where macroscopic rearrangements take place. In these two regimes the shape of the $`P_{t,t_w}(\mathrm{ln}\tau )`$ is different (this has been verified in many simulations with different $`t_w`$ values, even if here we report the results concerning only one waiting time). In Fig. 2 we display the pdf for different choices of the ratio $`t/t_w`$: the lowest curves correspond to the quasi-equilibrium regime, while the upper ones have values of $`t`$ of the order of $`t_w`$ and therefore include the effects coming from the aging regime.

It is clear that the time-scales separation is sharper in the quasi-equilibrium regime, while approaching the aging regime the gap between fast and slow peaks is partially filled. Observing that the slow spins create the environment where the fast spins move in (the so called cage effect), we may assume that the meta-stable state is identified by the spatial structure of slow spins which act as a rigid backbone on time scales typical of the quasi-equilibrium regime. Therefore we may also observe that the basic excitations that bring the system out of quasi-equilibrium are given by those spins which fill the gap between the peaks. The study of the spatial structure of these spins would help towards the understanding of basic excitation in 3D spin glasses.

The same set of simulations with Gaussian couplings gives no evidence of a so simple time scale separation. In this case the distribution of microscopic time scales is different and much longer simulations would be needed in order to unveil the separation. Moreover at $`T=0`$ in the Gaussian case the system rapidly get stuck in a local minimum and its configuration is fully frozen, consistently with the zero-temperature classification given in . Only the $`\pm J`$ model has a rich aging behavior at zero temperature.

In the $`T=0`$ limit activated processes disappear and the slow peak goes to infinity. Then the gap between the peaks can no longer be filled when the system leaves the quasi-equilibrium regime. In other words, the $`P_{t,t_w}(\mathrm{ln}\tau )`$ becomes $`t_w`$-independent (we have checked numerically this fact) and it describes only the fast spin component. Remarkably enough, the $`T=0`$ pdf is already present in the finite temperature data, in their quasi-equilibrium regime (see Fig. 2) and this would suggest that $`T=0`$ dynamical features could be present at low temperatures too, at least in the quasi-equilibrium regime. Hereafter we will consider only the following updating rule. A site is randomly chosen and it is oriented in the direction of the neighbors majority or at random if the local field on it is zero.

We find that the $`T=0`$ flipping time distributions can be well fitted by a power law over many decades, with an exponent $`\lambda +1=1.76(2)`$ in 3D and $`\lambda +1=2.08(3)`$ in 2D. This means that in 2D after a large but finite time (given by the average flipping time) an extensive fraction of the system can be flipped with zero-energy costs. On the contrary in 3D an average flipping time can not be defined and this divergence is maybe related to the existence of a finite temperature phase transition.

Persistence  is one the most studied properties of zero temperature dynamics , nevertheless it has been never measured in spin glasses. Here we fill this gap.

The persistence $`U(t,t_w)`$ is defined as the number of unflipped spins in the time interval $`[t_w,t_w+t]`$. In pure ferromagnetic models (at least for $`D=2,3`$) the persistence $`U(t,t_w)`$ decays to zero, in the large times limit, as $`U(t,t_w)(t/t_w)^{\theta (D)}`$, the exponent being $`t_w`$-independent. On the contrary, in disordered models the persistence is expected to remain finite  in the $`t\mathrm{}`$ limit and its asymptotic value $`U_{\mathrm{}}(t_w)`$ may depend on the waiting time. The physical mechanism underlying such effect is the slow freezing of a rigid component. It follows that the proper way of estimating the persistence in disordered systems is by taking a sufficiently large value of $`t_w`$, differently from what has been done in previous studies where the choice $`t_w=0`$ was used.

In Fig. 3 we show the results for 3D $`\pm J`$ spin glasses (very similar results have been found in the 2D case). While the system relaxes towards the stationary state the number of frozen spins grows and $`U_{\mathrm{}}(t_w)`$ monotonically increases with $`t_w`$. In the large times limit we can extract the $`\theta `$ exponent from the decay $`U(t,t_w)U_{\mathrm{}}(t_w)t^\theta `$. This fitting procedure is not an easy one, because of the slow persistence decay, and so our results are affected by large errors. Our best estimations are $`\theta =0.46(3)`$ and $`U_{\mathrm{}}(\mathrm{})=0.862(2)`$ in 3D and $`\theta =0.64(3)`$ and $`U_{\mathrm{}}(\mathrm{})=0.778(2)`$ in 2D. The $`\theta `$ value in 3D is confirmed by the local magnetizations pdf which can be well fitted by the law $`(1m^2)^{\theta 1}`$ at low and zero temperatures.

These numbers imply that in the large times limit a $`\pm J`$ spin glass is not completely frozen. While a large fraction $`U_{\mathrm{}}(\mathrm{})`$ of the spins get blocked forever, a finite fraction $`1U_{\mathrm{}}(\mathrm{})`$ of the volume is composed by spins that flip infinitely often , with their flipping times following power law distributions. The fast spins produce the persistence phenomenon.

In Fig. 3 we report with dashed lines the results of a microcanonical experiment. As before, we let the system relax for $`t_w`$ MCS and then we measure the persistence. The difference is that now the dynamics after time $`t_w`$ is performed keeping the energy constant to the value it had at time $`t_w`$. For $`tt_w`$ the data perfectly coincide with those for the usual persistence. Then the system leaves the configurational space region (state ) where it was confined by the dynamics and it goes far away (we check this also measuring correlations functions, that behave qualitatively like persistence). Thus we can conclude that the freezing phenomenon is completely dynamical . In other words the system seems to be stuck only because it always finds a path towards a lower energy configuration belonging to the same state before finding a path towards a different state.

For both dimensionalities we find that fast spins are organized in clusters. At any time and in each cluster there must be at least one “free” spin that can be flipped without any energy cost ($`\mathrm{\Delta }E0`$). The number of free spins equals the number of flat directions in the configuration space and decreases lowering the energy. Free spins usually appear in pairs, which wander inside their cluster .

Fast spins clusters sizes follow the pdf shown in Fig. 4. Measurements have been taken after a waiting time of $`t_w=10^4`$. Note that a finite $`t_w`$ can overestimate a little the clusters size, however for both dimensions we have that $`U_{\mathrm{}}(10^4)`$ is very near to $`U_{\mathrm{}}(\mathrm{})`$ and we did not find any dependence on $`t_w`$. Up to our accuracy, 3D data can be well fitted by a power law with an exponent $`1.94(3)`$. On the contrary 2D data clearly show a cut-off. We can not exclude that also 3D data would show a cut-off on larger scales.

For any finite cluster we can construct the set of allowed configurations and then we can evaluate the time a random walker needs to visit all of them. This time is, by construction, greater than the larger flipping time of any spin in the cluster. In this way we can link clusters sizes and flipping times, and one would expect to find a cut-off $`\tau _{max}`$ in the 2D distribution of flipping times as a consequence of the one present at $`c_{max}`$ in Fig. 4. However times may be exponential in the cluster size (like e.g. the recurrence time) and $`\tau _{max}`$ may be considered infinite for all practical purposes. Moreover, if one assume from percolation arguments that $`c_{max}`$ grows logarithmically with the system size $`N`$, then $`\tau _{max}`$ would grow as a power of $`N`$, like usually in any finite size systems.

In conclusion we have performed a numerical study of the low temperature dynamics of a finite-dimensional frustrated model with discrete couplings, namely the EA model in 2D and 3D, which can be viewed as a prototype model for glassy systems. The onset of a spontaneous time/space scale separation allows for a classification of spins in distinct classes characterized by their own dynamics as well as the identification of low-energy excitation responsible for structural rearrangements during the aging.

In the $`T=0`$ limit, while a part of the system completely freeze, the fast spins dynamics can be quantitatively described in terms of power law exponents which distinguish the 2D from the 3D case. The dynamics inside a cluster of fast spins, although diffusive-like, has strong constraints and can be interpreted along the lines of .

Given a proper definition of persistence in aging systems, we have been able to estimate the persistence exponent in the EA model, which turns out to be of kind $``$ , i.e. partially frozen and partially infinitely flipping. Finally, with a microcanonical experiment, we have shown that blocking phenomena are completely dynamical, even in disordered systems.

In spite of the central role played by space and time heterogeneities in glassy dynamics, such aspects have been rather poorly studied from the analytical point of view. Scope of this paper has been to provide some reference numerical result which hopefully could be captured by theoretical arguments. A first step in this direction can be done within the framework of finitely connected long range models (like e.g. the Viana-Bray one ). In these models the underlying graph generates non trivial heterogeneities, driven by the connectivity pattern.


text:
# Optimal packing of polydisperse hard-sphere fluids II

## I Introduction

Colloidal suspensions are never truly monodisperse, but are in general polydisperse. This might lead to unfavorable behavior, because it will for example not be possible to form a high quality crystal. Polydispersity will also have consequences on other properties like the viscosity of the system. A better understanding of these systems could therefore possibly lead to a handle to tune some of the relevant properties of colloidal suspensions.

In a recent paper  we considered a system of $`N`$ hard spheres of which the total volume is fixed. The spheres are, however, allowed to exchange volume under this constraint. Theory and Monte Carlo simulations showed that the Percus-Yevick (PY) theory, which normally is well suited for describing a polydisperse hard sphere system, has a breakdown at a low volume fraction ($`\eta 0.260`$). According to the theoretical description, the size-distribution function of the particles can not be normalized above this threshold. In simulations it is found that above this packing fraction a few of the particles will grow bigger on increasing pressure. They will contain most of the available volume and will be surrounded by a sea of small particles. The size of these big particles becomes of the order of the simulation box, therefore finite size effects do not allow to conclude how to describe the system at high pressures.

In this article we will focus on a similar system, but now the total surface of the particles is fixed. In the Monte Carlo simulations we performed, we allowed particles to exchange amounts of surface under the constraint that the total surface remains constant. An experimental system with this type of behavior, would be a system of surfactant molecules, which are free to aggregate into spherical micelles. For this model system we derive the Helmholtz free energy, which is given by the free energy of a polydisperse system with an arbitrary size-distribution, subject to two constraints, both the number of particles and the total surface of the particles are fixed. We define the optimal size-distribution of this system as the one that minimizes this restricted Helmholtz free energy.

The remainder of this paper is organized as follows: in Sec. II we derive the Helmholtz free energy of this system and calculate the optimal size-distribution. In Sec. III we show the results of our computer simulations and make a comparison with the theoretical predictions. In Sec. IV we discuss the possibilities to form a high density crystalline phases and conclude in Sec. V with a discussion of our results.

## II theory

In a multicomponent system the ideal entropy is given by a simple expression

$$Nk_B\underset{i}{}w_i\mathrm{ln}(\mathrm{\Lambda }_i^3\rho w_i),$$
(1)
where $`w_i`$ is the molar fraction of species $`i`$ and $`\mathrm{\Lambda }_i`$ its thermal wavelength. In a true polydisperse system, however, this entropy is infinite . Normally one would describe such a system by distributing the particles over boxes according to some property, e.g. diameter or volume, which would enable us to distinguish them. In the case of a fixed particle-size distribution this is perfectly allowed for whatever property one wishes to choose. In our case, however, the particle-size distribution is allowed to change and as a consequence the labeling property will influence the theoretical description of the system . The choice of how to label the different boxes is in fact determined by an a priori probability assumption, which in this case is dictated by the way polydispersity is sampled.

The way we sample the polydispersity can be thought of to describe a system with a large, constant number of tiny particles, which form $`N`$ spherical aggregates. Exchanging surface between the aggregates would then imply exchanging some of these tiny particles. The number of these particles forming the aggregate is proportional to the surface, which forms a natural labeling of our system.

The fact that the total surface of the $`N`$ spherical particles in conserved leads to a natural length scale $`\sigma ^2^{(1/2)}`$, where $`\sigma `$ is the diameter of a particles and $``$ denotes an average over particle-size distributions. We introduce the reduced pressure $`P^{}=(k_BT)^1P\sigma ^2^{(3/2)}`$, where $`k_BT`$, the Boltzmann constant times the temperature $`T`$, is our unit of energy.

For an $`n`$-component hard sphere mixture, the compressibility pressure in the Percus-Yevick approximation yields :

$$\frac{\pi }{6}P^{}=\frac{\xi _0}{1\xi _3}+\frac{3\xi _1\xi _2}{(1\xi _3)^2}+\frac{3\xi _2^3}{(1\xi _3)^3},$$
(2)
where the $`j`$-moment of the particle-size distribution $`\xi _j`$ is defined as

$$\xi _j=\frac{\pi }{6}\underset{i}{}\rho _i\sigma _i^j,$$
(3)
where the index $`i`$ is used to denote the different particle species, $`\rho _i=N_i/V`$ is their number density, and $`\sigma _i`$ is their diameter.

Equation (2) is also valid for a continuous size-distribution, in which case the sum in Eq. (3) is replaced by an integral. The corresponding expression for the chemical potential of a species with diameter $`\sigma `$ is

$$\mu ^{}=\mathrm{ln}\left[\rho \mathrm{\Lambda }^3W(s)\right]\mathrm{ln}(1\xi _3)+\frac{3\xi _2\sigma }{(1\xi _3)}+\frac{3\xi _1\sigma ^2}{(1\xi _3)}+\frac{9\xi _2^2\sigma ^2}{2(1\xi _3)^2}+\frac{\pi }{6}P^{}\sigma ^3,$$
(4)
where $`\mathrm{\Lambda }`$ is the de-Broglie thermal wavelength, $`\sqrt{h^2/(2\pi mkT)}`$, and $`W(s)`$ is the probability density to find a particle with a surface around $`s=\pi \sigma ^2`$. The pressure $`P^{}`$ is given by Eq. (2).

In an (NPT) description, the Gibbs free energy $`𝒢`$ of the system fulfilling the constraints, must be at a minimum, and is given by the following functional

$$𝒢[W(s)]=\mu ^{}W(s)𝑑s_0W(s)𝑑s_1sW(s)𝑑s.$$
(5)
The Lagrange multipliers $`_0`$ and $`_1`$ ensure that the conservation of the number of particles and of the total amount of surface respectively are satisfied. For a minimum of the Gibbs free energy the functional derivative of $`𝒢`$ with respect to $`W(s)`$ should be zero, which implies that $`W(s)`$ must be of the form:

$$W(s)=\mathrm{exp}\left\{\underset{i=0}{\overset{3}{}}\alpha _i\sigma ^i\right\},$$
(6)
where the values of $`\alpha _1`$ and $`\alpha _3`$ are explicitly given by

$`\alpha _1`$ $`=`$ $`{\displaystyle \frac{3\xi _2}{1\xi _3}},`$ (7)
$`\alpha _3`$ $`=`$ $`{\displaystyle \frac{\pi }{6}}P^{}.`$ (8)
The coefficients $`\alpha _0`$ and $`\alpha _2`$ are determined by the constraints on the number of particles and the total of available surface. Note that all $`\xi _i`$ $`(i=1,2,3)`$ are positive. Moreover, $`\xi _3`$ is equal to the volume fraction $`\eta `$, and is therefore necessarily less than one. Hence, $`\alpha _1`$ and $`\alpha _3`$ are always negative. This also implies that the particle-size distribution can always be normalized. This is different from our previous model , where the constraint on the total volume implied that $`\alpha _3`$ could be zero and gave rise to a critical endpoint on the equation of state, beyond which the PY-theory is not capable of describing that system. In this case, however, it is $`\alpha _2`$ which can change sign and become positive due to the constraint on the surface. This will lead to the formation of bimodal particle-size distributions in this model. An argument similar to the one used in Ref.  to show the break down of the Percus-Yevick equation of state, can also be used here to show that no break down will occur in this case, where the total surface of the particles is constrained .

By self-consistently solving Eqns. (7), (8) and the two equations imposed by the constraints, we obtain the equation of state of this model (Figs. 1 and 2) and the particle-size distributions in equilibrium according to the PY-theory. The equation of state is presented here both as function of the packing fraction $`\eta `$ and reduced density $`\rho ^{}=(N/V)\sigma ^2^{(3/2)}`$, because in contrast with most simulations the volume of the particles involved is not constant.

In the limit where the density goes to zero $`\alpha _1=\alpha _3=0`$ and the particle-size distribution function $`W(s)`$ is an exponential decaying function. For increasing densities the constraint on the total surface will lead to an increasing value of $`\alpha _2`$, which becomes zero at a packing fraction $`\eta =0.2169`$ and reduced pressure $`P^{}=0.7025`$. In order for the particle-size distribution function to become bimodal, however, it should have a local minimum. This requires its derivative with respect to the surface to be zero

$$\frac{dW(s)}{ds}=W(s)\left(\frac{\alpha _1}{2\sigma }+\alpha _2+\frac{3\alpha _3\sigma }{2}\right)=0.$$
(9)
This point, which is characterized by the equation $`\alpha _2^2=3\alpha _1\alpha _3`$, can be numerically evaluated and leads to a critical value of the density $`\eta _C=0.3889`$ and pressure $`P^{}=3.4019`$ above which bimodal behavior can be observed.

For the more accurate equation of state provided by Mansoori et al. , the pressure should be replaced by

$$\frac{\pi }{6}P_{HS}^{}=\frac{\xi _0}{1\xi _3}+\frac{3\xi _1\xi _2}{(1\xi _3)^2}+\frac{3\xi _2^3}{(1\xi _3)^3}\frac{\xi _3\xi _2^3}{(1\xi _3)^3}.$$
(10)
This equation of state is obtained by $`P_{HS}^{}=\frac{2}{3}P^{}+\frac{1}{3}P_{Vir}^{}`$, where $`P_{Vir}^{}`$ is the virial pressure. The analysis is identical to that of the PY equation of state, but leads to slightly different values. $`\alpha _2`$ becomes zero at $`\eta =0.2174`$ and $`P_{HS}^{}=0.7016`$, and the critical point is found at $`\eta _C=0.3944`$ and $`P_{HS,C}^{}=3.4362`$.

The different equations of state, lead here not only to a different pressure value for given density, but also to slightly different particle-size distributions. This is in contrast with the fixed volume model, where the different equations of state for given density would only influence the value of the pressure.

## III Simulations

For the Monte Carlo simulations we used the isothermal-isobaric ensemble or sometimes referred to as $`NPT`$-simulations . This means that the number of particles $`N`$, the pressure $`P`$ and the temperature $`T`$ are fixed. In the simulations there are three types of trial moves. The positions of the particles are allowed to change by small amounts and we allow the simulation box to fluctuate, in order to equilibrate with respect to the applied pressure. In the case of the simulations of crystals we use a rectangular box for which the three boxlengths are allowed to change independently, while for the other simulations the shape of the simulation box remains cubic.

A third type of trial move is required in order to sample the polydispersity. To this end we select at random two particles and attempt to exchange an amount of surface between them, which is uniformly drawn from the interval $`[\mathrm{\Delta }S_{max},\mathrm{\Delta }S_{max}]`$. Here $`\mathrm{\Delta }S_{max}`$ is the maximum amount of surface we allow to be exchanged, and is adjusted such that the acceptance of this move is between 35 and 50%. Note that although the total amount of surface of the particles is fixed the volume they will occupy will change.

The equation of state, resulting from our computer simulations of the isotropic fluid performed on 1000 particles, is shown in Figs. 1 and 2. It compares nicely with the equation of state of the Percus-Yevick approximation (2). The equation of state of Mansoori et al (10), performs even better.

In Fig. 3 we show several particle-size distributions obtained from our simulations of the isotropic fluid. For low pressure we find an exponential decay. For the reduced pressure $`P^{}=4.0`$ we observe that it is bimodal, which is in agreement with the theoretical predictions. For increasing pressure the fraction of small particles decreases. Theoretically it remains bimodal, since $`\alpha _1`$ in (6) is always negative, however it becomes very small and in our simulations has disappeared at $`P^{}=20.0`$. In Fig. 4, a comparison is made between the measured and predicted distributions. As can be seen they coincide nicely up to $`P^{}=20.0`$ and $`\eta =0.568`$, and upon increasing pressure the width of the unimodal distribution decreases.

This can be partially understood by the following argument. Upon increasing pressure the system will try to adapt by forming a higher density. This can be achieved by decreasing the volume of the simulation box. One way to accomplish this would be by crystallizing. In our simulation however there is another way in which this can be achieved: decreasing the volume occupied by the particles. This is possible, because the total surface of the particles is fixed, but their total volume is free to change. If a special case of the Power Mean inequality is applied to the particle diameters we obtain

$${}_{}{}^{3}\sqrt{\frac{1}{N}\sigma _i^3}\sqrt{\frac{1}{N}\sigma _i^2},$$
(11)
where equality only holds if all diameters are the same. Therefore the volume occupied by the particles can be lowered by forming a more uniform size-distribution.

If the pressure is increased even more the size-distribution starts to deviate again from the single peaked distribution (Fig. 3), and starts to develop second peaks. This suggest that the limit to what extend the system can equilibrate under higher pressures by becoming more uniform is reached. Although it is not observed within the duration of our simulations, one possibility is that it will try to stabilize by forming a crystalline structure with possible different sized particles.

## IV Crystal phases

As mentioned in the previous section this system can equilibrate in two ways to higher pressures, by forming a more uniform size-distribution or by forming a more ordered structure. In order to explore this last possibility we tried several possibilities.

For a monodisperse system the highest density which can be obtained, would be by forming a face-centered-cubic crystal (FCC) or the hexagonal close-packed crystal (HCP). The highest packing fraction $`\eta `$ which can be obtained by those crystals is $`\eta =\pi /\sqrt{18}=0.74048`$.

In our case however we can also explore the possibilities to form crystals with two or more particle sizes. One of these possibilities is the formation of a simple cubic crystal (SC). By alternating two particle sizes, as found in a NaCl-crystal, it can be easily shown that the highest packing fraction which this $`AB`$ type crystal can reach is $`\eta =0.79311`$. In this case the diameters in terms of the reference particle are given by 1.30656 and 0.541196. Note that for this calculation we have the additional constraint that the average particle surface is fixed.

We also considered the $`AB_2`$ structure as observed for an AlB2 crystal, where hexagonal closed packed layers of larger particles are separated by hexagonal rings of smaller particles. The highest packing fraction, which can be reached by this $`AB_2`$ type crystal, is $`\eta =0.78211`$ in the case that the particles have diameters 1.38828 and 0.73235.

We performed simulations starting with these three perfect crystals, using 720 particles for the FCC crystal, 1000 for the $`AB`$ crystal and 864 for the $`AB_2`$ crystal. For high pressures these crystals are all at least mechanically stable. The equations of state of these branches are also shown in Figs. 1 and 2. Note that due to the non-fixed occupied volume of the particles a crystal with higher packing fraction does not guarantee a higher number density. The particle-size distribution (Figs. 5, 6,7) shows the broadening of the peaks if the pressure is lowered, and a shift of the maximum of the later two. The lowest pressure points are close to the mechanical instabilities of these crystals and give an lower bound to the coexistence pressure, if any of these crystals is thermodynamically stable.

Compression on any of these crystal branches leaves the structure intact, and the same equation of state is followed. However the compression from the isotropic fluid does not show a transition to a crystal branch. This might be due to the length of our simulation, and the slow process of forming such a crystal from an unordered system.

Although the particle-size distribution functions for the high density crystals are almost monodisperse for the FCC crystal and bidisperse for the $`AB`$ and $`AB_2`$ structures, it is not possible to use this in order to simplify a cell theory to only one or two particle sizes . A comparison with simulation results of pure mono- and bidisperse systems with the “optimal” sizes, reveals a substantial difference in the equation of state (Fig. 8). Therefore polydispersity should be taken into account by incorporating the particle-size distribution function in the cell theory, which falls outside the scope of this article. As a consequence we can at this moment only conclude that these crystals are mechanically stable, but not which of them is or are thermodynamically stable.

We only considered the $`AB`$ and $`AB_2`$ compounds described above, because they can lead to packing fractions beyond that of a monodisperse FCC or HCP crystal. There are however more possible candidates, like the CsCl and CaF2 structures, or other compound types like $`AB_3`$, $`AB_4`$, $`AB_5`$ and $`AB_{13}`$ . Even more complex crystals could be obtained by using three or more particles sizes. A number of these structures are probably not stable, because either the the packing fraction can not be high enough, or the smaller particles are able to escape from their lattice positions .

## V Conclusion

We have investigated a polydisperse hard sphere system in which particles can exchange surface under the constraint that the total surface of the particles is fixed. In contrast with a similar model , the equation of state does not break down, and predictions made by the approximation of Mansoori et al works very well.

For low pressures the particle-size distribution is an exponential decaying function. On increasing pressure it becomes bimodal and eventually it forms a single peaked distribution by eliminating the smallest particles. If the pressure is increased even more there are signs that it might become bimodal again, e.g. by forming a crystal with two particle sizes. However, such a transition is not observed.

We performed simulations on three types of polydisperse crystals: a FCC-like crystal, an $`AB`$-crystal as found for NaCl, and an $`AB_2`$-crystal as found for AlB2. These crystals are all mechanically stable, and are just the most basic types of crystals possible. The polydispersity, however, influences the equation of state such, that a simple cell theory can not predict the true coexistence or thermodynamically stable crystal for these approximately mono- and bidisperse crystals.

The formation of an “Appolonian” packing, which in principle can lead to a space-filling structure, is not likely to happen, as the constraint on the total surface prevents this. The study of this and a previous model  have shown that these simple models lead to surprising results. Although they are both somewhat artificial we hope that they will stimulate the research towards more realistic, related models.

## Acknowledgments

We thank Daan Frenkel, José Cuesta and Emmanuel Trizac for a critical reading of the manuscript, as well as for their comments on and discussions about the results presented here.

FIGURE CAPTIONS

1. The equation of state of the polydisperse system. The circles are the simulation results of the isotropic fluid, the squares of an FCC crystal, the diamonds of an $`AB`$ crystal, and the triangles of an $`AB_2`$ crystal. The theoretical predictions, Percus-Yevick (dashed) and the Mansoori et al (solid), are also shown.
2. The equation of state of the polydisperse system as in Fig. 1, but now the reduced pressure $`P^{}`$ is given as function of the reduced density $`\rho ^{}`$. Note that the differences are caused by the fact that the occupied volume of the particles is not fixed.
3. The particle distribution function for several reduced pressures obtained by simulations of the isotropic fluid. At reduced pressure $`P^{}=4.0`$ the distribution is just bimodal. The highest pressure shown ($`P^{}=50.`$) suggest the formation of a second peak.
4. A comparison between the particle distribution function obtained from simulation (points) and the Mansoori et al approximation (lines), for reduced pressures $`P^{}=4.0`$, $`P^{}=20.`$, and $`P^{}=50.`$. There is a very good agreement, except for the highest pressure where the predicted curve remains nicely single-peaked, while the simulation suggest the formation of bimodal distribution.
5. The particle distribution function for several reduced pressures on the FCC branch. $`P^{}=22.`$ is the lowest pressure for which the crystal did not melt during the duration of our simulation.
6. The particle distribution function for several reduced pressures on the $`AB`$-crystal branch. $`P^{}=42.`$ is the lowest pressure for which the crystal did not melt during the duration of our simulation.
7. The particle distribution function for several reduced pressures on the $`AB_2`$-crystal branch. $`P^{}=32.`$ is the lowest pressure for which the crystal did not melt during the duration of our simulation.
8. A comparison of the equations of state of polydisperse and mono- or bidisperse systems. The squares correspond to a FCC crystal, the diamonds to an $`AB`$ crystal, and the triangles to an $`AB_2`$ crystal. The open symbols correspond to the polydisperse systems, while the filled symbols denote the results of true mono- or bidisperse systems with fixed “optimal” size ratios.


text:
# Superconductivity as a Bose-Einstein condensation?

## Abstract

Bose-Einstein condensation (BEC) in two dimensions (2D) (e.g., to describe the quasi-2D cuprates) is suggested as the possible mechanism widely believed to underlie superconductivity in general. A crucial role is played by nonzero center-of-mass momentum Cooper pairs (CPs) usually neglected in BCS theory. Also vital is the unique linear dispersion relation appropriate to weakly-coupled “bosonic” CPs moving in the Fermi sea—rather than in vacuum where the dispersion would be quadratic but only for very strong coupling, and for which BEC is known to be impossible in 2D.

Corresponding author:  M. de Llano dellano@servidor.unam.mx

$`^\text{a}`$Instituto de Física Teórica, Universidade Estadual Paulista, 01405-900 São Paulo, SP, Brazil
$`^\text{b}`$Deparment de Física, Universitat de les Illes Balears, 07071 Palma de Mallorca, Spain
$`^\text{c}`$Instituto de Física, Universidad Nacional Autónoma de México,
Apdo. Postal 20-364, 01000 México, DF, México
$`^\text{d}`$Instituto de Investigaciones en Materiales, Universidad Nacional Autónoma de México,
Apdo. Postal 70-360, 04510 México, DF, México
$`^\text{e}`$PESTIC, Secretaría Académica & CINVESTAV, IPN, 04430 México DF, México

Bose-Einstein condensation (BEC) of Cooper pairs (CPs) leads to a phase transition even in 2D in any many-fermion system dynamically capable of forming CPs. This transition could be germane to superconductivity in the quasi-2D cuprates. In the weak coupling limit one finds a nearly linear dispersion relation for the CP that suggests very high, even diverging, critical temperatures $`T_c`$. On the other hand, in the strong coupling limit a nearly quadratic dispersion relation gives vanishingly small $`T_c`$’s. For intermediate coupling one gets the finite $`T_c`$’s appropriate for real quasi-2D superconductors.

The single-Cooper pair problem may appear academic at first but has significant consequences. The familiar BEC formula for the transition temperature $`T_c3.31\mathrm{}^2n_B^{2/3}/m_Bk_B`$, with $`n_B`$ the number density of bosons of mass $`m_B`$ and $`k_B`$ the Boltzmann constant, is a special case of the more general expression  valid for any space dimensionality $`d>0`$ and any boson dispersion relation $`\epsilon _K=C_sK^s`$ with $`s>0`$ and $`C_s`$ a constant, given by

$$T_c=\frac{C_s}{k_B}\left[\frac{s\mathrm{\Gamma }(d/2)(2\pi )^dn_B}{2\pi ^{d/2}\mathrm{\Gamma }(d/s)g_{d/s}(1)}\right]^{s/d}.$$
(1)
If $`\mu _B(T)`$ is the boson chemical potential and $`e^{\mu _B(T)/k_BT}z`$ the fugacity, $`g_\sigma (z)_{l=1}^{\mathrm{}}z^l/l^\sigma `$. For $`z=1`$ and $`\sigma 1`$ this is just $`\zeta (\sigma )`$, the Riemann Zeta-function of order $`\sigma `$ which is finite for $`\sigma >1`$ and infinite for $`\sigma =1`$, while the series $`g_\sigma (1)`$ diverges for all $`\sigma 1`$. For $`s=2`$, $`C_2=\mathrm{}^2/2m_B`$, and since $`\zeta (3/2)2.612`$, this leads to the BEC $`T_c`$-formula cited above. Since $`g_{d/2}(1)`$ diverges for all $`d/21`$, $`T_c=0`$ for all $`d2`$. Eq. (1) follows from the boson number equation

$$N=N_0(T)+\underset{K0}{}\left[e^{\{\epsilon _K\mu _B(T)\}/k_BT}1\right]^1$$
(2)
where $`N_0(T)`$ is the number of bosons in the $`K=0`$ state. At $`T=T_c`$ both $`N_0(T_c)`$ and $`\mu _B(T_c)`$ virtually vanish so that using $`_k\left(L/2\pi \right)^dd^dk`$ in (2) then yields (1).

The fact that a CP can have a linear, as opposed to the usual quadratic, dispersion relation was mentioned as far back as 1964 by Schrieffer , p. 33. We have found it to be so for weak coupling; it becomes quadratic only for extremely strong coupling. In weak coupling $`\epsilon _K\mathrm{\Delta }_0\mathrm{\Delta }_Ka(d)\mathrm{}v_FK`$, where $`\mathrm{\Delta }_K`$ is the (positive) binding energy of a CP of center-of-mass momenta (CMM) $`\mathrm{}K`$, $`v_F\mathrm{}k_F/m`$ is the Fermi velocity, $`m`$ the fermion effective mass, while $`a(d)`$ $`2/\pi `$ and $`1/2`$ in 2D and 3D, respectively . For linear dispersion $`s=1`$, $`C_1=`$ $`a(d)\mathrm{}v_F`$, another special case of (1) is

$$T_c=\frac{a(d)\mathrm{}v_F}{k_B}\left[\frac{\pi ^{(d+1)/2}n_B}{\mathrm{\Gamma }[(d+1)/2]g_d(1)}\right]^{1/d}.$$
(3)
Now $`T_c=0`$ for all $`d1`$ only—and $`T_c>0`$ for all $`d>1`$, which is precisely the range of dimensionalities for all known superconductors if one includes the quasi-1D organo-metallic Bechgaard salts . Using the interpolation $`a(d)=(7/26/\pi )+(8/\pi 13/4)d+(3/42/\pi )d^2`$, which reduces to 1, $`2/\pi `$ and $`1/2`$ in 1D, 2D and 3D, respectively, Fig. 1 graphs (1) and (3) vs $`d`$ if one imagines all fermions in the many-fermion system paired into CPs of mass $`m_B=2m`$. The particle number density of the original fermions is $`nk_F^d/2^{d2}\pi ^{d/2}d\mathrm{\Gamma }(d/2)=2n_B`$. Here $`k_BT_FE_F=\mathrm{}^2k_F^2/2m`$. These curves are an upper bound to a more realistic $`T_c`$ where only a fraction of all fermions are actually bound.

The above holds for “unbreakable” CPs, i.e., bosons existing for all CMM so that the integrals over momenta go from 0 to $`\mathrm{}`$. However, CPs actually break up  into two “unpaired fermions” above a certain CMM $`K_0`$. This increases in value as coupling is increased and tends to infinity (as will be seen below) for infinite coupling. The CP breakup $`K_0`$, defined by $`\mathrm{\Delta }_{K_0}=0`$, for very weak coupling is just $`K_0\mathrm{\Delta }_0/a(d)\mathrm{}v_F0`$ as $`\mathrm{\Delta }_00`$. Eq. (2) with $`K_0`$ rather than $`\mathrm{}`$ as an upper limit eventually yields

$$T_c\frac{(d1)[2\pi a(d)\mathrm{}v_F]^dn_B}{A_dk_B\mathrm{\Delta }_0^{d1}}\underset{\mathrm{\Delta }_00}{}\mathrm{},$$
(4)
where $`A_d=2\pi `$ or $`4\pi /3`$ in 2D or 3D, respectively. Hence, since in that limit any remaining CPs are condensed in the $`K=0`$ state at all temperatures, the critical temperature $`T_c`$ must be infinite. Note that $`T_c0`$ explicitly for $`d=1`$. Indeed, if the unpaired fermions resulting from broken CPs are taken into account in a more realistic binary boson-fermion mixture model this artificially infinite $`T_c`$ is “tamed” down  to a finite $`T_c`$—which remarkably enough practically coincides with that of the pure unbreakable boson gas (3). Detailed calculations using the BCS model interaction gives a $`T_c`$ still almost an order of magnitude too large compared with the empirical range  of $`T_c(0.010.1)T_F`$ for cuprates—even for a moderate coupling of $`\lambda =1/2`$. However, the calculated $`T_c`$ of about 800 K is about 17 times larger than that of the BCS theory $`T_c1.13\mathrm{\Theta }_D\mathrm{exp}(1/\lambda )`$ of about 46 K. For cuprates $`d2.03`$ has been suggested  as more realistic as it accounts for inter-layer couplings, but results would be very similar to those discussed here for $`d=2`$.

A far more general interfermion interaction is the $`S`$-wave attractive separable potential  whose double Fourier transform is

$$V_{pq}=(v_0/L^2)g_pg_q.$$
(5)
Here $`L`$ is the size of the “box” confining the many-fermion system, $`v_00`$ is the interaction strength and $`g__p(1+p^2/p_0^2)^{1/2}`$ where $`p_0`$ is the inverse range of the potential. Hence, e.g., $`p_0\mathrm{}`$ implies $`g__p=1`$ and corresponds to the contact or delta potential $`V(r)=v_0\delta (𝐫)`$, and $`p_0=k_F`$ to a range of order of the average interfermion spacing, etc. Such an interaction model may mimic a wide variety of possible dynamical mechanisms in superconductors: a force mediated by phonons, or plasmons, or excitons, or magnons, etc., or even a purely electronic interaction. In the first instance mentioned one may have, e.g., a (possibly singular) coulombic interfermion repulsion surrounded by a longer-ranged electron-phonon attraction. The former is central to the problem at hand; the latter (whether or not phonon-based) is indispensable to create the CPs detected experimentally  in both conventional as well as cuprate superconductors. In vacuum, a two-body bound state with (positive) binding energy $`B_20`$ is possible for any 2D potential with sufficient attraction. The vacuum $`t`$-matrix  associated with the Lippmann-Schwinger equation then develops a pole at an energy $`E=B_2`$ so that for any dimension $`d`$ (5) leads to

$$\frac{L^d}{v_0}=\underset{k}{}\frac{g_k^2}{B_2+\mathrm{}^2k^2/m}.$$
(6)
In 1D the integral on the rhs is exact and gives the well-known binding energy $`B_2=mv_0^2/4\mathrm{}^2`$ of the single bound state, but in 2D that integral diverges logarithmically for large $`k`$. On the other hand, the CP equation for two fermions above the Fermi surface (of energy $`E_F`$) with momenta wavevectors $`𝐤_1`$ and $`𝐤_2`$ (and CMM wavevector $`𝐊𝐤_1+𝐤_2`$) is given by

$$\underset{k}{}{}_{}{}^{^{}}\frac{g_k^2}{\mathrm{}^2k^2/m+\mathrm{\Delta }_K2E_F+\mathrm{}^2K^2/4m}=\frac{L^d}{v_0}$$
(7)
where $`𝐤\frac{1}{2}(𝐤_1𝐤_2)`$ is the relative momentum wavevector, $`2E_F\mathrm{\Delta }_K`$ the total pair energy, $`\mathrm{\Delta }_K0`$ as before, and the prime on the summation implies restriction to states above the Fermi surface, viz., $`|𝐤\pm 𝐊/2|>k_F`$. In principle, (7) can be solved numerically for $`\mathrm{\Delta }_K`$ as a function of coupling strength $`v_0`$ for all $`0v_0<+\mathrm{}`$. More conveniently, however, for an interaction of the form (5) one can solve for $`\mathrm{\Delta }_0`$ (analytically if $`g_k=1`$) and for $`\mathrm{\Delta }_K`$ numerically, both as functions of $`B_2`$ alone, after eliminating the interaction strength $`v_0`$ by combining (6) and (7). This yields a renormalized CP equation employing $`0B_2/E_F<\mathrm{}`$  as dimensionless interaction coupling parameter instead of the (possibly singular) potential strength term $`1/v_0`$ is now eliminated. A much greater generality for the equation with $`g_k=1`$, viz., for any interfermion interaction with binding energy $`B_2`$, is conceivable. In 2D one has the surprising result (traceable to the fact that unlike 1D or 3D the fermionic density of states is a constant independent of $`k_F`$) that $`\mathrm{\Delta }_0=B_2`$ for all coupling, at least for an attractive delta interaction. Fig. 2 shows exact numerical results for the nonnegative CP excitation energy $`\epsilon _K\mathrm{\Delta }_0\mathrm{\Delta }_K`$ mentioned earlier, for zero range (full curves). For very weak coupling the exact curves are virtually linear, i.e., $`\epsilon _K2\mathrm{}v_FK/\pi `$, and for very large coupling tend aymptotically to the exact quadratic $`\mathrm{}^2K^2/2(2m)`$ (short-dashed curves). The long-dashed loops (reminiscent of the “roton” part of the excitation spectrum of liquid 4He) emerge for finite range at sufficiently strong coupling.

When the zero-range results like those in Fig. 2 are used in (2) assuming $`n_B=n/2`$ and $`m_B=2m`$, one gets Fig. 3 for the BEC $`T_c`$’s of a pure gas of breakable or unbreakable CPs. Consequently, for a substantial range of intermediate couplings in 2D finite BEC $`T_c`$’s are possible that lie within the empirical range  of cuprate $`T_c`$’s even for a simple pure boson gas of CPs. Significantly, $`T_c`$ is no longer zero as would be predicted in a BEC picture by a quadratic relation appropriate for “local-pair” CPs in vacuum—a result that convinced many that BEC is irrelevant for quasi-2D cuprate superconductors. More accurate BEC $`T_c`$’s should include refinements such as allowing for unpaired fermions  in a more realistic binary boson-fermion mixture model, CP-fermion interactions, non-$`S`$-wave interactions, etc.

Finally, CPs are here considered “bosonic” even though they do not obey (Ref.  p. 38) Bose commutation relations. This is because for a given $`K`$ they have indefinite occupation number as (in the thermodynamic limit) there are an indefinitely large number of allowed (relative wavevector) $`𝐤`$ values corresponding to an indefinitely large number of possible pairs of vectors $`𝐤_1`$ and $`𝐤_2`$. Hence, for any coupling and thus any degree of overlap between them, CPs do in fact obey the BE distribution (2) from which BEC is determined.

To conclude, the single CP problem with non-zero CMM (usually neglected in BCS theory) is illustrated in 2D for all coupling (and/or normal-state fermion density). BEC then suggests itself as a possible mechanism for superconductivity—in fact, for all observed superconductors from the quasi-1D organo-metallics to the more familiar 3D ones. But one must employ the correct dispersion relation in the CMM of the CPs. This relation is complicated but obtainable numerically for arbitrary coupling, and is purely linear for very weak and purely quadratic for very strong coupling.

Partial support from UNAM-DGAPA-PAPIIT (México) # IN102198, CONACyT (México) # 27828 E, DGES (Spain) # PB95-0492 and FAPESP (Brazil) is acknowledged.


text:
# Boundary induced phase transitions in driven lattice gases with meta-stable states

\[

## Abstract

We study the effect of meta-stability onto boundary induced phase transitions in a driven lattice gas. The phase diagram for open systems, parameterized by the input and output rates, consists of two regions corresponding to the free flow and jammed phase. Both have been entirely characterized. The microscopic states in the high density phase are shown to have an interesting striped structure, which undergoes a coarsening process, and survives in the thermodynamic limit.

\]

Driven lattice gases (DLG) are characterized by a non-vanishing mass-flow in the stationary state which is generated by non-potential forces . This feature of DLG has far reaching consequences, because even the stationary state of the system is not described in the framework of standard equilibrium statistical mechanics . In recent years the study of various non-equilibrium lattice gas models, which have important applications, e.g. in biological context  or as models for traffic flow , led to a profound understanding of some key features of systems far from equilibrium.

One of the most fascinating effects observed for this kind of models are so-called boundary induced phase transitions, i.e. the change of the bulk properties due to a variation of the boundary conditions . Such phase transitions have been extensively studied for a class of 1D DLG with a continuous flow density relation for periodic boundary conditions (PBC). For some models even exact results for the stationary state exist . A phenomenological theory  was recently able to predict the phase diagram in the case of open boundary conditions (OBC) also for more complicated models . However, less is known for DLG which exhibit meta-stable states for PBC , although they are relevant in different fields. Here, we study one of the simplest DLG model which shows meta-stability for PBC. We first discuss briefly the case of PBC, but our main focus is to study the effect of meta-stability in the case of OBC. Using a combination of analytical and numerical methods we explore the structure of the phase diagram and characterize the typical behavior in the different phases.

We study a DLG model defined on a chain of length $`L`$. Each site of the chain may either be empty ($`x`$) or occupied with a particle of velocity zero ($`0`$) or one ($`1`$). In a first step all particles which moved in the previous time step keep their velocity if the next site is empty. With probability $`q_0`$, velocity one is assigned to particles with velocity zero and an empty site ahead. In all other cases velocity zero is assigned to the particle. In the second step particles with velocity one move to the next site. Both steps are applied synchronously to all particles.

This model may be interpreted as a special case of the recently introduced VDR model for traffic flow , i.e. we restrict the maximal velocity of the particles to one. For $`q_0=1`$ we recover the deterministic asymmetric exclusion process with parallel dynamics .

Following , the structure of the flow-density relation, i.e. the fundamental diagram, in the case of PBC is a key to understand the behavior of the open chain. Fig. 1 shows the numerically established fundamental diagram for our model, in comparison with:

$`J(\rho )`$ $`=`$ $`\rho \text{ if }0\rho 0.5`$
$`J(\rho )`$ $`=`$ $`q_0(1\rho )\text{ if }q_0/(1+q_0)\rho 1`$
The form of the fundamental diagram can be easily understood. For densities $`0\rho 0.5`$ and suitable initial configurations (the initial velocity of the particles is one and no particle is blocked) all particles move deterministically with velocity one. At larger densities $`\rho >0.5`$ the number of particles exceeds the number of holes, i.e. one cannot avoid the formation of clusters of particles. The simulation results indicate that in the steady state one observes phase separation, i.e. the formation of a large compact cluster of particles which coexists with a free-flow regime. The density in the free-flow regime is determined by the frequency of particles leaving the cluster. The typical time the first particle needs to separate from the cluster is simply given by $`T=q_0^1`$. This time already determines the average gap between two particles in the low-density regime. Therefore the average density, $`\rho _{free}`$ in the free-flow regime is given by $`\rho _{free}=1/(T+1)=q_0/(1+q_0)`$. The number of moving particles and therefore the flux follows from particle conservation (see also ).

Both solutions for $`J(\rho )`$ coexist for $`\rho _{free}\rho 1/2`$. For finite systems it is possible that the large cluster dissolves due to fluctuations. In free flow states, fluctuations are absent and therefore no spontaneous cluster formation takes place at densities $`\rho 0.5`$. However, in the thermodynamic limit, any random initial configuration leads to a jammed stationnary state.

After this short discussion of the system with PBC, we proceed with the open chain. OBC are implemented as follows : at the left boundary particles with velocity one are introduced with probability $`\alpha `$ if the first site of the chain is empty. At the right end the particles leave the chain from site $`L`$ with probability $`\beta `$, irrespective of their velocity.

The values of $`\alpha `$ and $`\beta `$ determine for given values of $`L`$ and $`q_0`$ the bulk behavior of the system. First we notice that for $`\beta =1`$ the only stochastic element is due to the left boundary, i.e. we recover the deterministic asymmetric exclusion process. Therefore the exact result for the flow  is given by $`J(\alpha )=\alpha /(1+\alpha )`$. This result for the flow is expected to be valid not only for $`\beta =1`$ but in the whole low density phase, because for low bulk densities the flow is controlled by the input of particles.

Next, we consider the case of large bulk densities, i.e. we study the system’s performance for large values of $`\alpha `$ and $`\beta 1`$. For this case we expect to find frequently a compact cluster of particles at the right boundary. Therefore the mean time interval for a particle leaving the system is the sum of the average waiting time $`T_L=1/\beta `$ of the last particle at site $`L`$ and the average time which passes until the last site is occupied again, $`T_J=1/q_0`$. We find immediately that for large $`\alpha 1`$ and $`\beta 1`$ the flow is given by

$$J(\beta )=\frac{1}{T_L+T_J}=\frac{q_0\beta }{q_0+\beta }.$$
(1)
This simple scenario is only valid if the density of holes fed into the system is very low. For larger values of $`\beta `$ we expect that it is also possible that the cluster separates from the right boundary. In this case the density close to the exit is determined by the outflow from a cluster. The temporal headway $`t`$ between two consecutive particles is distributed via $`P_h(t)=q_0(1q_0)^{t1}`$ where $`t1`$. Now we calculate the typical amount of time which is needed from the arrival of a particle (at $`t=0`$) to the arrival of the second particle at ($`t=T`$) at site $`L`$. Here we have to distinguish two cases: $`(i)`$ the particle arrives at site $`L1`$ before the first particle has left the system or $`(ii)`$ that it arrives later, such that it reaches the last site without being blocked by the particle ahead. The probability that a particle gets blocked if it arrives at $`t`$ is given by $`P_b(t)=(1\beta )^t`$ and that it can pass without blockage by $`P_f(t)=_{\tau =1}^t\beta (1\beta )^{t1}`$.

This leads to the average value of $`T`$ given by

$`T`$ $`=`$ $`{\displaystyle \underset{t=1}{\overset{\mathrm{}}{}}}P_h(t)\left[P_f(t)(t+1)+P_b(t)(t+{\displaystyle \frac{1}{\beta }}+{\displaystyle \frac{1}{q_0}})\right]`$ (2)
$`=`$ $`{\displaystyle \frac{1}{\beta }}+{\displaystyle \frac{1}{q_0}}=J(\beta )^1.`$ (3)
Also for this scenario we recover the same value (1) of the flow. These two kind of configurations are generic if the output probability controls the capacity of the system, i.e. in the presence of clusters. Therefore we expect to find a unique functional dependence of the flow, if the system is controlled by the capacity of the exit. The transition between the $`\alpha `$ or $`\beta `$ controlled parameter regime is expected to be located at

$$\beta _t(\alpha )=\frac{q_0\alpha }{(1+\alpha )q_0\alpha }$$
(4)
where both fluxes agree.

Comparison of our prediction for the flow with simulation results shows excellent agreement (see fig. 2). Simulations for different lengths of the chain show that the differences between our predictions and the simulations systematically vanish for larger system sizes. The origin of these finite size effects, i.e. the enhancement of the flow for large $`\beta `$, is the possibility to find occasionally cluster-free micro-states. We indicated this by the broken lines in the phase diagram (Fig. 4). The absence of maximal current phase in the thermodynamic limit is related to the non-analyticity of the fundamental diagram .

We proceed characterizing the typical microscopic structure of the stationary state, in different parameter regimes. We distinguish three different domains, i.e. a low density domain (LD) for $`\beta >\beta _t(\alpha )`$ (LD), and two high density domains for $`\beta <\beta _t(\alpha )`$. Both high density domains correspond to a striped phase (SP), but they exhibit different behaviors in the vicinity of the entrance, depending whether $`\alpha <q_0`$ (SPI) or $`\alpha q_0`$ (SPII). In the LD phase we observe the typical two domain structure known from the asymmetric exclusion process. Clusters which eventually separate from the high density domain at the right dissolve in a few time steps. This dissolution leads to the localization of the high density domain at the right. The system takes a bulk density $`\rho _b(\alpha )=J(\alpha )`$. At site $`L`$ the density is given $`\rho _L=J(\alpha )/\beta `$. Fits of several density profiles follow the pure exponential form $`\rho (n)=\rho _b(\alpha )+[\rho _L\rho _b(\alpha )]\mathrm{exp}((Ln)/\xi )`$.

In the high density region (SPI + SPII), clusters are able to reach entrance. This region can be described by considering that at the right boundary, some free-flow segments are injected within a compact cluster. Typically, the free-flow segments reduce to isolated holes for $`\beta <q_0`$, whereas for $`\beta >q_0`$, they contain both empty sites and particles with velocity one. Both the free-flow segments and the compact clusters separating them move backwards, producing a striped structure in spatio-temporal plots (see fig. 3). As long as it is far enough from the left boundary, a cluster has the same in- and output rates, thus its width follows a non-biased random walk, whose average is a constant. However, there is a non-vanishing probability that the width of the cluster crosses the zero value, leading to the coalescence of the two neighboring free-flow segments, which will never be able to separate again afterwards. The excess particles are distributed among the neighboring clusters. Thus, as the stripes move towards the entrance, their widths increase. However, the average density remains a constant equal to the density selected on the left boundary $`\rho =J(\beta )/\beta `$. When the clusters arrive near the left boundary and are exposed to the entrance flux, their width suddenly obeys a biased random walk. If $`\alpha >q_0`$ (phase SPII), the width increases on average and the high input rate stabilizes the cluster until it reaches the boundary (fig. 3-b). If clusters and free-flow segments at a given position were monodisperse, then the sudden rate change would induce a linear density profile on the left, with a crossover to the constant bulk value occurring at site $`\stackrel{~}{n}=\alpha t_f/(1+\alpha )`$ where $`t_f`$ is the average time interval during which free-flow is observed in the first site, i.e. the typical time between the disappearance of a jam and the arrival of the next. Due to the coarsening of the striped pattern in upstream direction one observes a subextensive growth of $`t_f`$ with increasing system size. The density in the first site is determined by noticing that the flux (1) is also equal to $`j=\alpha (1\rho _1)`$. If we measure $`t_f`$ numerically, we find a quite good agreement between our estimate and the measured density profile (see fig. 5), except that the crossover region is in fact very large due to size dispersion of clusters and free-flow segments.

If $`\alpha <q_0`$ (phase SPI), the width of the leftmost cluster decreases on average due to the reduced capacity of the entrance. Some clusters will even not be able to reach the entrance, as seen on fig. 3. As long as we are far from the transition line $`\beta _t(\alpha )`$, still, clusters are on average large enough so that most of them reach the entrance and the depletion region remains localized on the left. The calculation done in SPII for the density profile is still valid. The only difference is that now the linear part has a positive slope. However, as we approach the transition line $`\beta _t(\alpha )`$, more and more clusters do not reach the boundary. When a given cluster disappears before reaching the left end, the next one is exposed earlier to the low incoming flux, and thus the depletion region may invade the whole system. The line separating the depletion region from the constant density region becomes delocalized over the whole system as the transition line is approached.

To conclude, we analyzed a DLG exhibiting meta-stable states. The structure of the phase diagram was explored by simulations as well as by a phenomenological description. The analytical predictions for each phase agree extremely well with the simulation results. The most characteristic feature is the spontaneous formation of simultaneously existing clusters for restricted outflow. The coarsening of these clusters as they flow backwards is due to the stochasticity of the model. We stress that this striped pattern survives in the thermodynamic limit. We expect that this behavior and more generally the structure of the phase diagram are generic for DLG with branched fundamental diagrams \- this branched structure being a signature of meta-stability. Moreover we established, in the limit of large system size, the absence of a maximal current phase.

Besides a phenomenological similarity with the clogging phenomena in granular flow , the practical relevance can be illustrated by the example of traffic flow, where the importance of boundary induced phase transitions is well accepted . Although we chose a rather simplistic model and do not aim to give a complete description of traffic flow, our results can be related to experimental studies that revealed the existence of parallel existing clusters on highways . These experimental results now can be directly related to the existence of meta-stable states. Moreover, our result shows that jams can be observed far upstream even when boundary induced. Another interesting feature are the important finite size effects, which probably can be used for a systematic flow optimization .

Acknowledgments: We thank Robert Barlovic, Andreas Schadschneider, and J. Krug for useful discussions. L. S. acknowledges support from the Deutsche Forschungsgemeinschaft under Grant No. SA864/1-1.