text:
# Valence-bond crystal in a {111} slice of the pyrochlore antiferromagnet

## Abstract

We investigate theoretically the ordering effect of quantum spin fluctuations in a Heisenberg antiferromagnet on a two-dimensional network of corner sharing tetrahedra. This network is obtained as a $`\{111\}`$ slice of the highly frustrated pyrochlore lattice, from which it inherits the equivalence of all three pairs of opposite bonds of each tetrahedron. The lowest-order (in $`1/S`$) quantum corrections partially lift the huge degeneracy of the classical ground state and select an ensemble of states with long-range valence-bond order.

Magnets with strong geometrical frustration have become a focus of experimental studies and inspired a growing number of theoretical inquiries into their unusual properties.Schiffer96etal  For classical spins, strong frustration precludes simple Néel order and creates a vast ground state degeneracy. This degeneracy, in turn, renders the magnet susceptible to nominally small perturbations (dipolar interactions, anisotropies, quantum effects etc.).

An antiferromagnet on the pyrochlore lattice—a network of corner-sharing tetrahedra—is perhaps the ultimate example of strong frustration.oldpyro ; Moessner98  Monte Carlo simulations with classical Heisenberg spins show that thermal fluctuations alone are ineffective in restoring magnetic order: the magnet remains in a paramagnetic but correlated state at a temperature as low as $`10^4JS^2`$. The effect of quantum fluctuations on the pyrochlore antiferromagnet has not yet been determined.

Previous studies have focused on the case of maximal quantum fluctuations, $`S=1/2`$. Early on, a ground state in the form of a valence-bond crystal (VBC) was proposed,Harris91  which preserves the spin-rotational $`O(3)`$ symmetry but breaks a number of discrete lattice symmetries, including the equivalence of the three pairs of opposite bonds of a tetrahedron. The relevant local order parameter is constructed from the expectation values $`𝐒_i𝐒_j`$. This result was obtained, in the absence (to this day) of a demonstrably reliable analytical treatment, via an uncontrolled approximation; several methods have since yielded similar results.Tsunetsugu01etal

An alternative route of studying quantum effects begins in the classical limit $`S\mathrm{}`$. At large $`S`$, quantum fluctuations are weak and can be evaluated in the framework of a perturbation theory, wherein physical quantities are expanded in powers of $`1/S`$. The large-$`S`$ approach has been previously applied to a variety of strongly frustrated magnets.Chalker91 ; Yildirim96 ; Sobral97 ; Tch03

In this paper, we report an application of the large-$`S`$ method to a $`\{111\}`$ slice of the pyrochlore lattice shown in Fig. 1(a). The resulting lattice can be described as a kagome plane flanked by two triangular lattices, so that every triangle is promoted to a tetrahedron. (Note that this is different from the pyrochlore slab geometry of SCGO,SCGO  where a triangular lattice is sandwiched between two kagome planes.)

From the symmetry point of view, our $`\{111\}`$ slice has the following features important in connection with the pyrochlore in 3D. All its units are equivalent tetrahedra and it thus permits collinear classical ground states (unlike SCGO). Further, unlike the checkerboard version of the 2D pyrochlore,Tch03  it retains equivalence between all three pairs of opposing bonds of each tetrahedron, a symmetry spontaneously broken in the $`S=1/2`$ pyrochlore VBC mentioned above.Harris91

For this lattice, we have found and characterized the subset of classical ground states favored by quantum fluctuations. The problem of minimizing the zero-point energy is equivalent to a discrete gauge-like theory, which by its very nature has many degenerate ground states. A natural first guess is that the ground states are those with uniform flux; this ensemble leads to a semiclassical valence bond liquid. However, our numerical minimisations reveal vacua with a pattern of gauge fluxes violating translational symmetry; this is reflected in a spatial modulation of nearest-neighbor spin correlations $`𝐒_i𝐒_i`$. Along with the absence of Néel order, this implies that the spin system is a valence-bond crystal.

We study the Heisenberg spin Hamiltonian

$$H=J\underset{ij}{}𝐒_i𝐒_j=(J/2)\underset{\alpha }{}|𝐋_\alpha |^2+\mathrm{const}=𝒪(S^2),$$
(1)
where $`J>0`$ and the sum is taken over all bonds $`(ij)`$ shown in Fig. 1(a). In the alternative representation, the Greek index $`\alpha `$ enumerates tetrahedra; $`𝐋_\alpha `$ is the sum of the four spins residing on tetrahedron $`\alpha `$. The energy (1) is minimized when $`𝐋_\alpha =0`$ for every tetrahedron $`\alpha `$. Classically, this condition does not determine a unique ground state. (This is the hallmark of frustration.) Rather, it defines an enormous manifold of classical vacua. By explicit construction, or Maxwellian counting,Moessner98  the dimension of the manifold can be seen to equal twice the number of tetrahedra. Thermal order by disorder is absent.

Quantum fluctuations produce virtual excitations out of the ground state, thereby lowering the vacuum energy. For large $`S`$, the quantum correction can be expanded in powers of $`1/S`$ by expressing the spin operators in terms of Holstein–Primakoff bosons. (The virtual excitations are magnon pairs produced by the operators $`b_i^{}b_j^{}`$ in the bosonized Hamiltonian.) The lowest-order correction is the zero-point energy of harmonic spin waves,

$$E^{(1)}=2NJS+\underset{a}{}\mathrm{}|\omega _a|/2=𝒪(S),$$
(2)
where $`\{\omega _a\}`$ is the set of spin-wave frequencies in a given classical ground state and $`N`$ is the number of tetrahedra. The first quantum correction (2) depends implicitly (through the magnon spectrum) on the spin configuration and thus can be viewed as an effective spin potential breaking the degeneracy of the ground-state manifold.

Typically, quantum fluctuations select collinear ground states,Henley87  provided that such vacua exist at the classical level. Although we are unaware of any general proof of this proposition, a compelling argument in its favor can be stated nonetheless. Given that the reduction of the vacuum energy is due to generation of virtual magnon pairs and that magnons are transverse excitations, the zero-point energy should be lower when more adjacent spins share both of their transverse directions, i.e. when spins are collinear. In certain cases this heuristic argument can be made precise.Yildirim96 ; Tch03 ; Chayes03  We assume that it applies in this case as well and minimize the zero-point energy in the set of collinear classical vacua. As a quick check, we have computed zero-point energies of a few simple classical vacua with the smallest magnetic unit cell, and found them to be $`1.2126JS`$, $`1.1612JS`$, and $`1.1436JS`$ for the collinear, “perpendicular,” and “anticollinear” states, respectively (notation of Ref. Henley87, ).

We parametrize collinear states by means of Ising variables $`s_i=\pm 1`$ such that $`𝐒_i/S=s_i\widehat{𝐧}`$, where $`\widehat{𝐧}`$ is a unit vector. The constraint of zero total spin on every tetrahedron carries over to the Ising variables: $`s_0+s_1+s_2+s_3=0`$. We can eliminate the variable $`s_0`$ on the out-of-plane site at the expense of introducing another constraint:

$$s_1+s_2+s_3=\pm 1\text{ on every triangle}.$$
(3)
It follows from Eq. (3) that collinear classical vacua of the pyrochlore slice can be mapped onto the ground states of a kagome Ising antiferromagnet. Redundancy of the out-of-plane Ising variables allows us to speak interchangeably of triangles and tetrahedra when we deal with ground-state configurations.

The zero-point energy (2) is an implicit function of the Ising variables $`\{s_i\}`$. As outlined in Refs. Moessner98,  and Tch03, , the spin-wave spectrum can be obtained by linearizing the classical equations of motion for the spins, written in terms of the $`𝐋`$ variables: $`\mathrm{}\dot{𝐋}_\alpha =_\beta 𝐒_{\alpha \beta }\times 𝐋_\beta .`$ Here $`𝐒_{\alpha \beta }`$ denotes the spin shared by tetrahedra $`\alpha `$ and $`\beta `$. (All of these spins reside in the kagome plane.) For a collinear spin state polarized along $`\widehat{𝐧}=(0,0,1)`$, transverse fluctuations of the total spin on tetrahedron $`\alpha `$ can be parametrized by a complex number $`\sigma _\alpha `$ as $`𝐋_\alpha =(\mathrm{Re}\sigma _\alpha ,\mathrm{Im}\sigma _\alpha ,0)`$. The magnon frequencies are then determined by solving the eigenvalue problem

$$JS\underset{\beta }{}s_{\alpha \beta }\sigma _\beta =\mathrm{}\omega \sigma _\alpha .$$
(4)
The zero-point energy can now be minimized numerically, e.g. by using simulated annealing. However, it turns out that the states minimizing Eq. (2) retain a large residual degeneracy. Hence, finding one such state does not yet solve the problem. We must first understand the nature of the residual degeneracy. This is best done by discussing a hidden symmetry responsible for it.

A gauge-like symmetry. Observe that the triangles $`\{\alpha \}`$ occupy the sites of a honeycomb lattice dual to the kagome \[Fig. 1(b)\]. The in-plane sites $`\{\alpha \beta \}`$ live on links of the dual lattice. It was noted by HenleyHenley  that Eq. (4) reveals a special gauge-like symmetry of the zero-point energy (2). Namely, a transformation

$$\sigma _\alpha \sigma _\alpha ^{}=\mathrm{\Lambda }_\alpha \sigma _\alpha ,s_{\alpha \beta }s_{\alpha \beta }^{}=\mathrm{\Lambda }_\alpha s_{\alpha \beta }\mathrm{\Lambda }_\beta ^1$$
(5)
leaves Eq. (4) invariant. (Conservation of Ising spin length $`|s_{\alpha \beta }|`$ requires that $`\mathrm{\Lambda }_\alpha =\pm 1`$.) Therefore, two spin configurations $`\{s_{\alpha \beta }\}`$ and $`\{s_{\alpha \beta }^{}\}`$ related by the gauge symmetry (5) have identical magnon spectra and, consequently, equal zero-point energies. We therefore expect a large residual degeneracy of the ground state even after the inclusion of the first quantum correction (2).

A word of caution. Transformation (5) is not exactly a gauge symmetry, even though it strongly resembles one. The transformed state must respect the ground-state constraint (3). The importance of this restriction is evident when one ponders the meaning of the transformation: if $`\mathrm{\Lambda }_\alpha =1`$, we flip all Ising spins on triangle $`\alpha `$. Doing so may violate the ground-state constraint (3) on adjacent triangles, where only one spin is flipped. This reduces the number of possible gauge transformations but—crucially—does not eliminate them entirely.

The uncovering of a hidden symmetry greatly simplifies our task: as soon as we find a single ground state minimizing the zero-point energy (2), we can generate all of the others by applying the gauge transformations (5). But the gauge-like symmetry is useful in more ways than one. It also puts restrictions on the way in which the zero-point energy depends on the Ising spin variables.

From the viewpoint of the $`Z_2`$ symmetry (5), the Ising spins $`s_{\alpha \beta }`$ are gauge variables. In a true gauge theory, the energy can only depend on their gauge-invariant combinations, in this instance the $`Z_2`$ fluxes through closed loops $`\mathrm{\Gamma }_{\mathrm{}}`$ of various lengths $`\mathrm{}`$ on the hexagonal lattice, $`\varphi _\gamma _{\mathrm{}}=s_{\alpha _1\alpha _2}s_{\alpha _2\alpha _3}\mathrm{}s_{\alpha _{\mathrm{}}\alpha _1}`$. The flux through any such loop can be factored into a product of fluxes piercing elementary plaquettes of the dual lattice, the hexagons $`\gamma =\mathrm{\Gamma }_6`$ \[Fig. 1(b)\]. Therefore, the zero-point energy will be a function of $`Z_2`$ fluxes:

$$E_1=E_1(\{\varphi _\gamma \}).$$
(6)
While ours is not precisely a gauge theory, Eq. (6) is the right guess: the zero-point energy is entirely determined by the values of the $`Z_2`$ fluxes. To prove this, we demonstrate that any two collinear ground states with the same set of flux variables $`\{\varphi _\gamma \}`$ have equal zero-point energies (2). To accomplish that, we show that the two states are related by a gauge-like transformation (5).

Proof. Let the Ising link variables be $`\{s_{\alpha \beta }\}`$ and $`\{s_{\alpha \beta }^{}\}`$ in the two states. Define new link variables as their ratio, $`t_{\alpha \beta }=s_{\alpha \beta }^{}/s_{\alpha \beta }=\pm 1`$. The variables $`t_{\alpha \beta }`$ are pure gauge: their product along any closed loop equals $`+1`$. Therefore we can unambiguously define a new variable $`\mathrm{\Lambda }_\alpha `$ on sites of the dual lattice in such a way that $`\mathrm{\Lambda }_\alpha =t_{\alpha \beta }\mathrm{\Lambda }_\beta `$ everywhere. Hence $`s_{\alpha \beta }^{}=\mathrm{\Lambda }_\alpha s_{\alpha \beta }\mathrm{\Lambda }_\beta ^1`$. By Eq. (5), the zero-point energies of the two states are equal. Q.E.D.

The precise form of the flux potential is not fixed by Eq. (6) and can be rather complicated. One might hope that a series expansion could successfully mimic it:

$$\frac{E_1}{N}=a_0+\frac{1}{N}\underset{\gamma }{}a_1\varphi _\gamma +\frac{1}{2N^2}\underset{\gamma ,\gamma ^{}}{}a_2(\gamma ,\gamma ^{})\varphi _\gamma \varphi _\gamma ^{}+\mathrm{}$$
(7)
If the linear term dominates, the ground state has uniform flux $`\varphi _\gamma =\mathrm{sign}(a_1)`$ on all elementary plaquettes. On the checkerboard, $`\varphi _\gamma =+1`$ in the ground state.Tch03

Extending the same simplistic approach to the pyrochlore slice we can take a guess that the ground states will either have no flux ($`\varphi _\gamma =+1`$) or have the maximal flux $`\pi `$ ($`\varphi _\gamma =1`$). While neither of these is the correct answer, they are interesting in their own right.

The $`\pi `$-flux states. In these ground states, the product of six Ising spins around any hexagon (of the original kagome lattice) equals $`1`$ \[Fig. 1(a)\]. The $`\pi `$ flux on the hexagon is related to the number of frustrated bonds. (A bond is frustrated if it connects Ising spins $`s_i=s_j`$.) More precisely, each hexagon has 0 or 2 frustrated bonds belonging to the triangles pointing down \[Fig. 1(b)\]. This fact allows us to map any $`\pi `$-flux state onto a dimer covering. If a hexagon has 2 such frustrated bonds, draw a dimer connecting the centers of the two triangles contributing the frustrated bonds \[Fig. 1(c)\]. Because every triangle has exactly one frustrated bond, it shares a dimer with one other triangle. The same construction applies to the sublattice of triangles pointing up \[Fig. 1(d)\]. We note that gauge transformations on the down triangles are sufficient to take any of the dimer states on the up triangles into any other without changing the dimer state on the down triangles – and vice versa. Thus every $`\pi `$-flux ground state maps onto a pair of completely independent dimer coverings of the triangular lattices.

Now one can use the well-known properties of the dimer model on the triangular latticeFendley  to characterize the ensemble of the $`\pi `$-flux states. By construction, the number of such states equals the number of pairs of dimer coverings. The latter grows exponentially with the number of sites, so that there is a finite entropy of $`0.4286`$ per tetrahedron. A disordered state of classical dimers on a triangular lattice translates into the absence of valence-bond order in our spin model. Spin correlations are also short-ranged. On a lattice with periodic boundary conditions (on a torus), these ground states fall into four disconnected topological classes distinguished by the values of the global $`Z_2`$ fluxes piercing the torus.

In short, the ensemble of $`\pi `$-flux states describes a classical valence-bond liquid.

The actual ground states. We have performed a numerical minimization of the zero-point energy (2) in the space of collinear vacua (3) using the method of simulated annealing. A Monte Carlo step, conserving the classical energy (1), consisted of flipping a randomly chosen chain of alternating spins. The zero-point energy (2) was then determined numerically.

Initial runs were done on systems containing between 100 and 200 tetrahedra. Starting from a periodic state we randomized the system by flipping randomly chosen antiferromagnetic chains until spin correlations with the original state were lost. The lattice was then annealed to achieve a state with the lowest energy. Throughout a simulation, we monitored Ising spins $`s_i`$, valence-bond variables $`s_is_j`$, and $`Z_2`$ fluxes on hexagons. The low-energy states obtained by the annealing showed no apparent spin or valence-bond order \[Fig. 2(a)\]. However, the $`Z_2`$ fluxes showed a clear tendency to order: 1/4 quarter of the hexagons had the trivial $`Z_2`$ flux $`+1`$ and formed a regular pattern \[Fig. 2(a)\]. We denote such states as $`2\times 2`$ after the size of the unit cell.

Finite-size scaling studies (with up to $`10^4`$ tetrahedra) yielded the zero-point energies of $`1.2127JS`$, $`1.2462JS`$, and $`1.2497JS`$ per tetrahedron for the no-flux, $`\pi `$-flux, and $`2\times 2`$ states, respectively, with the $`2\times 2`$ states the lowest. On the basis of the energy differences, we estimate that quantum selection of ground states will be effective below the temperature scale of roughly $`10^2JS/k_B`$. Note that the naïve energy scale $`JS`$ comes with a numerically small prefactor.

To obtain a dimer model appropriate for the $`2\times 2`$ flux pattern, we start with the dimers on a triangular lattice, which yield the $`\pi `$-flux states \[Fig. 1(c)\]. To change the flux from $`1`$ to $`+1`$ on one quarter of the hexagons, we decorate the triangular lattice by modifying one quarter of triangles \[Fig. 2(a)\]. Again, every frustrated bond is intersected by a dimer.

The properties of classical dimers on the decorated lattice can be readily computed by using Grassmann variables. The number of $`2\times 2`$ ground states is exponentially large, with an entropy per tetrahedron of $`0.4248`$. The probability to find a dimer on a short link is $`p=0.4004>1/3`$. In spin language this means that $`s_is_j=p(+1)+(1p)(1)=0.1992>1/3`$ on links surrounding the hexagons marked $`+1`$. These bonds are therefore more frustrated than the lattice average. On the other hand, a tetrahedron placed symmetrically with respect to the three closest fluxes $`+1`$ has $`s_is_j=1/3`$ for all its bonds. The valence-bond averages $`s_is_j`$ are shown schematically in Fig. 2(b). They are nonuniform, therefore the ensemble of $`2\times 2`$ states has valence-bond order that enlarges the lattice unit cell. (We note in passing a similarity with the valence-bond order proposed earlier for the three-dimensional $`S=1/2`$ pyrochlore:Harris91 ; Tsunetsugu01etal  3/4 of tetrahedra have unequal valence-bond correlations; the other 1/4 remain symmetric.)

Spin correlations $`s_is_j`$ between more distant neighbors decay exponentially with distance. The correlation length is 1.15 (the unit of length is the side of a tetrahedron). This points to the absence of Néel order.

To summarize, we find that, to leading order in $`1/S`$, quantum fluctuations select a collinear ensemble of degenerate ground states with additional bond order and no spin order: we have a valence-bond crystal. This degeneracy may be lifted in higher orders in $`1/S`$, as found previously for the Heisenberg antiferromagnet on kagome.Chubukov92 ; Henley95

The authors thank A.G. Abanov, C. Broholm, C.L. Henley, S.L. Sondhi, and O.A. Starykh for useful discussions. R. M. was supported in part by the Ministère de la Recherche et des Nouvelles Technologies.


format: text

text:
# The formation of the double neutron star pulsar J0737–3039

## 1 Introduction

A new double neutron star (DNS) system, J0737–3039, was recently discovered (Burgay et al. 2003). This is the fifth galactic DNS; all share the same characteristics of having eccentric orbits ($`e0.08`$) and short pulse periods (22 – 104 ms) and low period derivatives ($`10^{18}\mathrm{s}\mathrm{s}^1`$). The two latter parameters imply a weak surface magnetic field ($`10^{10}`$ G) which indicates that the pulsar has been recycled (Smarr & Blandford 1976; Alpar et al. 1982; Radhakrishnan & Srinivasan 1982), and that the observed pulsar in each system is the older neutron star and has gone through a mass-transfer phase (Srinivasan & van den Heuvel 1982).

Dewi & Pols (2003) found that the wide-orbit DNSs J1518+4904 and J1811–1736 are produced from helium star-neutron star binaries that avoided Roche-lobe overflow (RLOF); while the close-orbit B1913+16 and B1534+12 can be formed either with or without a mass-transfer phase in the helium star-neutron star binary stage of evolution. Among the five galactic DNSs, J0737–3039 has the shortest spin period (22 ms), shortest orbital period ($`0\stackrel{d}{.}102`$), and lowest eccentricity (0.088). In this paper we examine the possible formation history of J0737–3039. We will show that unlike the other DNSs, J0737–3039 can only be formed from helium star-neutron star binaries which have gone through RLOF.

## 2 The formation of J0737–3039

### 2.1 Possible progenitors

We examine the formation history of the system is terms of the so-called standard scenario (Bhattacharya & van den Heuvel 1991), in which DNSs were formed from Be/X-ray binaries. When the Be star leaves the main sequence and becomes a giant, a runaway mass transfer takes place due to the large mass ratio of the system. The neutron star is engulfed by the giant’s envelope, and a common-envelope (CE) phase ensues in which the neutron star spirals down to the helium core of the giant, ejecting its hydrogen-rich envelope in this process. The outcome of this spiral-in is a very close binary consisting of a helium star and a neutron star (HeS-NS). For companions less massive than about 5 $`\mathrm{M}_{\mathrm{}}`$, this mass transfer can still be stabilized such that the formation of a CE can be avoided (e.g. King & Ritter 1999; Podsiadlowski & Rappaport 2000; Tauris, van den Heuvel, Savonije 2000). However, for more massive companions such as the Be stars considered here the formation of a CE cannot be avoided.

The helium star evolves further through core-helium burning and further burning stages. During this evolution it transfers mass to the neutron star by RLOF and also loses mass by stellar wind. Fig. 1 depicts the parameters of these HeS-NS binaries at the end of this RLOF and wind mass loss, resulting from detailed evolutionary calculations (Dewi et al. 2002; Dewi & Pols 2003). The number indicated at each curve represents the initial mass of the helium star at the beginning of its evolution. The stars indicate the outcome of the evolution of systems with a range of initial orbital periods (note that the lines connecting the stars are not evolutionary tracks). Also presented in the figure are the maximum and minimum possible orbital periods of the progenitor of J0737–3039 just prior to the supernova (SN) explosion. These lines are derived with the assumption that the pre-SN (circular) orbital separation, $`a_\mathrm{t}`$, lies between the periastron and apastron of the post-SN separation, i.e.

$`a_\mathrm{f}(1e_\mathrm{f})a_\mathrm{t}a_\mathrm{f}(1+e_\mathrm{f})`$ (1)
(Flannery & van den Heuvel 1975), where $`a_\mathrm{f}`$ and $`e_\mathrm{f}`$ are the post-SN separation and eccentricity. Here $`a_\mathrm{f}`$ and $`e_\mathrm{f}`$ are the parameters at the time when it was born, which is assumed to be a characteristic time ago, $`\tau =P/2\dot{P}`$. These parameters were derived by integrating the effects of gravitational-wave radiation on the observed ($`a`$, $`e`$) over $`\tau `$.

Dewi & Pols (2003) found that systems represented with the open star symbols undergo a common-envelope (CE) phase due to the development of a helium convective envelope. The consequence of this is that, after the CE phase these systems will have a lower mass ($`M_{\mathrm{He},\mathrm{t}}1.4`$ $`\mathrm{M}_{\mathrm{}}`$) and much shorter orbital period ($`P_\mathrm{t}0\stackrel{d}{.}01`$) than their original parameters shown in Fig. 1. Hence, only the helium stars represented by solid star symbols (i.e. with initial masses in the range of 4 – 6.5 $`\mathrm{M}_{\mathrm{}}`$) are possible progenitors of J0737–3039. This range of mass requires an initial progenitor main-sequence star of 14 – 20 $`\mathrm{M}_{\mathrm{}}`$.

The region of the possible pre-SN orbital period $`P_\mathrm{o}`$ limited by the two thick lines in Fig. 1 does not cross the shaded area where the helium stars never fill their Roche lobes. This implies that J0737–3039 can only have been formed through Roche-lobe mass transfer from the helium star progenitor – as opposed to the other four galactic DNSs which could be formed either via RLOF or without a mass-transfer phase (Dewi & Pols 2003).

Fig. 2 depicts an enlargement of the region of allowed final HeS-NS systems from Fig. 1 in an $`a_\mathrm{t}`$ vs. $`M_{\mathrm{He},\mathrm{t}}`$ diagram (the stars are the same as the allowed systems in Fig. 1). For a given combination of pre-SN helium star mass $`M_{\mathrm{He},\mathrm{t}}`$ and orbital separation $`a_\mathrm{t}`$ we calculated the kick velocity required to produce the post-SN orbital period and eccentricity. Fig. 2 shows the ranges of the required kick velocities $`V_\mathrm{k}`$ limited by the curves of $`V_\mathrm{k}=`$ 70, 125, 175 and 230 $`\mathrm{km}\mathrm{s}^1`$. The region where the systems undergo a CE phase and evolve to extremely short orbital periods is to the left of the dash-dotted line in Fig. 2. Assuming that only systems which do not go through the CE phase are the possible progenitors of J0737–3039 (i.e. the systems on the right-hand side of the line), we find that to explain the orbital parameters of J0737–3039 when it was born (i.e. $`a=1.688`$ $`\mathrm{R}_{\mathrm{}}`$ and $`e=0.138`$; assuming that the pulsar’s age equals its characteristic age), a minimum kick velocity of 70 $`\mathrm{km}\mathrm{s}^1`$ is required at the birth of the pulsar’s companion.

The minimum separation, below which mass-transfer phase from the helium star to a neutron star is dynamically unstable resulting in a merger of the two stars (Dewi et al. 2002), is marked by the dotted line in Fig. 2. This line gives a constraint on the maximum kick velocity, i.e. 230 $`\mathrm{km}\mathrm{s}^1`$.

### 2.2 Magnetic field and spin history of the pulsar

All five double neutron star systems (DNSs) in eccentric orbits in the galactic disk share the characteristic of a combination of a weak dipole magnetic field (in the range $`10^9`$ to $`2\times 10^{10}`$ G) and rapid spin, indicating that we are dealing here with recycled pulsars that have been spun up during an accretion phase (Smarr & Blandford 1976; Alpar et al. 1982; Radhakrishnan & Srinivasan 1982; Srinivasan & van den Heuvel 1982). During such a phase with a non-degenerate donor companion the orbit will have been tidally circularized. Therefore the present eccentricity of these systems can only have been induced in the second supernova explosion when the companion collapsed to a neutron star (Srinivasan & van den Heuvel 1982). The reduction of the magnetic field strength of the pulsar to a value around $`10^9`$ to $`10^{10}`$ G is thought to somehow be the result of the accretion, as there is an observed correlation between the amount of matter accreted and the field strength (Taam & van den Heuvel 1986). A simple relation for this case was derived semi-empirically, by Shibazaki et al. (1989):

$`B(t)={\displaystyle \frac{B_0}{1+\mathrm{\Delta }M/m_\mathrm{B}}}`$ (2)
where $`B(t)`$ and $`B_0`$ are the present and initial magnetic field strengths, respectively, $`\mathrm{\Delta }M`$ is the amount of matter accreted and $`m_\mathrm{B}`$ is a constant to be derived from the observations. Francischelli et al. (2002) found $`m_\mathrm{B}=`$ $`1.25\times 10^5`$ $`\mathrm{M}_{\mathrm{}}`$ to agree well with the observed data for Low Mass X-ray Binaries. Adopting this value and a duration of the Roche-lobe overflow phase of the helium star of $`2\times 10^4`$ yrs (this is about the average value for helium stars in the mass range 4.0 to 6.3 $`\mathrm{M}_{\mathrm{}}`$ (Dewi et al. 2002; Dewi & Pols 2003)), with accretion at an Eddington rate for helium ($`3\times 10^8`$ $`\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$) one finds that $`B(t)`$ can have decreased by a factor of about 50 during this phase alone. Furthermore, a 15 $`\mathrm{M}_{\mathrm{}}`$ Be star resides $`10^7`$ yrs on the main-sequence. Be/X-ray binaries tend to be transients (Rappaport & van den Heuvel 1982; van den Heuvel & Rappaport 1987; Ziółkowski 2002) and are ”on” as an X-ray source only part of the time, during which they often reach a luminosity close to the Eddington limit. Assuming them to be ”on” only one per cent of the time and then to reach the Eddington limit, the neutron star will during the Be/X-ray binary phase have accreted about $`1.5\times 10^3`$ $`\mathrm{M}_{\mathrm{}}`$. With eq. (2) this is already sufficient to reduce its magnetic field by a factor of about 100. So, adding the amount of matter accreted during the Be/X-ray binary phase and the RLOF of the helium star, one finds with eq. (2) that these phases together provide enough accretion to reduce the magnetic field by a factor of about 200.

Magnetic fields of young pulsars on average typically have a dipole strength a few times $`10^{12}`$ G, and after reduction by this factor will typically be around $`10^{10}`$ G, as observed in the five DNSs. As to the spin-up of the ”old” neutron star: during the Roche-lobe overflow from the helium star accretion to the neutron star will have proceeded through a disk. The amount of angular momentum accreted during this phase can be calculated as follows. The mass transfer rate during Roche-lobe overflow is highly super-Eddington, in the range $`12\times 10^4`$ $`\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$ for helium stars in the mass range 4.0 – 6.3 $`\mathrm{M}_{\mathrm{}}`$ (Dewi et al. 2002; Dewi & Pols 2003). Dewi et al. (2002) and Dewi & Pols (2003) in their calculations of the evolution of the systems during this phase assumed the accretion rate onto the neutron star to be the Eddington rate ($`3\times 10^8`$ $`\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$) and the super-Eddington fraction of the transferred matter to be blown out of the system by the radiation pressure, carrying with it the specific angular momentum equal to the specific orbital angular momentum of the neutron star. Following this assumption of accretion at the Eddington rate – as radiation pressure will prevent accretion at a higher rate – and using the dipole strength of the magnetic field of J0737–3039 of $`6.9\times 10^9`$ G, one can calulate the value of the Alfven-radius $`R_\mathrm{a}`$, where the matter becomes coupled to the magnetic field and starts to co-rotate with the neutron star magnetosphere. $`R_\mathrm{a}`$ is given by (e.g. see Bhattacharya & van den Heuvel 1991)):

$`R_\mathrm{a}={\displaystyle \frac{(B_9R_6^3)^{4/7}}{(\dot{M}/\dot{M}_{\mathrm{E},\mathrm{H}})^{2/7}(M/\mathrm{M}_{\mathrm{}})^{1/7}}}\mathrm{.32.9}\mathrm{km}`$ (3)
where $`B_9`$ is the dipole field strength in units of $`10^9`$ G, $`R_6`$ is the neutron star radius in units of $`10^6`$ cm, $`\dot{M}`$ and $`\dot{M}_{\mathrm{E},\mathrm{H}}`$ are the accretion rate and the Eddington accretion rate for hydrogen, respectively and $`M`$ is the mass of the neutron star. With the dipole field strength of $`6.9\times 10^9`$ G of PSR J0737–3039, the Eddington rate for helium accretion $`\dot{M}=\dot{M}_{\mathrm{E},\mathrm{He}}=`$ $`3\times 10^8`$ $`\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$ and a $`M`$ = 1.3 $`\mathrm{M}_{\mathrm{}}`$ one obtains $`R_\mathrm{a}`$ = 78.4 km.

The rate at which angular momentum is accreted by the neutron star is then:

$`\dot{J}=\dot{M}\mathrm{\Omega }_\mathrm{k}(r=R_\mathrm{a})R_\mathrm{a}^2`$ (4)
where $`\mathrm{\Omega }_\mathrm{k}(r)`$ is the keplerian angular velocity of matter in the disk at distance $`r`$ from the center of the neutron star. With $`\dot{M}=\dot{M}_{\mathrm{E},\mathrm{He}}`$, $`M=1.3`$ $`\mathrm{M}_{\mathrm{}}`$ one finds $`\dot{J}=1.1\times 10^1\mathrm{M}_{\mathrm{}}\mathrm{km}^2\mathrm{s}^1\mathrm{yr}^1`$. The present spin angular momentum of the neutron star is

$`J=M(2\pi /P)R^2r_\mathrm{g}^2`$ (5)
where $`R`$ and $`P`$ are the radius and spin period of the neutron star and $`r_\mathrm{g}`$ is its radius of gyration, respectively. The value of $`r_\mathrm{g}^2`$ is typically about 0.075 for a polytropic index $`n=3`$ (e.g. Lattimer & Prakash 2001). With this value and $`R=10`$ km and $`P`$ = 22 ms one obtains $`J=2.8\times 10^3\mathrm{M}_{\mathrm{}}\mathrm{km}^2\mathrm{s}^1`$, which together with the above mentioned rate of accretion of angular momentum results in a required timescale for the spin-up of about $`2.5\times 10^4`$ years.

The duration of the phase of Roche-lobe overflow in the helium star plus neutron star binaries with helium star masses between 4.0 and 6.3 $`\mathrm{M}_{\mathrm{}}`$ ranges from 2.4 to $`0.7\times 10^4`$ years, which is of the same order as the value of $`2.5\times 10^4`$ years required for the spin-up. It thus appears that the magnetic field strength and spin rate of PSR J0737–3039 fit quite well with evolutionary model presented here in which the immediate progenitor system was a very close helium star plus neutron star binary which evolved through a phase of Roche lobe overflow with a duration of order $`10^4`$ years. As mentioned above, the progenitor of this system was so close that it cannot have avoided the phase of mass transfer by Roche-lobe overflow from the helium star (Dewi & Pols 2003), so here this phase must absolutely have taken place in the history of PSR J0737–3039. This may be the reason why the pulsar has the shortest spin period of all. The progenitors of the other systems may, in principle, have avoided RLOF, although it can certainly not be excluded that also in (some of) them Roche-lobe overflow from the helium star may have taken place.

## 3 Discussions and Conclusions

We have examined the formation of J0737–3039. We find that this system must have gone through a mass-transfer phase from a helium star with the initial mass in the range of 4 – 6.3 $`\mathrm{M}_{\mathrm{}}`$. We also find that an asymmetric SN explosion with a kick velocity of 70 – 230 $`\mathrm{km}\mathrm{s}^1`$ is necessary to explain the observed orbital parameters. However, these constraints may still be subject to some uncertainties. Whether the neutron stars in the systems that go through a CE phase (the open stars in Fig. 1) complete the spiral-in in the envelope of the helium star depends on the competition between the spiral-in timescale and the time left before the explosion of the helium star (Dewi & Pols 2003). If the SN explosion in systems with open star symbols in Fig. 1 occurs before the neutron star completes the spiral-in process, a lower kick velocity or even a symmetric SN explosion (shown as the solid circle) would then be possible. Furthermore, the constraints on the maximum pre-SN mass and the maximum kick velocity depend somewhat on the assumed wind mass-loss rate and the loss of angular momentum during the RLOF. As the stability of a mass-transfer phase depends on the mass ratio of the system and the period at the onset of RLOF, applying a higher wind mass-loss rate results in a lower mass helium star in wider orbit, and hence, more stable mass transfer than adopting wind mass loss with a lower rate.

## Acknowledgements

This work was sponsored by NWO Spinoza Grant 08-0 to E. P. J. van den Heuvel. We thank Dick Manchester for releasing the observational data of the system prior to publication, Onno Pols and Andrea Possenti for suggestions and discussions, and Hans Ritter for his referee comments.


text:
# Infrared behaviour and fixed points in Landau gauge QCD

## Abstract

We investigate the infrared behaviour of gluon and ghost propagators in Landau gauge QCD by means of an exact renormalisation group equation. We explain how, in general, the infrared momentum structure of Green functions can be extracted within this approach. An optimisation procedure is devised to remove residual regulator dependences. In Landau gauge QCD this framework is used to determine the infrared leading terms of the propagators. The results support the Kugo-Ojima confinement scenario. Possible extensions are discussed.

preprint: FAU-TP3-03-12, CERN-TH-2003-300, SHEP-03-33

Charged particles have long-range nature, whereas localised objects are neutral. For local quantum field theories one derives more generally that every gauge-invariant localised state is singlet with respect to the unbroken global charges of the gauge symmetry. Thus, in QCD without Higgs mechanism any localised physical state must be colourless. This extends to all physical states only if a mass gap is present. Then, colour-electric charge superselection sectors cannot arise, and every gauge-invariant state must also be a colour singlet. These are the signatures of confinement.

Within covariant linear gauges, the necessary conditions for confinement were formulated more than twenty years ago by Kugo and Ojima Kug79 . They express the necessity of a mass gap and the conditions for avoiding the Higgs mechanism in a covariant continuum formulation of QCD in terms of local fields, i.e., for the local measurement of (colour) quantum numbers. Both conditions are encoded in the critical infrared exponents of gluonic and ghost correlations in Landau gauge QCD Alkofer:2000wg . For the mass gap, the massless transverse gluon states of perturbation theory must be screened nonperturbatively, leading to an infrared suppression for the gluon propagator. In turn, the global gauge charges entail the infrared enhancement of the ghost correlations leading to their infrared dominance.

The infrared behaviour just described has been found in solutions of truncated Dyson-Schwinger equations Sme97 . Since then, the same qualitative behaviour has emerged in a variety of nonperturbative approaches including an increasing number of lattice simulations Cucchieri:1997fy  and investigations based on Fokker-Planck type diffusion equation of stochastic quantisation Zwa02 . Quantitatively, the latter leads to the same infrared exponents as the Dyson-Schwinger equations, which reflects an equivalence between the two formulations under certain conditions Zwa02 ; Lerche:2002ep . Interesting in their own right, such interrelations are especially useful for finding ways to go beyond the particular approximations employed in each case.

Here we present an investigation of the infrared sector of Landau gauge QCD by means of a nonperturbative flow equation Wegner ; Wetterich:yh . Central to the approach is the effective action $`\mathrm{\Gamma }_k`$, where quantum fluctuations with momenta $`p^2>k^2`$ are already integrated out, and $`k`$ denotes an infrared cut-off scale. The variation with $`k`$ leads to a flow equation for $`\mathrm{\Gamma }_k`$, thereby interpolating between the classical action in the ultraviolet and the full quantum effective action in the infrared where the cut-off is removed. Moreover, it has been shown that truncations can be systematically improved by means of an optimisation procedure Litim:2001up . The formalism has its particular merits in gauge theories as it is amiable to approximations being non-local in momenta and fields, precisely the behaviour we expect in QCD. So far it has been applied to Landau gauge QCD for a determination of the heavy quark effective potential Ellwanger:1996qf  and effective quark interactions above the confinement scale Bergerhoff:1997cv ; for a brief review and further applications in Yang-Mills theories see Litim:1998nf  and Pawlowski:1996ch  respectively.

Conceptually, the flow comprises a system of coupled differential equations for all Green functions which only involves dressed vertices and propagators. Hence, the correct RG-scaling is displayed by the full flow and truncations with the correct symmetry properties. This is one of the main contradistinctions to Dyson-Schwinger equations or stochastic quantisation which involve both bare and dressed Green functions. This makes it often difficult in the latter two approaches to implement the correct RG-scaling behaviour within truncations.

In this letter we explain how the full infrared momentum behaviour of Green functions is extracted within a fixed point regime. In Landau gauge QCD, this allows for a simple determination of infrared coefficients for gluon and ghost propagator as well as the running coupling. Results for general regulators are presented and interpreted in the context of an optimisation of the flow. The systematic extension of the truncation is discussed.

We begin by introducing the classical gauge fixed action for QCD in a general covariant gauge in four dimensions. Including also the ghost action, it is given by

$`S_{\mathrm{cl}}=\frac{1}{2}{\displaystyle \mathrm{tr}F^2}+\frac{1}{2\xi }{\displaystyle (D_\mu A^\mu )^2}+{\displaystyle \overline{C}_\mu D_\mu C}.`$ (1)
The flow equation approach relies on a momentum cut-off for the propagating degrees of freedom. Here, the momentum regularisation for both, the gluon and the ghost, is achieved by $`S_{\mathrm{cl}}S_{\mathrm{cl}}+\mathrm{\Delta }S_k`$ with

$`\mathrm{\Delta }S_k=\frac{1}{2}{\displaystyle A_a^\mu R_{\mu \nu }^{ab}A_b^\nu }+{\displaystyle \overline{C}_aR^{ab}C_b}.`$ (2)
The functions $`R`$ implement the infrared momentum cut-off at the momentum scale $`p^2k^2`$. The scale-dependence of the regularisation in (2) induces a flow equation for the effective action. With $`\varphi =(A,C,\overline{C})`$ and $`t=\mathrm{ln}k`$, the flow equation takes the explicit form

$`_t\mathrm{\Gamma }_k[\varphi ]={\displaystyle \frac{1}{2}}{\displaystyle d^4pG_{ab}^{\mu \nu }[\varphi ](p,p)_tR_{\mu \nu }^{ba}(p)}`$ (3)
$`{\displaystyle d^4pG^{ab}[\varphi ](p,p)_tR^{ba}(p)},`$
where $`G[\varphi ](p,q)=\left(\mathrm{\Gamma }_k^{(2)}+R\right)^1(p,q)`$ denotes the full regularised propagator. The flow is finite in both, the infrared and the ultraviolet, by construction. Effectively, the momentum integration in (3) only receives contributions for momenta in the vicinity of $`p^2k^2`$. Consequently it has a remarkable numerical stability. The flow solely depends on dressed vertices and propagators leading to consistent RG-scaling on either side of (3). The flow equations for propagators and vertices are obtained from (3) by functional derivatives.

To verify the Kugo-Ojima confinement criterion we have to determine the momentum behaviour of the propagators in the deep infrared region

$$k^2p^2\mathrm{\Lambda }_{\mathrm{QCD}}^2.$$
(4)
In this domain, quantum fluctuations are already integrated out and the physical quantities as well as general vertex functions will show no $`k`$ dependence. Moreover, (4) entails that we are already deep in the infrared regime. No new physics arises and the flow equation (3) only constitutes the shifting of the transition regime at momenta about $`k^2`$ towards smaller scales. Accordingly, Green functions $`\mathrm{\Gamma }_k^{(n)}\delta ^{(n)}\mathrm{\Gamma }_k/\delta \varphi _1\mathrm{}\delta \varphi _n`$ can be parameterised as

$`\mathrm{\Gamma }_k^{(n)}(p_i^2\mathrm{\Lambda }_{\mathrm{QCD}}^2)=z_n\widehat{\mathrm{\Gamma }}^{(n)}(p_i,p_i^2/k^2),`$ (5)
where $`p_i`$ stands for the momenta $`p_1,\mathrm{},p_n`$. Here $`z_n`$ are possibly $`k`$-dependent prefactors accounting for an RG-scaling that leaves the full action invariant. This reduces the remaining $`k`$-dependence in $`\widehat{\mathrm{\Gamma }}^{(n)}`$ to dimensionless ratios $`x=p^2/k^2`$. Hence, the interpolation between the physical infrared regime (4) and the regularised behaviour for $`p^2k^2`$ is $`k`$-independent. This implies a fixed point behaviour. So far, fixed point investigations within the exact renormalisation group have been put forward within derivative expansions. Here we provide a new truncation scheme for accessing the full momentum dependence of Green functions in the infrared limit.

Taking advantage of their simplified structure in the physical regime (4), we introduce the following parametrisation for the gluon and ghost two point function

$`\mathrm{\Gamma }_{k,A}^{(2)}(p^2)`$ $`=`$ $`z_AZ_A(x)p^2\mathrm{\Pi }(p)1\text{l}+\mathrm{longitudinal},`$
$`\mathrm{\Gamma }_{k,C}^{(2)}(p^2)`$ $`=`$ $`z_CZ_C(x)p^21\text{l},`$ (6)
where

$`Z(x)=x^\kappa \left(1+\delta Z(x)\right).`$
Here $`\mathrm{\Pi }_{\mu \nu }(p)=\delta _{\mu \nu }p_\mu p_\nu /p^2`$, $`1\text{l}_{ab}=\delta _{ab}`$ and $`Z,\delta Z,\kappa `$ without subscripts stand for $`Z_{A/C},\delta Z_{A/C},\kappa _{A/C}`$. The longitudinal part of the gluon two point function in (6) is not given here as it does not contribute to the Landau gauge flows considered below. The momentum dependence of (6) in the asymptotic regimes $`x1`$ and $`x1`$ fixes the functions $`\delta Z`$ in both regimes: the prefactor $`x^\kappa `$ encodes the leading behaviour for the momentum regime (4) where $`1x\mathrm{\Lambda }_{\mathrm{QCD}}^2/k^2`$, hence $`\delta Z`$ tends to zero. The $`k`$-independence of (6) at $`k=0`$ implies $`_t\kappa =0`$ and $`_tz=2\kappa z`$. In the limit of small $`x`$, we recover the trivial momentum behaviour $`\mathrm{\Gamma }_k^{(2)}(p^2)=m_k^2+z_1p^2+O(p^3)`$ as a consequence of the cut-off. Accordingly the functions $`\delta Z`$ have the limits

$`\delta Z(x\mathrm{})`$ $``$ $`0,`$
$`\delta Z(x0)`$ $``$ $`1+\mu x^{(1+\kappa )}+O(x^\kappa ).`$ (7)
We emphasise that the limit $`x\mathrm{}`$ implies $`k0`$ as we demand $`x\mathrm{\Lambda }_{\mathrm{QCD}}^2/k^2`$. In the intermediate region $`x1`$ the propagators show a regulator-dependent interpolation between the physical infrared behaviour and the trivial cut-off regime. Also

$`_t\delta Z(x)=2x_x\delta Z`$ (8)
reflects that the dependence of $`\delta Z`$ on $`k`$ stems from the interpolation property while its shape is independent of it. The vertices $`\mathrm{\Gamma }_k^{(n)}`$ in the present truncation are schematically given by

$`\mathrm{\Gamma }_k^{(n)}=z_nS_{\mathrm{cl}}^{(n)}`$ (9)
for $`n3`$. In the above parametrisation, the non-renormalisation of the ghost-gluon vertex at $`k=0`$ leads to

$`\kappa _A=2\kappa _C,\alpha _s={\displaystyle \frac{g^2}{4\pi }}{\displaystyle \frac{1}{z_Az_C^2}}.`$ (10)
It follows from (10) and $`_tz=2\kappa z`$ that the strong coupling constant achieves a fixed point in the infrared limit, $`_t\alpha _s=0`$. In particular, the $`z`$-factors for the three- and four-gluon vertices are proportional to $`\alpha _s^{1/2}z_A^{3/2}`$ and $`\alpha _sz_A^2`$, respectively. Moreover it follows within the truncation (9) that no mass term for the ghost is present corresponding to $`\mu _C=0`$ in (7). In summary, the truncation (6),(9) with the properties (7),(10) satisfy the RG properties of Landau gauge QCD as well as the truncated Slavnov-Taylor identities valid in the presence of the regulator.

We proceed with determining the strong coupling fixed point $`\alpha _s`$ and the infrared exponents $`\kappa `$ from (3) in the truncations put down in (6) . We fix the tensor structure of the regulators as $`R_{\mu \nu }^{ab}=R_A\delta _{\mu \nu }\delta ^{ab}`$ and $`R^{ab}=R_C\delta ^{ab}`$, leaving the scalar part $`R=zp^2r(x)`$ at our disposal. The integrated flow leads to an integral equation for $`\delta Z`$,

$`\delta Z_{A/C}(x)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{\pi ^2}}N{\displaystyle _x^{\mathrm{}}}{\displaystyle \frac{dx^{}}{x_{}^{}{}_{}{}^{2+\kappa }}}f_{A/C}(x^{}),`$ (11)
where

$`f_A(x)`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle _0^{\mathrm{}}}𝑑yy{\displaystyle _1^1}𝑑t{\displaystyle \frac{y}{u}}{\displaystyle \frac{(1t^2)^{3/2}𝒯_C(y)}{r_C(u)+Z_C(u)}},`$
$`f_C(x)`$ $`=`$ $`{\displaystyle \frac{x}{2}}{\displaystyle _0^{\mathrm{}}}𝑑yy{\displaystyle _1^1}𝑑t(1t^2)^{3/2}`$
$`\times {\displaystyle \frac{1}{u}}\left({\displaystyle \frac{y}{u}}{\displaystyle \frac{𝒯_C(y)}{r_A(u)+Z_A(u)}}+{\displaystyle \frac{𝒯_A(y)}{r_C(u)+Z_C(u)}}\right).`$
and $`u=x+y+2t\sqrt{xy}`$. We have also introduced the abbreviation $`2𝒯=k^2G(_tR)G`$, which reads explicitly

$`𝒯(y)={\displaystyle \frac{1}{y}}{\displaystyle \frac{\kappa r(y)yr^{}(y)}{(r(y)+Z(y))^2}}.`$
In (11), the integration over $`k`$ was traded for an integration over $`x`$, using $`z=\widehat{z}x^\kappa (\mathrm{\Lambda }_{\mathrm{QCD}}^2/p^2)^\kappa `$ with $`_t\widehat{z}=0`$. The purely gluonic diagrams in (3) are suppressed by positive powers of $`p^2/\mathrm{\Lambda }_{\mathrm{QCD}}^2`$ in the fixed point regime and do not contribute to (11).

For the sharp cut-off the kernel $`𝒯(x)`$ simplifies and becomes $`𝒯(x)=\delta (x1)/(zZ(1))`$. Then the $`t`$ and $`y`$ integrals in (11) can be performed analytically leading to a lengthy sum of hyper-geometric functions. For generic smooth cut-off, $`\kappa `$ and $`\alpha _s`$ are computed from (11) in the limit $`x0`$ as follows. On the right hand side of (11), the integrands $`f_{A/C}`$ are suppressed for small momenta $`x`$ and $`y`$ due to $`𝒯_{A/C}`$. This implies that $`\delta Z`$ can only contribute for sufficiently large arguments where it tends to zero due to (7). Therefore we can safely neglect $`\delta Z`$ under the momentum trace on the right hand side of (11) in a first step. Then, we proceed to fully consistent solutions $`\delta Z`$ of (11) by iteration. On the left hand side of the gluon equation the term $`\delta Z_A(x0)`$ is determined by (7). In the ghost equation, we have to carefully identify the terms. For $`\kappa _C>0`$ we read-off from (7) that $`\delta Z_C(x0)`$ diverges. However, $`\mathrm{\Gamma }_{k,C}^{(2)}(p^2)=z_Cx^{\kappa _C}(1+\delta Z_C(x))`$ has a Taylor expansion about $`x=0`$ due to the IR regularisation. Hence the coefficient of the constant term in $`\delta Z_C`$ is $`1`$.

The integral equation (11) can be solved numerically for $`\kappa `$ and $`\alpha _s`$. We have restricted ourselves to the domain $`\kappa _C(\frac{1}{2},1)`$ footnote . For general regulator $`r(x)`$ we obtain

$`\kappa _C=0.59535\mathrm{},\alpha _s=2.9717\mathrm{},`$ (12)
which agrees with the analytic DS result Lerche:2002ep . This might seem surprising given the completely different representations. However, the equivalence can be proven as follows: in the derivation of (12) we employed the truncation (6) and (9) with the additional constraint $`\delta Z0`$ under the momentum trace in (3). Integrating the flow (3) in this setting from $`0`$ to $`k`$ leads to

$`z_Ap^2x^{\kappa _A}\delta Z_A`$ $`=`$ $`\mathrm{Tr}\left[GS^{(3)}GS^{(3)}\frac{1}{2}GS^{(4)}G\right]_0^k,`$
$`z_Cp^2x^{\kappa _C}\delta Z_C`$ $`=`$ $`\mathrm{Tr}\left[GS^{(3)}GS^{(3)}\right]_0^k,`$ (13)
for the gluon and ghost two point function. Eq. (13) already implies an UV renormalisation at $`p^2=k^2`$. This is reflected by (11) being finite. Within the present truncation the renormalised DS equations used in Lerche:2002ep  are equivalent to (13). Beyond this truncation the DS equation and the flow equation differ. Apart from additional $`t`$-derivative terms, all vertices in (13) get dressed as distinguished from the DS equation, where both bare and dressed vertices are present. Consequently, the flow equation (3) carries the correct RG properties if the truncation for $`\mathrm{\Gamma }_k`$ does. In particular, equation (13) with the truncation (6) and (9) allows for a fully self-consistent renormalisation. Both the left and right hand sides of (13) transform the same way under RG scalings, including the purely gluonic terms dropped in (11). Indeed, it is the latter fact which makes their neglection fully consistent. The RG properties are most obvious for regulators proportional to the wave function renormalisations $`z`$, $`Rz\stackrel{~}{R}`$. Then, the propagators can be written as $`G=z^1G_0`$ where $`G_0`$ does not involve further wave function renormalisation factors. Also, the vertex functions involve the appropriate multiplicative powers of $`z`$ following from (10) and (9).

We proceed by iterating the above solution, also taking into account $`\delta Z0`$ on the right hand side of the flow. Then (11) can be used for a numerical iterative determination of $`\delta Z`$ and $`\kappa ,\alpha _s`$. We have computed $`\alpha _s`$ and $`\kappa `$ for various classes of regulators. The exponent $`\kappa _C`$ takes values in an interval bounded by the values (12) from above and the iterated sharp cut-off result $`\kappa _C0.539`$ from below. However, the values (12) are singled out by their physical relevance: $`\kappa `$ and $`\alpha _s`$ depend on the regulator. This reflects the fact that the truncation (6),(9) and the choice of $`R`$ define the approximation to the full problem. We expect that global extrema of $`\kappa (R),\alpha _s(R)`$ represent an optimised approximation based on (6),(9), see Litim:2001up . In general it is very difficult to translate this idea into properties of optimal regulators or even in statements about the values of optimal observables. Here we present a new idea how to extract this information: let us assume that we have identified a regulator $`R_0`$ which is an extremum in the space of regulators in the sense that

$`\left[\mathrm{Tr}\left(\delta R{\displaystyle \frac{\delta }{\delta R}}\right)\mathrm{\Gamma }_k^{(n)}\right]_{R=R_0}=0`$ (14)
for all $`\delta R`$. Eq. (14) implies that $`\kappa (R_0),\alpha _s(R_0)`$ are extrema, reflecting an optimisation of physical quantities. Note that the $`z_n(R_0)`$ by themselves are not physical quantities but they can be extrema up to RG scalings. However, (14) is more restrictive. For $`n=2`$, (14) is the additional demand, that $`\delta Z(x)`$ is an extremal curve. Hence the question arises whether solutions $`R_0`$ to (14) exist at all: the fixed point assumption underlying the whole analysis is valid for regulators $`R`$ that only lead to a suppression of infrared modes. This implies a monotonous interpolation between the physical infrared regime $`k^2p^2\mathrm{\Lambda }_{\mathrm{QCD}}^2`$ and the infrared regularised regime $`p^2k^2`$. Then the extremal curves exist and are given by vanishing $`\delta Z`$ footnote3 .

The above considerations allow us to identify the optimal values for $`\kappa _C(R)`$ and $`\alpha _s(R)`$ without further computation. The correction terms proportional to $`\delta Z`$ under the momentum trace are negligible near optimal regulators and we are left with the non-iterated equation (11) with an optimal regulator $`R_0`$. However, we have already shown that the non-iterated solutions are independent of the regulator and we conclude that $`\kappa _C(R_0),\alpha _s(R_0)`$ are given by (12). We are also in a position to device optimal regulators. For these regulators the effects of $`\delta Z`$ on the right hand side of (11) have to disappear. They need not to satisfy (14) which only applies to a small sub-set of optimal regulators. In Fig. 1 we show $`\delta Z_{A/C}`$ for the class of regulators

$`r_C(x)=z_C\theta ^1(1x),r_A(x)=\gamma {\displaystyle \frac{z_A}{x(1+x)}}.`$ (15)
The $`\gamma `$-dependence of $`\delta Z_A`$ is negligible on the scale of Fig. 1. The ghost regulator in (15) is a sharp cut-off while the gluon regulator is a rather soft mass-like cut-off. In the limit $`\gamma \mathrm{}`$ the solution of (11) approaches the optimal values (12). One verifies that all contributions from $`\delta Z`$ to the right hand side of (11) tend to zero for $`\gamma \mathrm{}`$ and we are left with the non-iterated result (12). Hence (15) with $`\gamma \mathrm{}`$ belongs to the class of optimal regulators. It is not in the extremal subset as only the ghost propagator satisfies (14) with $`\delta Z_C(x0)=0`$.

In summary, we have shown how the non-trivial momentum dependence of Green functions in the infrared sector is extracted within the exact renormalisation group. To that end we have developed new fixed point and optimisation techniques applicable to this situation. This reasoning has been applied to the infrared behaviour of ghost and gluon propagators in Landau gauge QCD. Our results provide further evidence for the Kugo-Ojima scenario of confinement: the ghost propagator is infrared enhanced while the gluon propagator develops a mass gap. In a truncation with dressed propagators, we have recovered previous results from Dyson-Schwinger equations in a very simple manner. This is noteworthy also in view of the qualitative differences between these approaches. Our method allows for a straightforward improvement of the truncation by adding vertex corrections and further terms with the correct RG properties. The same applies to dynamical quarks. In both extensions the correct RG scaling as well as the inherent finiteness of the integrated flow are most important. Further results and details will be presented elsewhere.

Acknowledgements

S.N. is supported by the DFG, contract SM70/1-1.


text:
# Segregation phases in a vibrated binary granular layer

## 1 Introduction

In segregation, excitation via flow or shaking can, counter-intuitively, cause an initially homogeneous mixture of grains to de-mix . Intriguingly, it does not always happen and the conditions for its occurrence are difficult to predict. Hence, a better understanding would have a major economic impact in the pharmaceutical, chemical processing and civil engineering industries . However, despite decades of research, a predictive model of the process has yet to emerge. More recently, segregation has received considerable attention from the physics community as an example of a challenging far from equilibrium system with pattern formation . Although many small scale laboratory studies involving vibration , avalanching in partially filled horizontal rotating drums  and stratification in vertically poured mixtures  have been investigated and a variety of mechanisms proposed , an understanding of the fundamental principles involved remains incomplete.

In our experiments  we have focused on quasi-2-dimensional layers of binary mixtures of particles which are horizontally driven by the stick–slip frictional interaction with the surface of a horizontal tray. Vibration of an initially homogeneous mixed binary layer can give rise to robust segregation of the constituents, depending on the total filling fraction of the layer. We call this parameter the *compacity*. As the compacity was increased the evolution of distinct phases was observed. We identified them as a *binary gas*, a *segregation liquid* and a *segregation crystal* using both microscopic and macroscopic measures to identify their properties. This horizontal and quasi-2-dimensional set up gives the practical advantage that any collective behavior is readily visualized, vertical compaction effects are effectively eliminated and the material is in contact with the forcing throughout the drive cycle.

For this class of quasi-2D binary granular systems a qualitative segregation mechanism has been suggested  using the idea of *excluded volume depletion* as in colloidal systems and binary alloys . This is in the spirit of the physics of complex systems where attempts have been made to unify descriptions of granular materials, colloids, gels and foams . We discuss the relevance of these ideas to our observations in some concluding remarks.

## 2 The experiment

A schematic diagram of the top view of the apparatus is presented in Fig. 1(a). It consisted of an accurately levelled and horizontal smooth rectangular tray, of dimensions $`(x,y)=180\times 90mm`$ with a flatness of less than $`\pm 5\mu m`$, on which particles were vibrated longitudinally. The tray was made from aluminum tool plate for increased rigidity and was mounted on a horizontal platform which was connected to an electro-mechanical shaker. Its motion was a close approximation to one-dimensional as it was constrained by four lateral high precision linear bush bearings. The dynamic displacement and acceleration of the shaking bed was monitored by a Linear Displacement Variable Transformer (LVDT) and a piezoelectric accelerometer. The tray’s acceleration was checked to be sinusoidal, to 0.1% in power, using the power spectrum of the accelerometer trace.

We have investigated segregation using a variety of mixtures but here focus on results for two specific types: *mixture A* (polystyrene spheres + sugar particles) and *mixture B* (phosphor-bronze spheres + poppy seeds). In particular, we have used mixture B for the detailed quantitative results. The details for the particles’ material properties is presented in Table 1.

At this point we define the *layer compacity* to be the total filling fraction of the system,

$$C=\frac{N_1A_1+N_2A_2}{xy},$$
(1)
where $`N_1`$ and $`N_2`$ are the numbers of species $`a`$ and $`b`$ in the layer, $`A_1`$ and $`A_2`$ are the two dimensional projected areas of the respective individual particles ($`A_1=(0.87\pm 0.15)mm^2`$ for poppy seeds and $`A_2=(1.77\pm 0.06)mm^2`$ for phosphor-bronze spheres) and $`x`$ and $`y`$ are the longitudinal and transverse dimensions of the rectangular tray. We have chosen to fix the forcing at a single frequency $`f=12Hz`$ and amplitude $`A=\pm 1.74mm`$ and explore the effect of the layer compacity.

Throughout this study, experiments were performed in a quasi-monolayer regime so that the smaller particles were rarely on top of the larger ones. This was particularly the case for the mixtures of poppy seeds and phosphor bronze spheres since the size ratio for these two types of particle is $`2:3`$. However, some overlapping of particles did occur so that compacities up to $`C1.1`$ were achieved in the high packing limit.

## 3 Single particle motion

As discussed above, the forcing was strictly sinusoidal but both types of particles were driven by a stick-and-slip interaction with the oscillatory surface of the container. These interactions induced stochastic motion and trajectories of single particles with the characteristics of a 2D random walk were obtained. Typical one minute long trajectories for a poppy seed and a phosphor-bronze sphere are shown in Fig. 1(b). The statistics of the motion were anisotropic with a preference for motion in the direction of the forcing. Note that the two types of particles responded to the driving differently, due to their mass, shape and surface properties. When many particles were present, collective motion dominated with diffusion biased in the direction of the forcing.

## 4 Granular segregation

In Fig. 1(c) and (d) we present two examples of segregation patterns for the two mixtures considered. The snapshots correspond to the segregated states which formed after the mixture was vibrated for a period of 15 minutes. The experimental runs were systematically started from a homogeneous mixed layer using the following method. Firstly, $`N_2`$ small particles were vibrated at large amplitudes, $`A\pm 5mm`$, creating an homogeneous and isotropic sub-monolayer. The large particles were then suspended above the layer, on a horizontal perforated plate with ($`57\times 28`$) $`2mm`$ diameter holes arranged in a triangular lattice. A shutter was then opened and the 1596 particles fell onto the layer of small particles, creating a near homogeneous binary mixture. Domains such as those shown in Fig. 1(c) and (d) then formed. The differences in the physical properties of the two mixtures illustrates that, in our system, patterned segregation is a robust phenomena.

In order to quantify the dynamics of growth of the segregated domains we focus on mixture B and establish two macroscopic measures: the average longitudinal width and the average number of domains of phosphor-bronze spheres. At this stage, the images of the whole layer ($`300\times 600`$ pixels) were processed such that regions of poppy seeds were given a value of 0 (black) and regions of spheres a value of 1 (white). For each acquired frame this was calculated by scanning through each of the 300 horizontal lines and the length, $`L_i`$, of each of the *white steps* was measured. We therefore define, $`\varphi =L`$, to be the average value of the domain width and the average number of domains, $`\eta `$, to be the number of such steps, per pixel line.

A typical time evolution of $`\varphi `$ for $`C=1.057`$, over a period of 20min, is presented in Fig. 2(a). It is characterised by a fast initial segregation growth, over timescales of the order of $`1min`$. During this period the mechanism of the segregation process is as follows. Initially, single large particles diffuse in a sea of the smaller ones. When two large particles come close together, the smaller particles cannot fit between them, and hence the pair is subjected to an asymmetric pressure that will keep it together. Subsequently, pairs may encounter others so that progressively larger clusters form. The unidirectionality of the driving induces an asymmetry in the segregated domains such that elongated domains of the larger particles (phosphor-bronze spheres or sugar particles) develop. The segregation domains were elongated in a direction perpendicular to the drive. The growth of the domains then exhibits an intermediate saturation level as can be seen in the plateau labelled $`P`$ in Fig. 2(a).

Further coarsening of the domains can take place, albeit at longer timescales. Two examples of the patterns at later times are presented in Fig. 2(b) and (c) at $`t_A=300s`$ and $`t_B=702s`$, respectively. During this period coarsening has occurred via the break-up of one of the stripes and consequent merging with its neighbours. In contrast to the fast initial segregation, this slow coarsening may be associated with large scale collective motion in the layer of poppy seeds.

Since there is a distinct separation of the timescales we will now focus on the *fast* initial segregation growth regime. In Fig. 3(a) and (b), we present typical time-series for $`\varphi `$ and $`\eta `$, respectively, for four different values of $`C`$, over the first $`4min`$ of vibration. The superimposed lines are fits of,

$$\varphi (C,t)=\mathrm{\Sigma }^\varphi b^\varphi .\mathrm{exp}\left(\frac{t}{t_s^\varphi }\right)\mathrm{and}\eta (C,t)=\mathrm{\Sigma }^\eta +b^\eta .\mathrm{exp}\left(\frac{t}{t_s^\eta }\right),$$
(2)
to the experimental time series of both $`\varphi `$ and $`\eta `$, where $`\mathrm{\Sigma }^\varphi `$ and $`\mathrm{\Sigma }^\eta `$ are the respective saturation levels and $`t^\varphi `$ and $`t_s^\eta `$, the *segregation times*, are timescales associated with the saturation of the domain growth. We denote $`\mathrm{\Sigma }^\varphi `$ by the *segregation level*. A result which may be inferred from both sets of time sequences in Fig. 3 is that, when the domains form, there is a progressive coalescence of clusters with an increase of the average width and corresponding decrease in the number of domains in the system.

## 5 Segregation Phase transition

We next investigated the dependence of the saturation levels, $`\mathrm{\Sigma }`$, and segregation times, $`t_s`$, on the compacity of the layer for both $`\varphi `$ and $`\eta `$. This was done by incrementally increasing the number of poppy seeds, $`N_1`$, in the layer, in measured steps, while the number of phosphor-bronze spheres, $`N_2=1596`$, held constant, i.e. $`C=C(N_1)`$. We chose to perform the experiments by changing the numbers of poppy seeds. This can be regarded analogous to the experimental protocol in colloid-polymer mixtures of increasing the concentration of the polymer in solution (see Section 8 and Ref. ). The qualitative behaviour presented here was found to be robust over a range of $`N_2`$.

The $`C`$-dependence of the segregation level is presented in Fig. 4(a). For low compacities, up to a value, $`𝒞_c0.647`$, the segregation level remains low and approximately constant around $`(3.25\pm 0.24)mm`$. This corresponds to approximately two sphere diameters and we do not classify these as domains; the system is deemed to be in a mixed or unsegregated state. This is essentially a collisional regime since there is enough free volume in the system so that particles collide and randomly diffuse, without aggregation. We denote this regime by a *binary gas*.

As $`C`$ is increased, past $`C_c`$, clusters of the larger particles form mobile segregation domains, of increasing $`\varphi `$, emerge in a nonlinear way. Above $`C_c`$, the clusters of larger particles were mobile and their movement was reminiscent of oil drops on water as they flowed, merged and split and the motion of the particles within the clusters was highly agitated. As $`C`$ was incrementally increased there was a gradual decrease of the overall motion of segregation domains so that they become less mobile and increasingly more stable. Also, for larger values of $`C`$, there was a significant decrease of the agitation of the spheres within the domains.

The average number of domains, shown in Fig. 4 remains approximately constant up to the same value $`C_c`$, after which it rapidly decreases as $`C`$ is increased. This is expected for increasingly wider domains since the number of phosphor-bronze spheres is constant. This scenario is consistent with the existence of a continuous phase transition for granular segregation of the binary mixture, as the layer compacity is incrementally increased, with a critical point for segregation, $`C_c0.647`$.

In Fig. 4(b) we present the segregation times obtained from the fits of Eq. (2) to $`\varphi (t)`$ and $`\eta (t)`$ as a function of $`C`$. Note that we only define $`t_s`$ for $`C>C_c`$, since this is the parameter range over which segregation occurs. As expected for a continuous phase transition, $`t_s`$ diverges, as $`C`$ is decreased from above, near $`C_c`$, a phenomena usually referred to as *critical slowing down*. This adds further evidence for the existence of a phase transition.

As the compacity is increased, past $`C_c`$, this slowing down of the segregation persists up to $`C=0.893`$. At this point there is a minimum in the segregation time, as shown in Fig. 3(c), with $`t_s=11.41s`$. For greater values of $`C`$ there is a qualitative change in the segregation behaviour and the layer takes increasingly longer to segregate. At this large values of $`C`$, stable crystalline stripes form which are separated by dense regions of poppy seeds.

The above scenario suggests that the segregated phase, for $`C>C_c`$, can be sub-divided into two further regimes: one in which the speed of segregation increases with $`C`$ ($`0.647<C<0.893`$) and another in which the opposite occurs ($`C>0.893`$). We denote the first by *segregation liquid* and the second by *segregation crystal*.

## 6 Binary gas, segregation liquid and segregation crystal phases

To further justify this idea of three segregation phases we turn to the structural information given by the radial pair distribution function for the large spheres, $`g(r)`$,

$$g(r)=A(r)\underset{i}{}\underset{ji}{}\delta (rr_{ij}),$$
(3)
where $`r_{ij}`$ is the separation of the $`i`$ and $`j`$th particles and the angled brackets denote a time average. $`A(r)`$ is a normalization constant such that $`g(r)=1`$ for a uniform distribution of particles. The data for the calculation of $`g(r)`$ was obtained 3min after the driving was switched on. Time averaging was then performed over 1min, for 1500 temporal realisations of the layer’s configuration.

In Fig. 5(a) we show three curves of $`g(r)`$ where each represents the segregation phases discussed above. Typical configurations of the spheres, for each of the phases, are presented in the zoomed photographs in Fig. 5(b),(c) and (d). Note that the photographs have been image processed such that only the phosphor-bronze spheres are visible. For the binary gas ($`C<0.647`$), $`g(r)`$ has a peak at $`r/d=1`$ which quickly decays at large distances, as expected for a disordered gas. Increasing $`C`$ results in a monotonic increase of the height of the first peak, $`g(d)`$. At intermediate compacities ($`0.647<C<0.893`$), *liquid-like* behavior is observed, with $`g(r)`$ peaked at 1, 2 and 3 particle diameters. The positions of neighbouring spheres are correlated and the maxima may be associated with shells of neighbours. The oscillations are rapidly damped, showing the decay of short-range order. This behaviour is commonly seen in hard sphere liquids and was first observed experimentally by Bernal . In the segregation crystal phase ($`C>0.893`$) two further peaks emerge near $`r/d=\sqrt{3}`$ and $`r/d=\sqrt{7}`$, characteristic of a hexagonally packed crystal, in two dimensions. At these high compacities, individual spheres in the segregation domains rattle inside the cage formed by its neighbours.

## 7 Oscillatory states

At the highest values of the compacity where robust crystalline stripes form, we have observed intriguing well defined oscillations of the segregation patterns. The phosphor-bronze stripes bend, periodically, backwards and forwards in a reproducible way. For each of the regions of poppy seeds in between two phosphor-bronze stripes, there is a complementary collective motion in the form of a single large vortex. The direction of rotation of the vortex is reversed every half-cycle of the oscillation. In Fig. 6(a)-(c) we show three photographs of the granular layer which were captured 40 min after having started from a homogenous mixture and correspond to the beginning of the cycle, half-cycle and full cycle points of one oscillation, respectively.

We illustrate the periodic nature of the oscillations, in Fig. 6(b) with a space time diagram of the process. This was constructed by sampling along a single line in the $`x`$-dimension, and progressively stacking the lines over a period of $`1.5hrs`$. The sampling line was positioned at 1/3 of the y-length of the tray as indicated by the solid white line superposed in Fig. 6(c). A well defined cycle is clearly observed and it has a period of $`T=286s`$. This suggests the existence of a strong coupling between the oscillation of the stripes through bending and the vortex inversion in the poppy seed’s layer. It is interesting to note that the timescale associated with these oscillations is $`3500`$ times slower than the driving (at $`12Hz`$) excluding the possibility of a simple resonance in the driving/granular-layer system. Hence, the mechanism which gives rise to this oscillatory phenomena remains, as yet, unexplained.

## 8 Discussion and conclusion

It is interesting to note that there are some analogies between segregation in our binary granular system and aggregation in binary colloids  and emulsions . In particular, when the size ratio between the large colloidal particles and the radius of gyration of the polymer is larger than $`0.3`$, colloid-polymer mixtures display colloidal gas, liquid and crystal phases of the large spheres . Aggregation in such systems is usually explained using an excluded volume depletion argument introduced by Asakura and Oosawa  in the context of binary hard spheres. This entropic argument is equivalent to the mechanistic view that if two large particles are close enough so that no other particle (or polymer) may fit between them, they will be subjected to an asymmetric osmotic pressure that leads to an effective attractive inter-particle force. In colloid-polymer mixtures the strength of the interaction can be tuned by changing the concentration of the polymer in solution, therefore inducing gas-to-liquid and liquid-to-crystal transitions. These resemble the granular phases we have observed by incremental increase of the compacity. However, an important difference is that our granular system is both driven and dissipative and therefore far from equilibrium. In addition, Brownian motion is an intrinsic part of the equilibrium dynamics of colloidal systems. Hence a direct connection between the two systems remains an open question.

In conclusion, our results provide quantitative evidence for the existence of three phases - binary gas, segregation liquid and segregation crystal - in granular segregation of horizontally excited binary mixtures. At high compacities, intriguing regular oscillatory states were observed. We show that segregation undergoes a continuous phase transition that occurs for compacities above a critical value and exhibits critical slowing down. We have presented both macroscopic and microscopic measures which are self-consistent. Analogies of this phase behaviour may be drawn with other depletion driven self-assembling binary systems. This raises the possibility of importing ideas from binary colloids in equilibrium to formulate new models for granular segregation.

The authors would like to thank A. Stephenson for valuable discussions. P.M.R. was supported by the Portuguese Foundation of Science and Technology. G.E. was funded by a studentship from the EPSRC. The research of T.M. is supported by an EPSRC Senior Fellowship. This work was completed while visiting the Isaac Newton Institute, Cambridge.