text:
# Numerically generated quasi-equilibrium orbits of black holes: Circular or eccentric?

## I Introduction and Summary

The late stage of inspiral of binary systems of neutron stars or black holes is of great current interest, both as a challenge for numerical relativity, and as a possible source of gravitational waves detectable by laser interferometric antennas. Because this stage, corresponding to the final few orbits and ultimate merger of the two objects into one, is highly dynamical and involves strong gravitational fields, it must be handled by numerical relativity, which attempts to solve the full Einstein equations on computers (see Refs. cook ; seidel  for reviews).

The early stage of inspiral can be handled accurately using post-Newtonian techniques, which involve an expansion of solutions of Einstein’s equations in powers of $`ϵ(v/c)^2Gm/rc^2`$, where $`v`$, $`m`$, and $`r`$ are typical velocity, mass and separation in the system, respectively. By expanding to very high powers of $`ϵ`$, increasingly accurate formulae can be derived to describe both the orbital motion and the gravitational waveform. Currently, results accurate through 3.5PN order ($`O(ϵ^{7/2})`$ beyond Newtonian gravity) are known jaraschafer98 ; jaraschafer99 ; djs00 ; bf00 ; bf01 ; abf01 ; dire2 .

An important issue in understanding the full inspiral of compact binaries is how to connect the PN regime to the numerical regime. This is a non-trivial issue because the PN approximation gets worse the smaller the separation between the bodies. On the other hand, because of limited computational resources, numerical simulations cannot always be started with separations sufficiently large to overlap the PN regime where it is believed to be reliable.

Numerical simulations of compact binary inspiral start with a solution of the initial value (I-V) equations of Einstein’s theory; these provide the initial data for the evolution equations (some simulations ggb02  solve in addition one of the six dynamical field equations). The initial state is assumed to consist of two compact objects (neutron stars or black holes) in an initially circular orbit. The circular-orbit condition is imposed by demanding that $`dr/dt=0`$ initially, where $`r`$ is a measure of the orbital separation. More precisely, the system is presumed to have an initial “helical Killing vector” (HKV), which corresponds to a kind of rigid rotation of the binary system. This amounts to ignoring initially the effects of gravitational radiation reaction. It also implies that the black holes are co-rotating, a condition which is astrophysically unlikely, albeit computationally advantageous. To make the problem tractable numerically, it is also generally assumed that the spatial metric is conformally flat. This approximation is usually justified by the neglect of radiation reaction in the initial state. For black hole binaries, suitable horizon boundary conditions must be imposed, while for neutron star binaries, hydrodynamical equations and an equation of state must be provided.

One important product of these initial value solutions is a relationship between the energy $`E`$ and angular momentum $`J`$ of the system as measured at infinity, and the orbital frequency $`\mathrm{\Omega }`$. As all quantities are well-defined and gauge invariant, they are useful variables for making comparisons with PN methods.

We have developed a formula for $`E(\mathrm{\Omega })`$ and $`J(\mathrm{\Omega })`$ using PN methods. Our analytic formula includes point-mass terms through 3PN order, but ignores radiation reaction. It also includes rotational energy and spin-orbit terms for the case in which the bodies are rotating. For black holes, tidal effects can be ignored. In contrast to previous work luc02 ; dgg02 , our formulae apply to general eccentric orbits, not just to circular orbits.

We then compare this formula with HKV numerical solutions for corotating binary black holes obtained by Grandclément et al. ggb02 , for the regime where the black holes are separated from the location of the innermost circular orbit by a factor of around two, where PN results might be expected to work well ($`Gm/rc^20.1`$). We find the following two results:

1. When we assume circular PN orbits, our PN formulae for $`E(\mathrm{\Omega })`$ and $`J(\mathrm{\Omega })`$ agree well with other PN methods, including those using resummation or Padé techniques. However all PN methods consistently and systematically underestimate the binding energy and overestimate the angular momentum, compared to the values derived from the numerical HKV simulations, by amounts that are larger than the spread among the PN methods.

2. When we relax the assumption of a circular orbit and demand only that $`dr/dt=0`$, our PN formula agrees extremely well with the numerical data. In this case the system is found to be initially at the apocenter of a slightly eccentric orbit. For values of $`Gm\mathrm{\Omega }/c^3`$ ranging from 0.03 to 0.06, corresponding to orbital $`v/c`$ between 0.3 and 0.4, or orbital separation between 10 and 6 $`Gm/c^2`$, perfect agreement with the binding energy can be obtained with eccentricities that range from 0.03 to 0.05.

If this apparent eccentricity is a real effect, there are two possible explanations. One is that the imposition of a helical Killing vector, which is equivalent to assuming a stationary binary star configuration as seen in a rotating frame, is not sufficient to guarantee a circular orbit in the relativistic regime. In Newtonian gravity, this assumption would be equivalent to imposing $`dr/dt=0`$, which only puts the orbit at a turning point. An additional condition $`d^2r/dt^2=0`$ is needed to enforce a circular orbit. But this is a dynamical statement, which is outside the realm of the I-V equations of Einstein’s theory. It may be that, in solving one of the dynamical field equations (the trace of the equation for $`K^{ij}/t`$) in addition to the four I-V equations, Grandclément et al. take care of this automatically.

The more likely possibility is that the approximations made in most numerical solutions, the main one being that the spatial metric is conformally flat, somehow introduce a systematic eccentricity into the orbit.

At present, these results can only be considered a hint of a possible eccentricity, however. It is entirely possible that circular PN orbits are consistent with the HKV results within the errors of the numerical simulations. This can only be decided if and when the numerical groups that carry out these simulations publish quantitative error bars determined by studying the convergence properties of the solutions as a function of grid size, domain size and other computational assumptions.

The remainder of this paper sketches the arguments that lead to these conclusions. Detailed derivations and formulae, and applications to neutron-star simulations, will be the subject of a future publication.

## II Orbits at the turning point in post-Newtonian gravity

In Newtonian gravity, the orbit of a pair of point masses may be described by the set of equations

$`p/r`$ $`=`$ $`1+e\mathrm{cos}(\varphi \omega ),`$
$`r^2\mathrm{\Omega }`$ $``$ $`r^2d\varphi /dt=(mp)^{1/2},`$
$`E`$ $`=`$ $`\mu (\dot{r}^2+r^2\mathrm{\Omega }^2)/2\mu m/r,`$
$`J`$ $`=`$ $`\mu |𝐱\times 𝐯|,`$ (1)
where $`p=a(1e^2)`$ is the semi-latus rectum ($`a`$ is the semi-major axis), $`\omega `$ is the angle of pericenter, $`m=m_1+m_2`$ is the total mass, $`\mu =m_1m_2/m`$ is the reduced mass, and $`E`$ is the total orbital energy (henceforth we use units in which $`G=c=1`$). A circular orbit corresponds to $`e=0`$, with $`r=a=`$ constant, $`\mathrm{\Omega }^2=m/a^3`$, $`E/\mu =a^2\mathrm{\Omega }^2/2m/a=(m\mathrm{\Omega })^{2/3}/2`$, and $`J/\mu =\sqrt{ma}=(m/\mathrm{\Omega })^{1/3}`$. However, if we demand only that the orbit be at apocenter, so that $`\dot{r}=0`$ only, we have $`\varphi =\omega +\pi `$, $`r=p/(1e)`$, $`\mathrm{\Omega }^2=(m/p^3)(1e)^4`$, and

$`E/\mu `$ $`=`$ $`{\displaystyle \frac{1}{2}}(1e^2)\left[{\displaystyle \frac{m\mathrm{\Omega }_a}{(1e)^2}}\right]^{2/3},`$
$`J/\mu m`$ $`=`$ $`\left[{\displaystyle \frac{m\mathrm{\Omega }_a}{(1e)^2}}\right]^{1/3}.`$ (2)
where $`\mathrm{\Omega }_a`$ is the angular velocity at apocenter. The energy and angular momentum expressions when the orbit is at pericenter may be obtained by letting $`\mathrm{\Omega }_a\mathrm{\Omega }_p`$ and $`ee`$.

At 1PN order, the equations of motion can be solved using either the direct approach of Wagoner and Will wagwill , or the “perturbed osculating orbit elements approach” of Lincoln and Will lincwill  to yield the orbit equations

$`\stackrel{~}{p}/r`$ $`=`$ $`1+\stackrel{~}{e}\mathrm{cos}(\varphi \omega )`$
$`+\stackrel{~}{\zeta }[(3\eta )+(1+9\eta /4)\stackrel{~}{e}^2`$
$`+(7/2\eta )\stackrel{~}{e}\mathrm{cos}(\varphi \omega )+3\stackrel{~}{e}\varphi \mathrm{sin}(\varphi \omega )`$
$`(\eta /4)\stackrel{~}{e}^2\mathrm{cos}(2\varphi 2\omega )]+O(\stackrel{~}{\zeta })^2,`$
$`r^2\mathrm{\Omega }`$ $`=`$ $`(m\stackrel{~}{p})^{1/2}[1\stackrel{~}{\zeta }(42\eta )\stackrel{~}{e}\mathrm{cos}(\varphi \omega )+O(\stackrel{~}{\zeta })^2],`$
$`E`$ $`=`$ $`\mu (\dot{r}^2+r^2\mathrm{\Omega }^2)/2\mu m/r`$
$`+{\displaystyle \frac{3}{8}}\mu v^4(13\eta )+{\displaystyle \frac{1}{2}}\mu (m/r)^2+{\displaystyle \frac{1}{2}}\mu \eta \dot{r}^2m/r`$
$`+{\displaystyle \frac{3}{2}}\mu (v^2m/r)(1+\eta /3)+O(m/r)^4,`$
$`J`$ $`=`$ $`\mu |𝐱\times 𝐯|\{1+{\displaystyle \frac{1}{2}}v^2(13\eta )+{\displaystyle \frac{m}{r}}(3+\eta )`$ (3)
$`+O(m/r)^2\},`$
where $`v^2=\dot{r}^2+r^2\mathrm{\Omega }^2`$, $`\eta =\mu /m`$, and $`\stackrel{~}{\zeta }=m/\stackrel{~}{p}`$. The limit $`\stackrel{~}{e}0`$ corresponds to a circular PN orbit. However, at higher PN orders, neither the orbital eccentricity $`\stackrel{~}{e}`$ nor the semilatus rectum $`\stackrel{~}{p}`$ is uniquely or invariantly defined. One definition of eccentricity used by Lincoln and Will lincwill  was that of a Newtonian orbit momentarily tangent to the true orbit (the “osculating” eccentricity); other authors ds88  define multiple “eccentricities” to encapsulate different aspects of non-circular orbits at PN order.

We define an alternative eccentricity and semilatus rectum according to:

$`e`$ $``$ $`{\displaystyle \frac{\sqrt{\mathrm{\Omega }_p}\sqrt{\mathrm{\Omega }_a}}{\sqrt{\mathrm{\Omega }_p}+\sqrt{\mathrm{\Omega }_a}}},`$
$`\zeta {\displaystyle \frac{m}{p}}`$ $``$ $`\left({\displaystyle \frac{\sqrt{m\mathrm{\Omega }_p}+\sqrt{m\mathrm{\Omega }_a}}{2}}\right)^{4/3},`$ (4)
where $`\mathrm{\Omega }_p`$ is the value of $`\mathrm{\Omega }`$ where it passes through a local maximum (pericenter), and $`\mathrm{\Omega }_a`$ is the value of $`\mathrm{\Omega }`$ where it passes through the next local minimum (apocenter). These definitions have the following virtues: (1) they reduce precisely to the normal eccentricity $`e`$ and semilatus rectum $`p`$ in the Newtonian limit, as can be seen from Eqs. (1); (2) they are constant in the absence of radiation reaction; (3) they are somewhat more directly connected to measurable quantities, since $`\mathrm{\Omega }`$ is the angular velocity as seen from infinity (eg. as measured in the gravitational-wave signal) and one calculates only maximum and minimum values, without concern for the coordinate location in the orbit; and (4) they are straightforward to calculate in a numerical simulation of orbits without resorting to complicated definitions of “distance” between bodies. They have the defect that, when radiation reaction is included, they are not local, continuously evolving variables, but rather are some kind of orbit-averaged quantities (for this reason, they may not be as “covariant” as they seem – this issue is under investigation). Nevertheless, when an eccentric orbit decays and circularizes under radiation reaction the definition of $`e`$ has the virtue that it tends naturally to zero when the orbital frequency turns from ocillatory behavior to monotonically increasing behavior (i.e. the maxima and minima merge).

By virtue of these definitions, $`\zeta `$ has the further property that

$$\zeta =\left(\frac{m\mathrm{\Omega }_p}{(1+e)^2}\right)^{2/3}=\left(\frac{m\mathrm{\Omega }_a}{(1e)^2}\right)^{2/3}.$$
(5)
It is then simple to show that the relation between $`e`$ and $`\zeta `$ and the corresponding quantities used in the Wagoner-Will solution (3) is $`e=\stackrel{~}{e}\{1+\stackrel{~}{\zeta }[9/2\eta +(13\eta )\stackrel{~}{e}^2]\}`$, and $`\zeta =\stackrel{~}{\zeta }\{14\stackrel{~}{\zeta }[1\eta /3+(1/3\eta )\stackrel{~}{e}^2]\}`$. Applying these definitions to Eqs. (3), we find the 1PN expression for $`E(\mathrm{\Omega })`$ and $`J(\mathrm{\Omega })`$ for a non-circular orbit, expressed in terms of $`\mathrm{\Omega }`$ at the apocenter or pericenter:

$`{\displaystyle \frac{E}{\mu }}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1e^2)\zeta +[{\displaystyle \frac{3}{8}}+{\displaystyle \frac{\eta }{24}}+({\displaystyle \frac{5}{12}}+{\displaystyle \frac{\eta }{12}})e^2`$
$`+({\displaystyle \frac{1}{24}}{\displaystyle \frac{\eta }{8}})e^4]\zeta ^2,`$
$`{\displaystyle \frac{J}{\mu m}}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{\zeta }}}\left\{1+\left[{\displaystyle \frac{3}{2}}+{\displaystyle \frac{1}{6}}\eta \left({\displaystyle \frac{1}{6}}{\displaystyle \frac{1}{2}}\eta \right)e^2\right]\zeta \right\}.`$ (6)
Note that in the circular orbit limit ($`e0`$), $`E`$ and $`J`$ satisfy the expected relation $`dE/d\mathrm{\Omega }=\mathrm{\Omega }dJ/d\mathrm{\Omega }`$. We have extended this result to 3PN order using the 3PN orbit equations of Blanchet and Faye bf01 , using both harmonic-gauge and ADM gauge expressions, but neglecting radiation-reaction terms at 2.5PN order (details will be given in a future publication).

For those numerical simulations in which the bodies are assumed to be corotating, we must include a number of rotational effects. First is the energy of rotation of the individual corotating bodies, of order $`\delta E_{\mathrm{rot}}/\mu I\mathrm{\Omega }^2/\mu (m/r)(R/r)^2`$, where $`R`$ is the characteristic size of the body. For compact bodies, $`RM`$, so these effects are equivalent to 2PN order and higher. For black holes, we will use the standard procedure of splitting the mass of each body into its irreducible mass and its rotational energy using the Kerr formulae $`M=M_{\mathrm{irr}}/[14(M_{\mathrm{irr}}\mathrm{\Omega })^2]^{1/2}`$ and $`S=4M_{\mathrm{irr}}^3\mathrm{\Omega }/[14(M_{\mathrm{irr}}\mathrm{\Omega })^2]^{1/2}`$ luc02 ; dgg02 . Also, since the sequences of numerical simulations hold $`m_{\mathrm{irr}}=(M_{\mathrm{irr}})_1+(M_{\mathrm{irr}})_2`$ fixed as they vary $`\mathrm{\Omega }`$, we will expand the masses that appear in the Newtonian orbital contribution to the energy and angular momentum about $`m_{\mathrm{irr}}`$. These contributions together yield

$`\delta E_{\mathrm{rot}}/\mu `$ $`=`$ $`(2/\eta )(13\eta )(m_{\mathrm{irr}}\mathrm{\Omega })^2`$
$`+(6/\eta )(15\eta +5\eta ^2)(m_{\mathrm{irr}}\mathrm{\Omega })^4`$
$`(1e^2)({\displaystyle \frac{2}{3}}\eta )\zeta (m_{\mathrm{irr}}\mathrm{\Omega })^2,`$
$`\delta J_{\mathrm{rot}}/\mu m_{\mathrm{irr}}`$ $`=`$ $`(4/\eta )(13\eta )(m_{\mathrm{irr}}\mathrm{\Omega })`$ (7)
$`+(8/\eta )(15\eta +5\eta ^2)(m_{\mathrm{irr}}\mathrm{\Omega })^3`$
$`+{\displaystyle \frac{1}{\sqrt{\zeta }}}(m_{\mathrm{irr}}\mathrm{\Omega })^2({\displaystyle \frac{4}{3}}2\eta ),`$
where henceforth, $`\zeta =[m_{\mathrm{irr}}\mathrm{\Omega }_a/(1e)^2]^{2/3}`$, and $`\mu `$ and $`\eta `$ are expressed in terms of irreducible masses. The first two terms in each of Eqs. (7) are the rotational terms, expanded to second order, while the third comes from expanding the Newtonian orbital term.

Similarly spin-orbit and spin-spin effects are of order $`\delta E_{\mathrm{S}.\mathrm{O}.}/\mu LS/\mu r^3`$ and $`\delta E_{\mathrm{S}.\mathrm{S}.}/\mu S_1S_2/\mu r^3`$ where $`L`$ denotes the orbital angular momentum and $`S`$ denotes the bodies’ spin. For co-rotating bodies, these effects can be shown to be of order $`\delta E_{\mathrm{S}.\mathrm{O}.}/\mu (m/r)(mR^2/r^3)`$ and $`\delta E_{\mathrm{S}.\mathrm{S}.}/\mu (m/r)(mR^4/r^5)`$, which for compact bodies are equivalent to 3PN and 5PN order, respectively. Henceforth, we will ignore the spin-spin terms. Assuming co-rotating bodies with spins aligned with the orbital angular momentum, including both the direct and orbital effects of the spin-orbit terms kww ; kidderspin  as well as their effect on our definitions of $`e`$ and $`\zeta `$, we find

$`{\displaystyle \frac{\delta E_{\mathrm{S}.\mathrm{O}.}}{\mu }}`$ $`=`$ $`\left({\displaystyle \frac{16}{3}}+12\eta +{\displaystyle \frac{4}{3}}(511\eta )e^2{\displaystyle \frac{4}{3}}(12\eta )e^4\right)`$
$`\times (m_{\mathrm{irr}}\mathrm{\Omega })\zeta ^{5/2}.`$
$`{\displaystyle \frac{\delta J_{\mathrm{S}.\mathrm{O}.}}{\mu m_{\mathrm{irr}}}}`$ $`=`$ $`\left({\displaystyle \frac{40}{3}}+30\eta {\displaystyle \frac{2}{3}}(411\eta )e^2\right)(m_{\mathrm{irr}}\mathrm{\Omega })\zeta .`$ (8)
Apart from the generalization to eccentric orbits, these methods parallel exactly those of Blanchet luc02 .

For tidal interactions, the contribution to the orbital energy is of order $`\delta E_{\mathrm{tidal}}/\mu (m/r)(R/r)^5`$. For black holes, tidal effects are thus equivalent to 5PN order, and will be neglected. For neutron stars, $`R5M`$, and tidal effects must be included, however it is sufficient to use Newtonian gravity to calculate them. These will be discussed elsewhere.

## III Comparison with numerical simulations of corotating binary black holes

We combine the PN, rotational and spin-orbit contributions to $`E`$ and $`J`$, set $`\mathrm{\Omega }_a=\mathrm{\Omega }`$ so that the stars corotate at the orbital angular frequency at apocenter, set $`\eta =1/4`$ for equal-mass black holes, and plot $`(Em_{\mathrm{irr}})/m_{\mathrm{irr}}`$ and $`J/m_{\mathrm{irr}}^2`$ versus $`m_{\mathrm{irr}}\mathrm{\Omega }`$, where $`m_{\mathrm{irr}}`$ denotes the total irreducible mass. The results are shown in Figures 1 and 2. For $`e=0`$, we show results both from the 3PN orbital expressions, and from truncated 2PN orbital expressions, while for $`e=0.03`$ we plot the full 3PN results. The crosses denote the corresponding numerical data of Grandclément et al.ggb02 ; the other curve is the 3PN resummation result of Damour et al.dgg02  for circular orbits. We note that all circular-orbit PN estimates, including conventional 2PN and 3PN approaches luc02 , Padé and resummation approaches dgg02 , and our approach, agree well among themselves, but are systematically displaced from the numerical points, both for $`E`$ and for $`J`$. This can also be seen in Figures 5 and 6 of dgg02 . By contrast, with $`e0.03`$, our PN estimate fits the curves perfectly for small values of $`m_{\mathrm{irr}}\mathrm{\Omega }`$. This fit is merely illustrative: if the numerical simulations do have some eccentricity, there is no reason why the eccentricity should be the same for each numerical point. Figure 3 shows the values of $`e`$ that best fit the numerical data over the range of $`\mathrm{\Omega }`$ considered.

The conclusions are not substantially different if we use the Wagoner-Will type eccentricity $`\stackrel{~}{e}`$. Whatever definition one uses, $`e=0`$ gives a worse fit than a finite eccentricity.

## IV Discussion

High-order post-Newtonian approximations for circular orbits appear to give results for the energy and angular momentum of corotating binary black holes well away from the innermost orbit that are in excellent agreement – with each other. This is illustrated in Figure 4, which plots the energy from four PN results quoted by Damour et al dgg02 , and from our 2PN and 3PN results, in the range around $`\mathrm{\Omega }m_{\mathrm{irr}}0.32`$. Apart from the one 2PN result quoted by Damour et al., all are within 0.5 percent of each other. This could be viewed as an estimate of the accuracy of the PN approximation in this regime. But all results are four percent displaced from the numerical HKV data. In the absence of a quantitative estimate of the accuracy of the numerical simulations in this regime, it is difficult to decide if this difference is a signal of a physical effect, such as the small eccentricity suggested by our work.

One way to test whether the orbits represented by the HKV numerical simulations are really eccentric would be to evolve the orbits numerically for a short period of time, say 1/4 of an orbit ($`\delta \varphi =\pi /2`$). For a nearly circular orbit of equal mass bodies, gravitational radiation reaction should cause the parameter $`\zeta =m/p=(m\mathrm{\Omega })^{2/3}`$ to evolve according to $`d\zeta /d\varphi =(16/5)\zeta ^{7/2}`$ (see, for example Eq. (3.11a) of lincwill ). Thus, over a quarter of an orbit, $`\zeta `$ should change by approximately $`\delta \zeta /\zeta (8\pi /5)\zeta ^{5/2}(8\pi /5)(m\mathrm{\Omega })^{5/3}0.016`$, for $`m\mathrm{\Omega }0.032`$. Hence the orbital separation should decrease by about 1.6 percent or the angular velocity should increase by 2.4 percent. By contrast, if the orbit has an eccentricity of 0.03 and is at apocenter, then in a quarter of an orbit, the separation should decrease by 3 percent, while the angular velocity should increase by 6 percent. Even with radiation damping, the orbit should pass through a distinct pericenter when $`\delta \varphi \pi `$, and the orbital separation should increase, while the angular velocity decreases.

###### Acknowledgements.

We are grateful to Luc Blanchet, Mark Miller and Wai-Mo Suen for useful discussions, and to Eric Gourgoulhon and Thibault Damour for providing the numerical and PN data from Figs. 5 and 6 of dgg02  for our comparisons. This work was supported in part by the National Science Foundation under grant No. PHY 00-96522. T.M. was supported in part by an internship from the Centre Nationale de la Recherche Scientifique.


format: text

text:
# A test of arm–induced star formation in spiral galaxies from near–IR and H𝛼 imaging

## 1 INTRODUCTION

Two of the leading theories of star formation in spiral galaxies use the concept of a density wave, either as the actual triggering mechanism, or as a means of organisation of star–forming material. The first of these theories has been termed the large scale galactic shock scenario (Roberts 1969; see also Shu et al. 1972; Tosa 1973; Woodward 1975; Nelson & Matsuda 1977). This model hypothesises that the gas settles into a quasi–stationary state, with a velocity and density distribution that is driven by the gravitational field of the galaxy. The gas response can be non–linear to an imposed azimuthal sinusoidal potential, if relative motion between the density wave and the cold interstellar medium (ISM) is supersonic (Binney & Tremaine 1987). This leads to the formation of a shock near the trailing edge of spiral arms, assuming that the region is inside the corotation radius, which compresses the gas to densities at which stars can form. In observations of spiral galaxies the shock is thought to be characterised by dust–lanes seen on the trailing edges of arms. The time delay needed for the onset of star formation after the compression of the gas implies that star–forming regions should be seen towards the leading edges of arms.

Many of the recent advances in this area have arisen from studies of the atomic and molecular gas in nearby spiral galaxies, through HI and CO line emission (e.g Nakai et al. 1994; Rand 1993, 1995). These studies reveal streaming velocities of gas through the spiral arms and find offsets between the peaks of the gas density and the old stellar population, in agreement with the predictions of the large scale shock scenario.

The alternative picture to this form of triggered star formation is stochastic star formation. Öpik (1953) first hypothesised that a supernova explosion could trigger star formation. A model proposed and developed by several authors (e.g. Gerola & Seiden 1978; Seiden & Gerola 1982; Seiden 1983; Jungwiert & Palous 1994; Sleath & Alexander 1995, 1996) suggests that the dominant process for forming stars is stochastic self–propagating star formation, and not density wave triggering. In this model, density waves are only responsible for the organisation of the ISM and stars, and for concentrating new HII regions along spiral arms (Elmegreen & Elmegreen 1986; Elmegreen 1993). Thus, star formation rate efficiency (i.e. normalised to unit mass of disc material) should be unaffected by location in arm or interarm regions, in this model. This is in clear distinction to the predictions of the large–scale shock model.

Previous H$`\alpha `$ studies of spiral galaxies (Kennicutt 1989, 1998a; Kennicutt, Tamblyn & Congdon 1994; and the review by Kennicutt 1998b) have mainly looked at global star–formation, particularly the form of the Schmidt law (Schmidt 1959, 1963). In this paper we are looking at localised star formation in disc galaxies as well as global properties of star formation, and comparing the distribution of star formation with the underlying old stellar population. The main aim of this paper is the analysis of star formation efficiencies in arm and interarm regions in spiral galaxies.

Other studies of star formation efficiencies include Lord & Young (1990) and Tacconi & Young (1986). Lord & Young (1990) looked at the molecular, neutral, and ionized hydrogen distributions in M51. They compared a ratio of massive star formation rates (MSFR) to gas surface density between arm regions and interarm regions and found a higher ratio in the arm regions. Tacconi & Young (1986) performed a similar analysis for NGC 6946 and also found that star formation is more efficient in arm regions than in interarm regions. Cepa & Beckman (1990) and Knapen et al. (1992) again looked at star formation efficiencies in arm and interarm regions, but with high spatial resolution. Cepa & Beckman (1990) compare the H$`\alpha `$/$`\sigma _{HI}`$ ratio (where $`\sigma _{HI}`$ is the HI surface brightness) in the arm and interarm regions of NGC 3992 and NGC 628. Knapen et al. (1992) compare the H$`\alpha `$/$`\sigma (H_2+HI)`$ ratio in arm and interarm regions. Both Knapen et al. (1992) and Cepa & Beckman (1990) explicitly find evidence for spiral arm triggering of star formation, via a non–linear dependence of SFR on gas density. Wyder, Dolphin & Hodge (1998) obtained HST WFPC2 V band (F555W) and I band (F702W) data of the spiral galaxy NGC 4321. From the I–V colour of the arm and interarm regions, they concluded that the SFR over the last 5 Myr has been approximately 4 times larger in the arm regions than in the surrounding interarm regions. Also, Knapen (1998) performed a study of HII regions in M100, using a new H$`\alpha `$ image. He found that the arms collect HII regions in a similar way to that described by Elmegreen & Elmegreen (1986). He concluded that spiral arms may trigger some star formation, but do not affect the HII region luminosity function or mass distribution. Finally, Alonso–Herrero & Knapen (2001) used the specific SFR or SFR per stellar mass (i.e. a comparison of Pa$`\alpha `$ flux with H band continuum emission), and found that while the specific SFR in the central region does not vary statistically with galaxy type, it is higher in barred than in non-barred galaxies.

The remainder of this paper is arranged as follows. Section 2 describes the observations and data analysis; section 3 describes local star–formation rates with respect to properties of spiral structure; and section 4 contains our conclusions.

## 2 OBSERVATIONS AND DATA REDUCTION

In this sample, 14 of the galaxies were selected as part of a program designed to understand the underlying structure of galaxies. The selection of these galaxies, and near–IR observations and data reduction, is described in Seigar & James (1998a). Briefly, the full galaxy sample is comprised of 45 spirals of types Sa–Sdm, with inclinations $``$ 45 and diameters $``$1$`\stackrel{}{.}`$5, the latter enabling them to be imaged in the field of the UKIRT near–IR camera IRCAM3. A further 5 galaxies have been observed with the UKIRT Fast Track Imager (UFTI), with the same selection criteria. The final galaxy was observed with INGRID on the William Herschel Telescope (WHT), again with the same selection criteria, except its diameter was chosen to be $``$4.

The optical data presented here are broad–band (R) and narrow–band (redshifted H$`\alpha `$) images of the 20 galaxies described above. The observations were made using the 1.0m Jacobus Kapteyn Telescope (JKT) on La Palma. Some of the galaxies were observed on the nights 1998 May 14–20, with some additional observations being taken during a later run, 1999 February 4–11, using the 1024$`\times `$1024 Tek4 CCD camera which has 0$`\stackrel{}{.}`$331 pixels and a field of 5$`\stackrel{}{.}`$65$`\times `$5$`\stackrel{}{.}`$65. The rest of the objects were observed between February 2000 and January 2002, using the 2048$`\times `$2048 SITe2 CCD camera, which has 0$`\stackrel{}{.}`$33 pixels and a field of 11$`\stackrel{}{.}`$3$`\times `$11$`\stackrel{}{.}`$3. The filters used were standard Harris R filter, and narrow–band filters with the following central wavelengths in Angstroms (full–width at half maximum in brackets): 6594(44), 6626(44), 6656(44), 6686(44), 6695(53), 6712(42), 6727(48), 6785(58). Measured transmission curves are available for all these filters, and were used in selecting the optimum filter for each galaxy redshift, and in flux calibrating the images, as described below.

All 20 galaxies were observed in redshifted H$`\alpha `$ filters, with between one and four 1200 second integrations for each galaxy. All were also observed at R with integration times of 300–900 seconds. In addition, the following standard stars, selected from Landolt (1992), were observed: PG1047+003 A, B & C; PG1323–086 A, B & C; PG1633+099 B, C & D; SA111–773; and SA111–775. Two spectrophotometric standards from Oke (1990), Feige 34 and G191B2B, were also observed.

Data reduction made use of standard STARLINK image reduction packages. All images were bias subtracted, and then flat fielded using twilight sky flats. There was no evidence for any fringing effects in broad– or narrow–band images. The H$`\alpha `$ images were continuum–subtracted using scaled R band images. The scaling factor applied to the R band images was calculated photometrically, using the ratio of fluxes of red standard stars from the Landolt (1992) list. These values were checked using the spectrophotometric standards from Oke (1990), to ensure that the red standards did not have strong spectral features in the passband of the H$`\alpha `$ filter, and by numerically integrating under the filter profiles. All methods gave very similar results, but the red standards were adopted as best representing the colour of the galaxy continuum light. The throughputs of the narrow band filters to continuum light were measured in this way to be between 36 and 64 times lower than that of the R band filter. The narrow–band images were continuum subtracted by normalising broad– and narrow–band images to the same integration time, dividing the former by the ratio of the throughputs, and subtracting this scaled broad–band image.

Calibration of the H$`\alpha `$ photometry required calculation of the effective throughput of the narrow–band filter to H$`\alpha `$ light. This was done by calculating the throughputs of both narrow– and broad–band filters at the wavelength corresponding to redshifted H$`\alpha `$, for each galaxy, from the measured transmission curves. The narrow–band transmission was then corrected by subtracting the broad–band transmission, scaled by the ratio of continuum throughputs (36–64, as given above). Given the large size of this factor, the continuum subtraction works well, even though the R band filter contains the H$`\alpha `$ line, and the only effect of this is to reduce the effective throughput for H$`\alpha `$ emission by a few percent. Effective throughputs calculated in this way were in the range 28–50%, typically $``$45%.

This continuum subtraction procedure worked extremely well, with the only problems being flat–fielding residuals which were apparent in some of the continuum (R band) images taken in bright moonlight, and a slight change in plate–scale between the broad– and narrow–band images, which was not significant over the size of the galaxies in this study. The accuracy of the continuum subtraction was confirmed by the almost complete removal of old–stellar light (non–star–forming bulges etc) in the resulting line images. Foreground stars were also removed fairly accurately, but residuals of a few percent were sometimes left due to the differing colours compared with integrated galaxy light. There was also a small residual sky level, which was subtracted before photometric measurements were made.

The photometry was obtained using the STARLINK software package GAIA. This extracted H$`\alpha `$ fluxes as raw counts, which were then corrected for airmass, using a correction based on the R band observations of the standard stars listed above. A correction was then applied for the effective throughput of the narrow–band filter to the redshifted H$`\alpha `$ radiation as described above.

Table 1 contains the following entries: Galaxy name (column 1), galaxy classification (column 2), heliocentric redshift in km s<sup>-1</sup> (column 3) and the telescope and instrument with which they were observed (column 4).

## 3 RESULTS

In this section we look at how density waves affect the local star formation rates within the discs of individual galaxies, even though other factors such as cold gas mass or galaxy environment may dominate the differences in star formation rate between galaxies. To test the local effect of arms, we compared the degree of concentration of H$`\alpha `$ and K band light in the arms of the present galaxy sample. We use the H$`\alpha `$ data as an indicator of star formation rate and the K band data as a tracer of where the density waves are (the K band largely reflects old stars, which is a good tracer of the overall stellar mass density, and hence of density waves within discs). If the simple model of Elmegreen (1993) is true, the density waves should simply gather up all types of disc material to the same extent, and the same arm–to–interarm contrast should be seen in all disc tracers. Then, arms would be just as prominent in the near–IR emission from old stars as in the H$`\alpha `$ emission from star–formation regions. If, however, spiral density waves trigger star formation to a significant extent, then the arm–to–interarm contrast could be substantially enhanced in the star–formation tracer relative to the old stellar population. Note that this rather naive argument neglects the effects of extinction, but since the optical light is likely to be more affected than the near–IR, any relative enhancement in H$`\alpha `$ arm contrasts compared to say the K band should be a strong indication of arm–induced star formation.

### 3.1 A comparison of H$`\alpha `$/K ratios for arm and interarm regions

In order to measure any enhancement in H$`\alpha `$ arm contrasts compared to K band arm contrasts, we came up with the following conservative test, which is based on the technique used to measure star–formation efficiency by Knapen et al. (1996), to identify arm–induced star formation in the present sample of galaxies. For those galaxies which showed well–defined spiral arms in their K band images, we rebinned the K images to the same pixel scale as the H$`\alpha `$ images, and shifted them to overlay the H$`\alpha `$ images, using stellar centroids as references. This procedure should be accurate to a fraction of a pixel, or better than 0$`\stackrel{}{.}`$1. We then used the ‘ARD region’ option within the STARLINK GAIA package to define a polygon lying around each K band arm, without reference to any of the optical images. Elliptical regions were also fitted by eye to define the disc of the galaxy, excluding the central bulge, and extending to the outer ends of the spiral arms. The images were sky subtracted, and then the number of detected counts within each of the arm polygons, and the overall disc, were calculated. GAIA then permits the apertures to be imported into the H$`\alpha `$ images, which had been aligned exactly, such that the apertures lie over the same physical regions of the galaxies at both K and H$`\alpha `$. Photometry was then obtained for all regions (arms, total disc and nucleus) in the H$`\alpha `$ images in exactly the same way as for the K images. Finally, ratios were calculated of H$`\alpha `$ flux divided by K band flux for all of the regions. The test is then to see whether this ratio is larger for the arm regions than for the disc generally, where the latter ratio is taken from the total disc minus the central region. Any obvious foreground stars are removed from from the arm and disc regions before this comparison is done.

If arms represent regions where all disc material is concentrated by an equal factor, then the H$`\alpha `$/K ratio would be the same in arms as in the disc overall. If there were no connection between star formation and K band structure, then this ratio would be lower in the selected arm regions, since these were chosen to have higher than average K surface brightnesses. However, in the majority of cases, we found H$`\alpha `$/K ratios to be significantly higher in the K band arms than in the discs generally. We measured a total of 49 arm regions in 20 galaxies (see Table 2); 38 (76%) of these arms had H$`\alpha `$/K ratios greater than that of the disc of the same galaxy. Averaging the ratios over all arms observed in a given galaxy, we find (Table 2 final column) that 19 of the 20 galaxies show a net enhancement in H$`\alpha `$/K ratio within their arm regions. The errors in table 2 were calculated by both adding and subtracting the uncertainty in the sky background from the H$`\alpha `$ images and performing the same analysis for the original H$`\alpha `$ image, the H$`\alpha `$ image subtracted by the sky background uncertainty, and the H$`\alpha `$ image added by the sky background uncertainty. The error in the sky background was calculated using both the formal uncertainty from pixel–pixel noise, and larger–scale systematic variations in the sky structure. Combining the results for all 49 arms, the median ratio was 1.47$`\pm `$0.11 times higher in the arms than in the corresponding discs. We have taken the median ratio in order to give less weight to galaxies which may be affected by starbursts. Given the conservative nature of the test, this represents a highly significant finding of triggering of star formation within spiral arms. We checked that the location of the arms would not have been changed had we used J rather than K images, and thus our results are not affected by the possibility that hot dust may contribute to galaxy light at K (James & Seigar 1999).

Two of the galaxies warrant further comment. UGC 3053 exhibits 3 well–defined arm segments in its K band image, but there is little or no H$`\alpha `$ luminosity from HII regions associated with these arms. This is reflected in the low H$`\alpha `$/K ratios for all three arm regions in this galaxy, which thus appears to have no arm–induced star formation according to this test. UGC 6332 has two tightly wound arms, which could be classified as a single ring, and which appear highly symmetrical in the K image. However, this symmetry breaks down completely in H$`\alpha `$ light, since one arm is clearly associated with a string of HII regions, and the other with no detectable HII regions at all. As a result, these two arms give rise respectively to the second–highest, and to easily the lowest, of the 49 measured H$`\alpha `$/K ratios. This asymmetry could be an extinction effect, due to the offset between the dust lane and the peak of the induced star formation.

In this analysis we have ignored the result that the K band light may have a contribution of up to 30% from star formation (James & Seigar 1999). An effect such as this will tend to increase the amplitude of the K band spiral arms, and hence decrease the effect that we find. Hence our result is conservative and can only be strengthened by making corrections for this.

### 3.2 The effect of arm strength on star formation rate

A further test of arm induced star formation is to see if the enhancement in H$`\alpha `$ flux in the arm regions relates to the strength of shocks in the arms. Seigar & James (1998b) used the ratio between arm equivalent angle (EA) and arm full width half maximum (FWHM) as measured from K band images, as a measure of the relative strength of shocks in spiral arms. From here on, we refer to this ratio as the ‘central arm contrast’. Arm EA is defined as that angle subtended by the disc that contains an amount of light equivalent to that in the spiral arm (see Seigar & James 1998a for a detailed discussion) and is therefore a measure of the amount of light contained in a spiral arm. However, it contains no information about how concentrated the light from the arm is, and this is the justification in dividing this quantity by arm FWHM (measured from the cross–sectional light profiles of K band arms) to estimate the ‘central arm concentration’. If shocks are necessary to drive star formation, it seems reasonable that an arm with a high central arm concentration will be more efficient at forming stars than an arm with a low central arm concentration. We have investigated the relationship between this quantity and the H$`\alpha `$/K ratios measured in Table 2. The result is shown in Figure 1. One galaxy is left out of this plot, UGC 5786, which is an HII galaxy, and may have processes other than arm induced star formation, contributing to its extraordinarily high star formation rate. Figure 1 shows a non–linear relationship between these two quantities, which is significant at the 99.93% confidence level. It therefore seems that the star formation in spiral galaxies is triggered by shocks in arms, up to a certain limiting threshold, where H$`\alpha `$/K$``$2.

### 3.3 The star formation rate in UGC 5786

Finally, UGC 5786 (NGC 3310) has a strong ring of star formation at an average radius of 0.42$`\pm `$0.06 kpc from its centre. The ring of star formation observed in this galaxy has been the subject of many previous studies (e.g. Elmegreen et al. 2002). Rings of nuclear star formation are typically associated with the Inner Lindblad Resonance (ILR) and have been found in other galaxies, e.g. M 74 (Wakker & Adler 1995; James & Seigar 1999). It should also be noted that UGC 5786 is classified as an HII galaxy. Bell & Kennicutt (2001) report an H$`\alpha `$ luminosity, $`L_{H\alpha }=3\times 10^{41}`$ ergs s<sup>-1</sup>. Using the following equation,

$$SFR(M_{}yr^1)=7.9\times 10^{42}L_{H\alpha }(ergss^1)$$
(1)
they calculated a SFR of 2.37 M$`{}_{}{}^{}/`$yr. We have measured an $`L_{H\alpha }=`$(1.336$`\pm `$0.084)$`\times `$10<sup>42</sup> ergs s<sup>-1</sup>. Also using equation 1, we have calculated a SFR of 10.55$`\pm `$0.66 M$`{}_{}{}^{}/`$yr. This is clearly in disagreement with the SFR reported in Bell & Kennicutt (2001). It should be noted that we calculated a distance of 18.2 Mpc from a Virgocentric inflow model (and using H<sub>0</sub>=75 kms<sup>-1</sup>Mpc<sup>-1</sup>) for UGC 5786, whereas Bell & Kennicutt (2001) use the distance calculated by Waller et al. (1997), which is 13.9 Mpc. However, even taking into account this difference in measured distances to UGC 5786, our SFR is still a factor of $``$2.6 larger than that reported by Bell & Kennicutt (2001).

## 4 CONCLUSIONS

We have found an increase in H$`\alpha `$ flux in the vicinity of K band arms in this sample. If the K band arms are dominated by light from the old stellar population, then this can be interpreted as star formation triggered by a density wave. This result is in agreement with Lord & Young (1990), who found a higher star formation efficiency in arm regions. It also agrees with Cepa & Beckman (1990) and Knapen et al. (1992) who found evidence for spiral arm triggering of star formation via a non–linear dependence of SFR on gas density. This result therefore agrees with the large scale shock scenario (Roberts 1969; Shu et al. 1972; Tosa 1973; Woodward 1975; Nelson & Matsuda 1977), which predicts a higher star formation efficiency in arm regions. It does not support the ideas of stochastic self–propagation star formation (e.g. Gerola & Seiden 1978; Seiden & Gerola 1982; Seiden 1983; Jungwiert & Palous 1994; Sleath & Alexander 1995, 1996) which predict no change in the star formation efficiency from interarm to arm regions.

## 5 ACKNOWLEDGMENTS

The United Kingdom Infrared Telescope is operated by the Joint Astronomy Centre on behalf of the U.K. Particle Physics and Astronomy Research Council. The Jacobus Kapteyn Telescope and the William Herschel Telescope are operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The authors thank Neville Shane for both helping out with the observations and calculations provided in this paper. We also thank the referee for comments which improved the content of this paper.


text:
# Vortex Glass and Vortex Liquid in Oscillatory Media

\[

## Abstract

We study the disordered, multi-spiral solutions of two-dimensional homogeneous oscillatory media for parameter values at which the single spiral/vortex solution is fully stable. In the framework of the complex Ginzburg-Landau (CGLE) equation, we show that these states, heretofore believed to be static, actually evolve on ultra-slow timescales. This is achieved via a reduction of the CGLE to the evolution of the sole vortex position and phase coordinates. This true defect-mediated turbulence occurs in two distinct phases, a vortex liquid characterized by normal diffusion of individual spirals, and a slowly relaxing, intermittent, “vortex glass”.

\]

Spiral waves are ubiquitous in oscillatory and excitable two-dimensional active media . Their cores are robust wave sources which determine the oscillating frequency of the entire system and may dominate the surrounding dynamics. Spirals (often called vortices) may spontaneously appear and annihilate in a typical manifestation of spatiotemporal chaos. When the single spiral solution is stable, one easily observes locked, quasi-frozen, multi-spiral disordered structures whose glassy character has been suggested . Despite being found in most models of excitable or oscillatory media as well as in experiments (Belousov-Zhabotinskii reaction , surface growth , and others) surprisingly very little is known about the properties of these disordered states and the transitions leading to their formation.

In this Letter we investigate, in the framework of the complex Ginzburg-Landau equation (CGLE), the general issue of these multi-spiral solutions and show that they actually evolve on ultra-slow timescales. This is achieved thanks both to long numerical simulations of the CGLE and to a quantitatively correct reduction of the dynamics to the evolution of the sole vortex position and phase coordinates. This true defect-mediated turbulence occurs in two distinct phases, a vortex liquid characterized by normal diffusion of individual spirals, and a slowly relaxing, intermittent, “vortex glass”.

The CGLE describes most properties of generic oscillatory media, at least at a qualitative level, even if one is not in the vicinity of a supercritical long-wavelength Hopf bifurcation, where it can be systematically derived (for reviews, see ). Under appropriate scaling of the physical variables, it takes the universal form

$$_tA=A(1+ic)|A|^2A+(1+ib)\mathrm{\Delta }A,$$
(1)
where $`A`$ is a complex amplitude, $`b`$ and $`c`$ are real parameters characterizing relative dispersion and nonlinear frequency shift, and $`\mathrm{\Delta }`$ is the Laplace operator. Intensive studies conducted over the last ten years revealed a wide variety of striking dynamical phenomena in one, two, and three space dimensions, many of which were also observed in various experimental contexts, sometimes up to a quantitative agreement with CGLE predictions.

A distinctive feature of the two-dimensional CGLE is the existence of nontrivial sources of spiral waves (vortices) which determine the oscillating frequency of the entire system . The single-spiral solution for a vortex centered at $`𝐫_0`$ reads as: :

$$A_\mathrm{s}(𝐫,t)=F(r)\mathrm{exp}[i(m\theta +\psi (r)+\omega t+\varphi )],$$
(2)
where $`r=|𝐫𝐫_0|`$, $`\theta `$ is the polar angle measured from the vortex core, $`\varphi `$ is an arbitrary phase, and $`m=\pm 1`$ is the topological charge. Far away from the core, the solution approaches a plane wave with $`\psi (r)kr`$, where the asymptotic wavenumber $`k`$ is related to the rotation frequency as $`\omega =c(bc)k^2`$. The dependence of $`k`$ on $`b`$ and $`c`$ is known analytically for $`|bc|1`$ and $`|bc|1`$, e.g. $`kc^1\mathrm{exp}(\pi /|2c|)`$ for $`b=0`$ and $`|c|1`$ . In the parameter region where the single-spiral solution is stable and $`cb`$, the interaction between two well-separated spirals falls off exponentially . The screening was attributed to the shock lines where the waves emitted by the cores collide. For relatively small $`b,c`$ the interaction is monotonic, and the spirals exhibit asymptotic repulsion irrespectively of their charge (“monotonic range”) . In contrast, for larger $`b,c`$ satisfying the condition $`(cb)/(1+bc)>c^{}0.845`$ the velocity vs distance dependence is modulated (“oscillatory range”), and the spirals become keen to form a variety of long-living (but unstable) bound states .

While the above interactions cannot account for the strongly chaotic regimes where many defects are spontaneously generated and undergo violent motion, they are expected to play a leading role in the occurrence of the quasi-static structures of large and small spirals surrounded by a complex network of shocks commonly observed in the large region of the $`(b,c)`$ parameter space where the single spiral solution is stable . Rather little is known about these cellular structures except that they seem to be frozen on usual observation timescales. Our results lead to a dramatically different picture.

We conducted detailed numerical studies of the two-dimensional CGLE restricted to the case $`b=0`$ and $`c>0`$ (for the regimes of interest scaling relations apply in parameter space ). The integration domain with periodic boundary conditions was typically of area $`S=512\times 512`$ and the integration time about $`10^7`$. Various integration schemes were used and all results checked against variation of the numerical resolution. We monitored the positions $`𝐫_j(t)`$ of all $`N`$ spirals. We measured their instant “activity” $`T=S^1_S𝑑𝐫\left|_t|A|\right|`$ and —when applicable— the spiral diffusion coefficient $`D=\frac{1}{Nt}_{j=1}^N|𝐫_j(t)𝐫_j(0)|^2`$. Since, for well-separated spirals, $`|A|`$ is constant everywhere except near the cores, $`T`$ is related to the velocity of spirals.

In the monotonic range $`c<c^{}`$ (Fig. 1), after a short transient during which the initial number of defects may decrease, the number of spirals remains constant and $`T`$ fluctuates around a well-defined mean value. The typical vortex velocity is very small (of the order of $`10^4`$ for $`c=0.7`$), and the velocity distribution is approximately Maxwellian (Fig. 1d). In large enough systems, the spiral cores perform normal diffusive motion (Fig. 1d). This vortex liquid can be characterized by a viscosity $`\nu D^1`$. Since all scales, already very large near $`c^{}`$, diverge exponentially when $`c0`$, an extensive investigation of the dependence of $`D`$ on $`c`$ and on the vortex density $`\rho N/S`$ is beyond reach of present-day computers. Nevertheless, our data indicate that $`D`$ increases with $`\rho `$.

In contrast, in the oscillatory range $`c>c^{}`$, after a very long, possibly infinite, transient (Fig. 2e), the population of spirals spontaneously segregates into two distinct phases: large and almost immobile spirals and droplets (clusters) of small vortices confined between them (Fig. 2a,b, see also ). One can thus define a “liquid fraction” of small spirals, whose sizes are typical of the vortex liquid observed in the monotonic range. When the liquid fraction is small, the resulting state exhibits slow intermittent dynamics (bursts of activity separated by long quiescent intervals, Fig. 2d) reminiscent of glassy dynamics, whereas for higher fractions the activity $`T`$ fluctuates around some mean value and the “liquid” vortices can be shown to perform normal diffusion in the space surrounding the big spirals. This picture is reminiscent of phase coexistence in first-order transitions of isolated equilibrium systems . In this respect, our results indicate that one can indeed distinguish a “vortex glass” phase. Unfortunately, because of the ever slower timescales over which the spatial structure of the system evolves, a precise characterization of this phase using simulations of the CGLE is currently impossible.

To obtain further insight into the problem, we reduced the CGLE to a set of ordinary differential equations describing the motion of the cores $`𝐫_j(t)`$ and the phases $`\varphi _j(t)`$ of individual spirals. Only pair interactions are taken into account, making use of the results of . The original partial differential equation is thus replaced by a set of $`3N`$ first-order equations governing the core positions $`𝐫_j(t)`$ and the spiral phases $`\varphi _j(t)`$. The phases determine the configuration of the shock lines, which, in turn, yield the velocity of the cores. One clear advantage of this approach, in addition to virtually suppress physical space, is that it allows to bypass the ultraslow timescale related to the vanishing value of the wavenumber $`k`$ as $`c0`$.

The problem of the interaction between two oppositely-charged spirals is equivalent to that of one spiral and a straight shock line, i.e. a half-plane with the no-flux boundary condition $`A/n=0`$. For $`c1`$, the equations governing the position and phase of the spiral are

$`{\displaystyle \frac{dz}{dt}}`$ $`=`$ $`2c^2\pi k(imcB^{}){\displaystyle \frac{\mathrm{exp}(2|ck|X)}{\sqrt{\pi |ck|X}}},`$ (3)
$`\dot{\varphi }`$ $`=`$ $`2ck^2\pi {\displaystyle \frac{\mathrm{exp}(2|ck|X)}{\sqrt{\pi |ck|X}}}`$ (4)
where $`z=x+iy`$ is the (complex) spiral position, $`dz/dt=v_x+iv_y`$, $`\dot{\varphi }`$ is the phase correction, $`m`$ is the topological charge, $`B^{}0.48`$, and $`X`$ is the distance from the core (located at $`(X,0)`$) to the boundary (the y-axis is chosen to be on the shock line) . Eqs. (4) imply the asymptotic repulsion of the spiral from the boundary, in agreement with simulations . This prevents spiral annihilation in the vortex liquid.

The equations in the oscillatory range $`c>c^{}`$ read:

$`C_xv_x+mC_yv_y={\displaystyle \frac{k\sqrt{1k^2}\mathrm{exp}(pX)}{\delta \sqrt{2\pi pX}}}X^\mu ,`$ (5)
$`\dot{\varphi }`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{Im}(C_{10}C_0^{})}}\mathrm{Im}\left({\displaystyle \frac{k\sqrt{1k^2}\mathrm{exp}(pX)}{\delta \sqrt{2\pi pX}}}X^\mu C_0^{}\right)`$ (6)
where the complex constants $`C_{x,y},C_0,C_{10}`$ are obtained numerically and the real parameters $`p,\mu ,\delta `$ are derived from the linearized CGLE . Eqs. (5) describe stable bound states of spirals moving along the plane boundary. However, the problem is complicated by the fact that the position of the shock line (where the waves emitted from the cores collide) depends on the relative phase of the spirals. From the condition that the total spiral phases $`\mathrm{\Phi }_j=\psi (r)+\varphi _j|k|r+\varphi _j`$ are equal, one finds that the distance to the shock is given by:

$$X=|𝐫_k𝐫_j|/2(\varphi _k\varphi _j)/2|k|$$
(7)
Substituting (7) into Eqs. (5) one finds that the symmetric bound states are unstable in the oscillatory range . Similar phenomena occurs for likely-charged spirals.

Eqs. (4), (5) and condition (7) were used to investigate the dynamics of many-spiral states. In this case, one considers each spiral in the local coordinate system associated with every other spiral and sum up their contributions (details will be published elsewhere). Substituting $`r2c|k|r,t4\sqrt{2\pi }c^3k^2`$ and introducing a rescaled “phase length” $`\zeta =c\varphi `$, one obtains from Eqs.(4):

$`{\displaystyle \frac{dz_j}{dt}}`$ $`=`$ $`{\displaystyle \underset{k}{}}(cB^{}+im_k){\displaystyle \frac{z_kz_j}{|z_kz_j|}}{\displaystyle \frac{\mathrm{exp}(X_k)}{\sqrt{X_k}}},`$ (8)
$`\dot{\zeta _j}`$ $`=`$ $`{\displaystyle \underset{k}{}}{\displaystyle \frac{\mathrm{exp}(X_k)}{2c\sqrt{X_k}}}`$ (9)
where $`X_k=|z_kz_j|/2+\zeta _j\zeta _k`$. A similar rescaling was performed in the oscillatory case.

We investigated numerically and analytically the reduced equations which, in fact, are rather different in the oscillatory range  and in the monotonic range . Overall quantitative agreement with the full CGLE dynamics has been found, which constitutes maybe the first success for the ideas underlying the concept of “defect-mediated turbulence” put forward in the 80’s .

In particular, we found, in the monotonic range, clear “liquid” behavior accompanied by normal diffusion for all $`c`$ values. Each spiral core moves chaotically and the velocity distribution is also very close to a Maxwellian. The reason underlying the chaotic behavior of this vortex liquid is that Eqs. (9) do not obey a variational principle for any $`c`$ value, due to the nontrivial form of the pairwise interaction. However, an ensemble of likely-charged spirals tend to form a stable hexagonal “Wigner” crystal due to mutual repulsion. With periodic boundary conditions (used here), one has an equal mixture of positive and negative spirals, and it can be proved that a square lattice of spirals with alternating charges is unstable in a large system with respect to long-wavelength perturbations. This explains the short-range crystalline order and local spatial charge separation observed (Fig. 1c). The faithfulness of the reduced equations to the CGLE is also testified by the quantitative agreement: found for the value of the diffusion constant $`D`$. for $`c=0.6`$, the domain size $`L=512`$ and vortex number $`N=46`$ the diffusion $`D`$ obtained from CGLE is $`D0.0033`$ and from Eqs. (9) is $`D0.0036`$. The reduced equations do allow for an extensive study of the asymptotic (i.e. large-size) properties of the vortex liquid. In particular, we find that $`D\mathrm{exp}(2c|k|/\sqrt{\rho })`$ at fixed $`c`$ value (Fig. 1f). Taking into account that the typical inter-spiral distance varies like $`1/\sqrt{\rho }`$, this scaling law is in agreement with the expression for the spiral velocity as given by Eq.(4).

In the oscillatory range also, the reduced equations faithfully reproduce the phenomenology of symmetry breaking, intermittent activity at low vortex densities, and “vortex glass” formation observed in the CGLE (Fig. 2c,f). The reduced equations provide an interesting framework to discuss the mechanisms leading to intermittent dynamics in the vortex glass containing liquid droplets. For a spiral embedded in a large domain of arbitrary shape, Eqs. (5) generally do not have a stationary solution $`𝐯_j=\dot{\varphi }_j=0`$ for all $`j`$ (because the complex constants $`C_{x,y}`$ and $`C_{10,0}`$ are typically not proportional). As a result, it is impossible to satisfy the conditions $`𝐯_j=0`$ and $`\dot{\varphi }_j=0`$ simultaneously. The drift of the spirals changes the phases and, therefore, the shape of their domains. Eventually, when the phase difference becomes large enough, the shock line boundaries rearrange and trigger rapid dynamics within the liquid droplets. Since the velocity of the cores is an exponential function of the inter-spiral separation, the typical timescale between such rapid rearrangement events is exponentially large in the spiral size. Thus, in very large systems, where the domains occupied by spirals may be comparable with system size, the dynamics will be very slow indeed.

We have found that the transition from vortex liquid to vortex glass possesses some features of a first order phase transition. Additional numerical studies will be needed to precise the structure of corresponding phase diagram. Returning to the question of the possibility of “true” glassy dynamics  (as testified, e.g., by some aging phenomena), several comments are now in order. Preliminary numerical results of the reduced equations in very large systems at low vortex density indicate that the intermittent space-time activity usually reaches a complex but self-averaging asymptotic regime, in contrast with spin-glass-type behavior  and in closer agreement with a structural glass, which can be viewed sometimes as a very viscous fluid. On the other hand, these regimes, as far as they can be studied for the full CGLE, do exhibit a slow, possibly aging-like, decay of the number of vortices in the system (Fig. 2e). This leaves the door open for an actual vortex glass in this fully-deterministic, noiseless, disorder-free context — a challenge to models of glassy behavior in statistical physics.

Our work has shown that the multi-spiral “frozen” states of the CGLE actually evolve on very long timescales. These ultra-slow regimes are relevant, we believe, to the many related models and experimental situations. A vortex liquid and a vortex glass can be distinguished depending on the effective interaction law between spirals. These “phases” are dynamical, and are remarkably well accounted for, in the case of the CGLE, by a reduction of the full problem to a finite set of ordinary differential equations governing the vortex cores. Given that the typical timescales of vortex motion are several orders of magnitude larger than the period of the basic oscillations, any experimental observation will require a very good control of the apparatus over long times. This work was supported by the U.S.Department of Energy under contract W-31-109-ENG-38.