text:
# Excited States in Staggered Meson Propagators

## 1 INTRODUCTION

Extracting ground state hadron masses from lattice operator correlation functions has become a standard exercise, and results generally agree with experimental measurements. (See, for example .) However, determining excited state masses on the lattice (e.g. ) is more difficult. Such determinations can provide tests of the particular lattice approach, can be used to study decays of unstable hadrons and, in principle, could help clarify quark model interpretations of higher mass hadrons.

## 2 SIMULATION

We measured meson propagators on lattices generated using the “Asqtad” action with $`2+1`$ flavors of dynamical quarks, that is, one flavor with mass approximately equal to the physical strange quark mass and two lighter flavors with degenerate masses. (See Table 1.) We describe the $`20^3\times 64`$ lattices as “coarse” lattices with lattice spacing of $`a0.12`$ fm, while the $`28^3\times 96`$ “fine” lattices have a spacing of $`a0.09`$ fm. More detailed descriptions of these simulations and analysis procedures can be found in Refs. . In this report we discuss excited states in the pseudoscalar ($`0^+`$) and scalar ($`0^{++}`$) isovector meson propagators.

## 3 RESULTS AND CONCLUSIONS

### 3.1 $`0^+`$ states

In the staggered fermion formulation propagators generally appear with an oscillating component that represents a parity partner state. There is no partner to the taste-pseudoscalar pion and $`s\overline{s}`$ states, hence these propagators contain contributions only from $`J^{PC}=0^+`$ states. Moreover, the signal to noise ratio is large at all time separations, so it is relatively easy to measure the contribution of excited states.

We fit $`0^+`$ propagators to the form:

$`P(t)`$ $`=A_0(e^{m_0t}+e^{m_0(Tt)})`$
$`+A_1(e^{m_1t}+e^{m_1(Tt)}),`$
where $`T`$ is the length of the time dimension of the lattice, $`t`$ is the time separation of the source and sink operators and $`m_0`$ and $`m_1`$ are the masses of the ground state and lowest excited states. We look for high-confidence, low-error fits that fall in plateaus of stability of the fit parameters with respect to the limits $`t_{\mathrm{max}}`$ and $`t_{\mathrm{min}}`$ of the propagator included in the fit. Figure 1 shows a set of pion fits with varying $`t_{\mathrm{min}}`$, and $`t_{\mathrm{max}}`$ fixed at 28.

Figure 2 shows the results of this fitting procedure for all of the data sets. The darker error bars are statistical, and lighter error bars represent systematic error due to fit choice. We plot the product of ground state and first excited state masses with $`r_1`$, a length scale derived from the static quark potential, against the ratio $`(m_\pi /m_\rho )^2`$, a measure of the quark mass. At the physical value of $`(m_\pi /m_\rho )^2=0.03`$ we plot the experimental mass of the states $`\pi ^0(135)`$ and the $`\pi (1300)`$ scaled with $`r_1=0.317(7)`$ fm. Also, an excited $`0^+`$ $`s\overline{s}`$ state, an average of measurements on all of the fine lattices, is compared to the experimental mass of the $`\eta (1440)`$, with which agreement is very marginal.

We feel that the consistency of the simulation results with the $`\pi (1300)`$ is an important indication that the fourth root of the determinant in the action does not introduce a non-locality with scale larger than the $`\pi (1300)`$ Compton wavelength: 0.15 fm.

### 3.2 $`0^{++}`$ states

The $`a_0`$ meson ($`J^{PC}=0^{++}`$) propagator occurs as the oscillating parity partner in the taste-scalar pion propagator. With states of both parities contributing to the propagator, it is more difficult to resolve excited states. We fit these $`a_0\pi `$ propagators to the form:

$`P(t)`$ $`=`$ $`A_0(e^{m_0t}+e^{m_0(Tt)})`$
$`+`$ $`A_1(e^{m_1t}+e^{m_1(Tt)})`$
$`+`$ $`A_2(1)^t(e^{m_2t}+e^{m_2(Tt)})`$
$`+`$ $`A_3(1)^t(e^{m_3t}+e^{m_3(Tt)}).`$
Here $`m_2`$ and $`m_3`$ are the $`0^{++}`$ masses. In one case, $`\beta =7.09`$, we were able to unambiguously discern an excited $`0^{++}`$ state. See Figure 3.

The $`0^{++}`$ states are interesting in that several of the lowest lying states are likely not pure $`q\overline{q}`$ states because the two meson states to which the $`a_0`$ decays are accessible in our simulations. In Figure 4, the straight line is an extrapolation of the heavier-$`m_q`$ quenched $`a_0`$ masses, and is an expectation of the mass of a $`q\overline{q}`$ state with $`0^{++}`$. The curved line, which crosses below the “$`a_0`$-line” at low $`m_q`$ traces the sum of the $`\pi `$ and $`\eta `$ ground state masses. The fitted full-QCD $`0^{++}`$ ground state masses roughly trace the $`\pi +\eta `$ line, possibly displaying the lower half of an avoided level crossing, an indication of the decay of the pure $`q\overline{q}`$ $`a_0`$. The measured $`0^{++}`$ excited state (burst) is not consistent with the straight line either, but relatively close to the experimental value of the mass of a $`K\overline{K}`$ molecule. Future efforts will try to fill out Figure 4 with more $`0^{++}`$ excited states.


format: text

text:
# KEK-CP-145, UT-CCP-P139, UTHEP-473 Study of finite volume effects in the non-perturbative determination of 𝑐_sw with the SF method in full three-flavor lattice QCDPresented by K-I.I.

## 1 Introduction

Last year we reported the existence of a sizable finite volume effect in the non-perturbative $`c_{\mathrm{sw}}`$ determined by the Schrödinger functional (SF) method at a fixed lattice size ($`L/a`$$`=`$const) when the RG-improved gauge action is adopted in full three-flavor QCD . This effect causes a constant deviation from the true non-perturbative value of $`c_{\mathrm{sw}}`$, and yields $`O(a)`$ errors in physical observables.

This year we explain our strategy for removing the unwanted finite volume effect from $`c_{\mathrm{sw}}`$ and report the bare coupling dependence of the non-perturbative $`c_{\mathrm{sw}}`$ for three-flavor dynamical QCD.

## 2 Method and Simulations

The basic strategy to obtain the non-perturbative $`c_{\mathrm{sw}}`$ with the SF setup is described in . We used Eq. (2.17) of  for the SF boundary coupling of the RG-improved gauge action. The polynomial Hybrid Monte Carlo (PHMC) algorithm  is employed for the numerical simulation of the three-flavor dynamical QCD.

To remove the finite lattice size effect we fix the physical lattice extent $`L`$ so that $`a/L`$ corrections vanish with $`a`$. This can be achieved by the following steps.
(1) Define the reference lattice size $`L^{}/a`$ at a bare coupling value $`g_{0}^{2}{}_{}{}^{}`$. This also defines the reference physical lattice extent $`L^{}`$.
(2) For other $`g_0^2`$’s, the lattice sizes $`L^{}/a`$ at each $`g_0^2`$ are estimated by the beta function.
(3) Determine $`c_{\mathrm{sw}}`$ and $`\kappa _c`$ via the PCAC condition with the SF setup at each $`g_0^2`$ on the corresponding lattice size $`L^{}/a`$.

In actual numerical simulations we employ the universal two-loop beta function to estimate $`L^{}/a`$. This causes a systematic error in the non-perturbative $`c_{\mathrm{sw}}`$, which could be large in the strong coupling region. This systematic error, however, can be partly avoided by defining $`L^{}`$ at the largest coupling value. We can also examine the magnitude of the error by investigating the $`L/a`$ dependences because the error only appears through the incorrect value of $`L^{}/a`$.

With the above consideration we take the following setup in the simulations.
(i) We define $`L^{}`$ so that $`L^{}/a=6`$ holds at $`\beta =6/g_0^2=1.9`$.
(ii) We employ both $`L/a=8`$ and $`6`$ at $`\beta =2.0`$, and $`L/a=8`$ for other weaker couplings.
(iii) We correct $`c_{\mathrm{sw}}`$ and $`\kappa _c`$ at $`L/a`$ to those at $`L^{}/a`$ using the one-loop correction calculated with the SF setup in finite lattice volume .
Although the perturbative error from setup (iii) appearers in strong coupling region, the magnitude can be reduced when $`L^{}/aL/a`$ holds and $`L/a`$ dependence is small. The small $`L/a`$ dependence is confirmed at $`\beta =2.0`$ as shown later. We will briefly discuss the possible resolution of the systematic error from the use of perturbative beta function in the last section.

We use the following non-perturbative tuning condition for $`c_{\mathrm{sw}}`$ and $`\kappa `$.

$$\begin{array}{c}M(g_0^2,L/a,\kappa ,c_{\mathrm{sw}})=0,\hfill \\ \mathrm{\Delta }M(g_0^2,L/a,\kappa ,c_{\mathrm{sw}})=0,\hfill \end{array}$$
(1)
where $`M(g_0^2,L/a,\kappa ,c_{\mathrm{sw}})`$ is the bare PCAC quark mass and $`\mathrm{\Delta }M(g_0^2,L/a,\kappa ,c_{\mathrm{sw}})`$ is the difference of the two bare PCAC quark masses with different definition. The actual definitions are explained in . The difference from the tuning condition we employed last year is the absence of the tree-level finite PCAC quark mass which arises from the SF boundary for a finite lattice extent. The tree-level finite volume effect is tuned to the corresponding value at $`L^{}/a`$ by setup (iii) together with the one-loop finite volume effect.

As shown in setup (iii), we use the following one-loop perturbative correction to convert $`c_{\mathrm{sw}}`$ and $`\kappa `$ at $`L/a`$ to those at $`L^{}/a`$ at each $`g_0^2`$;

$$\begin{array}{c}\kappa _c(L^{}/a)=\kappa _c(L/a)+\delta \kappa _c(L/a;L^{}/a),\hfill \\ c_{\mathrm{sw}}(L^{}/a)=c_{\mathrm{sw}}(L/a)+\delta c_{\mathrm{sw}}(L/a;L^{}/a).\hfill \end{array}$$
(2)
$`\delta \kappa _c(L/a;L^{}/a)`$ and $`\delta c_{\mathrm{sw}}(L/a;L^{}/a)`$ are defined as

$$\begin{array}{c}\delta \kappa _c(L/a;L^{}/a)=\left[\kappa _c(L^{}/a)\kappa _c(L/a)\right]_{\text{1-loop}},\hfill \\ \delta c_{\mathrm{sw}}(L/a;L^{}/a)=\left[c_{\mathrm{sw}}(L^{}/a)c_{\mathrm{sw}}(L/a)\right]_{\text{1-loop}},\hfill \end{array}$$
(3)
where the right hand side of Eq. (3) is calculated up to one-loop level.

## 3 Results

We have carried out the numerical simulation in the region $`\beta =12.01.9`$ with the setup explained in Sec. 2. It was found that below $`\beta =2.2`$ the PHMC algorithm failed at vanishing or negative PCAC quark masses due to large quantum fluctuations. We extrapolated them from $`M>0`$ to $`M=0`$ as shown in Fig. 1. It was also found that below $`\beta =1.9`$ with $`L/a=8`$ it became difficult to satisfy the condition $`\mathrm{\Delta }M=0`$, which is why we employ $`L^{}/a=6`$ at $`\beta =1.9`$ for the definition of the reference physical lattice extent.

Figures 2 and 3 are the result for $`c_{\mathrm{sw}}`$ and $`\kappa _c`$. Open circles are the result without $`L/a`$ matching. Triangles are after tuning to $`L^{}/a`$ by the one-loop correction (or interpolation at $`g_0^2=3.0`$). For $`c_{\mathrm{sw}}`$ we observe large one-loop corrections in the region $`1.0<g_0^2<2.7`$. This is due to rather strong coupling and a large discrepancy between $`L/a`$ and $`L^{}/a`$. We discard the data from the functional fit. At strong couplings, $`g_0^2=6/2.2`$ and $`6/2.1`$, the one-loop correction is rather small since $`L/aL^{}/a`$. The $`L/a`$ dependence at $`g_0^2=3.0`$ $`(\beta =2.0)`$ is negligible compared to the statistical error (upper circle: $`L/a=8`$, lower circle: $`L/a=6`$). We have $`6<L^{}/a<8`$ via the two-loop beta-function at $`g_0^2=3.0`$, we simply interpolate $`c_{\mathrm{sw}}`$ to $`L^{}/a`$ using both lattice size data. While uncorrected data undershoots the one-loop result for $`L/a=\mathrm{}`$ (dotted line) in the weak coupling region , those with one-loop correction gradually approaches it ($`L^{}/a`$ is essentially infinite in the region). For $`\kappa _c`$ the one-loop correction is rather small. We conclude that the non-perturbative results in the weak coupling region is consistent with the one-loop one for both $`c_{\mathrm{sw}}`$ and $`\kappa _c`$. We fit the filled triangles to obtain the functional form in $`g_0^2`$ for $`c_{\mathrm{sw}}`$ and $`\kappa _c`$.

After the curve fitting including the knowledge of the one-loop perturbative coefficient we obtain

$$\begin{array}{c}c_{\mathrm{sw}}(g_0^2)=1+0.113g_0^2+0.0209(72)g_0^4\hfill \\ +\mathrm{\hspace{0.17em}0.0047}(27)g_0^6,\hfill \end{array}$$
(4)
with $`\chi ^2/\mathrm{DOF}=0.5754`$, and

$$\begin{array}{c}\kappa _c(g_0^2)=0.125+\left(3.681192\times 10^3\right)g_0^2\hfill \\ +\left(0.141(43)\times 10^3\right)g_0^4\hfill \\ +\left(0.124(67)\times 10^3\right)g_0^6\hfill \\ \left(0.049(21)\times 10^3\right)g_0^8,\hfill \end{array}$$
(5)
with $`\chi ^2/\mathrm{DOF}=1.0014`$. All data are statistically independent and $`\chi ^2/\mathrm{DOF}`$ results in reasonable values. We show statistical errors only. The systematic error from the use of perturbative beta function is not estimated. We, however, consider that the systematic error is minimized since the slope in $`g_0^2=3.03.15`$ of $`c_{\mathrm{sw}}`$ should be correct because of the negligible $`L/a`$ dependence at $`g_0^2=3.0`$.

## 4 Summary

We have extracted the non-perturbative $`c_{\mathrm{sw}}`$ as a function of $`g_0^2`$ for $`N_f=3`$ with the RG-improved gauge action. The functional form is shown in Eq. (4). In deriving the non-perturbative values at a fixed physical lattice extent $`L^{}`$, we used the two-loop beta-function to obtain the corresponding lattice size $`L^{}/a`$ at a given $`g_0^2`$. This may contradict the word “non-perturbative”. In order to improve this point we should iteratively execute our procedure with the non-perturbative beta-function that is measured with the non-perturbative $`c_{\mathrm{sw}}`$ extracted from the previous iteration. We leave this iterative study for future studies. The first step to obtain $`N_f=2+1`$ hadron mass spectrum has been started using the non-perturbative $`c_{\mathrm{sw}}`$ of Eq. (4.

This work is supported by Large Scale Simulation Program No. 98 (FY2003) of High Energy Accelerator Research Organization (KEK), and also in part by Grants-in-Aid of the Ministry of Education (Nos. 13135204, 14046202, 15204015, 15540251, 15540279). N.Y. is supported by the JSPS Research Fellowship.


text:
# The spatial profile of the matter distribution in a dilute atomic-molecular Fermi cloud

## Abstract

A dilute atomic-molecular fermion in a trap acquires a rather complex layered spatial structure. At very low temperatures, typically at the center of the trap, a superfluid bosonic core is formed out of dimers, weakly bound pairs of fermionic atoms. Surrounding this core is a mantle, formed by a mixture of superfluid bosonic dimers and superfluid fermionic atoms. The outermost layer consists of only fermionic atoms, forming mostly a normal Fermi gas.

Jin’s group  has demonstrated that weakly bound dimers of fermionic atoms are easily created in a cloud of cold Fermi atoms, by ramping up and down the magnetic field near a Feshbach resonance. Several other groups have been able to quickly confirm these results, by using similar experimental techniques . In all these experiments a relatively cold Fermi gas of atoms, interacting with a negative scattering length, was driven past the Feshbach resonance into a mixture of fermionic atoms and bosonic dimers, all of them interacting with positive scattering lengths . The ratio of atoms to molecules in these experiments could be varied essentially at will, by tuning the ramping rate of the ambient magnetic field. It was also explicitly demonstrated that this mixed system of fermionic atoms and bosonic dimers has a rather long lifetime  and that the rate of various recombination processes is rather low, an effect qualitatively suggested in Ref. . In a series of related experiments it was shown that the weakly bound dimers form a Bose-Einstein condensate .

The natural question arises: ”What is the nature of this new quantum object created in these experiments?” As it was shown in Ref. , in the regime where the scattering length is large, but at the same time the various average inter-particle separations are even larger, it is possible to perform a rather straightforward theoretical analysis of its properties. If the system could be regarded as locally homogeneous, then by changing the scattering length from negative to positive values part of the system undergoes a BCS $``$ BEC crossover . At zero temperature, when the scattering length is negative, one would expect that the atomic Fermi system would be in a BCS state at sufficiently low temperatures. By changing adiabatically the scattering length to positive values, the size of the Cooper pairs would gradually decrease, until dimers, with sizes smaller than the average inter-particle separation, are eventually formed and the entire system turns into a pure BEC. At the finite temperatures achieved so far in atomic Fermi clouds, typically only a fraction of the atoms are converted into dimers and one ends up with a mixed system of fermions and bosons.

In all these experiments  the estimated final temperatures were well below the critical BEC temperature of the Bose component. One could have naively concluded that a novel mixture of a Bose superfluid mixed with a normal Fermi gas was thus realized. One reason to believe this would be the fact that the atomic scattering length is positive, which corresponds to an effective repulsion at these energies. Therefore the Fermi subsystem could be expected to be a normal Fermi system. The other reason would be the fact that the dimers were expected to repel each other  and thus they can form a BEC. This aspect was subsequently confirmed experimentally . A third reason would be the fact that atoms and dimers repel under these conditions too .

Nevertheless, a somewhat subtle quantum phenomenon could lead to rather dramatic changes of the properties of the fermionic component. Namely, even though two fermions effectively repel each other (since the scattering length is positive), there is a boson induced interaction, which is attractive and relatively long ranged. A fermion traveling in a Bose superfluid (fermions always move) will excite a Bogoliubov sound mode of the Bose subsystem, which in its turn will propagate to rather large distances, where it will ”shake” another fermion. This type of interaction between fermions is reminiscent of the phonon mediated electron-electron interaction in metals, which is responsible for the superconductivity. Technically, this phonon exchange between two fermions is a second order perturbation theory contribution to the ground state energy of the system, and as such is always attractive. As a result, if the temperatures are low enough and the induced attraction strong enough, a homogeneous mixture of such bosons and fermions could turn into a mixtures of a Bose superfluid in interaction with a Fermi superfluid. If the system is dilute and if all the scattering lengths are known  and the energy of such a system is easily calculated. Here, I shall primarily be concerned with the spatial profile and spatial structure of such a mixed atom-dimer system at low temperatures. Since the induced pairing correlations in the Fermi subsystem are in the weak coupling regime , the presence or absence of the pairing correlations plays essentially no role in determining the spatial distribution of such an atomic cloud.

The properties of these novel Bose-Fermi mixtures of quantum fluids could be affected by the shape of the trapping potential in a non-trivial manner. As it will be shown here, the trapping potential has a rather strong influence on the detailed spatial matter density profile and spatial composition of these new objects. Even though one can create either a pure Bose system or a pure Fermi system, however, such objects are likely an exception rather than the rule. Most of the time one would create an object, which at the center will be a pure Bose superfluid, followed by a layer of mixed Bose and Fermi superfluids, with an outer layer being a normal Fermi fluid. At the interface between the normal and the superfluid Fermi systems, another pure quantum phenomenon, known as Andreev reflection, leads to the formation of superfluid correlations in the outer normal Fermi fluid. These pairing correlations will be present in a layer of width comparable with the so called coherence length . When a fermion, impinging on a normal-superfluid interface from the normal part, undergoes an Andreev reflection, it is reflected essentially backward as a hole. This process leads to the appearance in a relatively narrow region of the normal Fermi fluid of pairing correlations, even though the atom-atom interaction is repulsive. Again, the presence of such pairing correlations will have an exceedingly small influence on the spatial distribution of the atoms .

In Ref.  it was shown that the energy density of a mixed fermion-dimer system in a trap is given by the following expression:

$`={\displaystyle \frac{3}{5}}{\displaystyle \frac{\mathrm{}^2k_F^2}{2m}}n_f+{\displaystyle \frac{1}{2}}U_{ff}n_f^2`$ (1)
$`+U_{fb}n_fn_b+{\displaystyle \frac{1}{2}}U_{bb}n_b^2+\epsilon _2n_b+Vn_f+2Vn_b,`$ (2)
where $`V`$ is the trap potential, $`n_b`$ is the dimer number density, $`n_f=k_F^3/3\pi ^2`$ is the total number density of both fermion species, assumed here to be present in equal amounts, $`k_F`$ is the Fermi wave vector and where

$`\epsilon _2={\displaystyle \frac{\mathrm{}^2}{ma^2}},U_{ff}=2{\displaystyle \frac{\pi \mathrm{}^2a}{m}}>0,`$ (3)
$`U_{fb}=3.537{\displaystyle \frac{\pi \mathrm{}^2a}{m}}>0,U_{bb}=1.2{\displaystyle \frac{\pi \mathrm{}^2a}{m}}>0`$ (4)
and $`a`$ is the atom-atom scattering length. The trapping potential experienced by a dimer has twice the value of the corresponding potential experienced by an atom at the same position. The position dependence of the number densities and trap potentials are not shown explicitly. Since various recombination rates are sufficiently small  one can consider that the number of fermions and bosons are conserved separately and one can account for the corresponding slow changes later-on when and if needed. As it was mentioned a few times above and it was shown in Ref. , the fermions, in the presence of superfluid bosons, become superfluid at sufficiently low temperatures. The gain in the energy density arising from the appearance of superfluid correlations could in principle be included in the above energy density. However, this gain is very modest and is not expected to lead to any major changes of the spatial profile of this system. This gain in energy is called condensation energy and per particle is given by $`3\mathrm{\Delta }^2/8\epsilon _F\epsilon _F`$, where $`\epsilon _F`$ is the Fermi energy.

In the case of a large number of atoms one can find the atom and dimer number densities with excellent accuracy in the Thomas-Fermi approximation , namely form the equations:

$`{\displaystyle \frac{\mathrm{}^2k_F^2}{2m}}+\left(U_{ff}{\displaystyle \frac{U_{fb}^2}{U_{bb}}}\right)n_f`$ (5)
$`=\mu _f{\displaystyle \frac{U_{fb}}{U_{bb}}}\mu _b\left(1{\displaystyle \frac{2U_{fb}}{U_{bb}}}\right)V,`$ (6)
$`n_b={\displaystyle \frac{\mu _b2VU_{fb}n_f}{U_{bb}}}.`$ (7)
I have absorbed in the definition of the chemical potential $`\mu _b`$ the energy of the dimer $`\epsilon _2`$. These two equations have been obtained by assuming that both $`n_{f,b}0`$. If either one of the number densities is vanishing, then one has to solve instead one of the following two equations:

$`{\displaystyle \frac{\mathrm{}^2k_F^2}{2m}}+U_{ff}n_f=\mu _fV,\mathrm{i}fn_f>0\mathrm{a}ndn_b0`$ (8)
$`n_b={\displaystyle \frac{\mu _b2V}{U_{bb}}},\mathrm{i}fn_b>0\mathrm{a}ndn_f0`$ (9)
As usual in the Thomas-Fermi approximation, in the region of space where densities are exponentially small (classically forbidden regions), the corresponding densities are set equal to zero. In principle one can rather straightforwardly improve on this approximation, but changes in distributions are quantitatively minor. Since $`12U_{fb}/U_{bb}5`$, see Eq. (6), the effective trap potential experienced by fermions changes its character in the presence of dimers and the fermions are pushed outward. In the spatial region where both $`n_f0`$ and $`n_b0`$ the fermion number density increases outwardly. By analyzing various possibilities one is quickly convinced that besides the rather rare cases, when there are either only fermions or only bosons, a finite cloud consisting of both Fermi atoms and Bose dimers has fermions only in the outside layer.

I shall describe now the profiles of the atom and dimer number densities, by constructing them from point to point. I assume that the trap potential is increasing monotonically from its ”center” outward. The ”center” of the trap is defined as the point where the trapping potential has a minimum. The shape of the potential could otherwise be arbitrary. I shall fix the origin of the coordinate system in the ”center” of the trap and construct the number densities along an arbitrary straight line. Let me denote by $`𝑹_3`$ the farthest point from origin along this line where the atom density vanishes, that is where $`n_f(𝑹_3)=0`$. At this point

$$\mu _f=V(𝑹_3),$$
(10)
and one can thus fix the value of the atom chemical potential. According to Eq. (10) the atom number density vanishes along an equipotential surface. Starting from the point $`𝑹_3`$ one moves inward and point by point solves Eq. (8) and thus determine the atom number density along this chosen line. Now I shall assume that there are dimers present in the system as well, and that their number density vanishes along the same line at the farthest point from origin $`𝑹_2`$ ($`R_2<R_3`$). According to Eq. (7) at this point the following condition is fulfilled

$$\mu _b=2V(𝑹_2)+U_{fb}n_f(𝑹_2),$$
(11)
which thus determines the value of the dimer chemical potential. Again, the locus of points where $`n_b(𝑹_2)=0`$ is another equipotential surface of the trapping potential. From this point inward, one has now to solve the equations (6) and (7). One can then find a point $`𝑹_1`$ ($`R_1<R_2`$) along this chosen line, where $`n_f(𝑹_1)=0`$. Such a point might or might not exist, and its presence is determined by the particular values of the chemical potentials $`\mu _{f,b}`$, and therefore by the particular total number of atoms and dimers. Naturally, if the number of dimers is sufficiently large, such a point $`𝑹_1`$, where $`n_f(𝑹_1)=0`$, exists. If such a point $`𝑹_1`$ exists, from that point on inward the atom number density vanishes identically and the dimer number density is determined from Eq. (9) alone. Moreover, the dimer number density continues to increase further on from the point $`𝑹_1`$ until it reaches the origin. The equipotential surface, which separates the bosons only region at the center of the trap, from the region with both fermions and bosons can be determined from the equation:

$$U_{fb}\mu _bU_{bb}\mu _f=(2U_{fb}U_{bb})V(𝑹_1).$$
(12)
One could have generated the same atom and dimer density profiles by starting from the origin instead. In that case one chooses a particular value for the dimer chemical potential $`\mu _b`$. One has however at that point to decide whether at the origin one has atoms as well or not. Depending on the presence of atoms at the origin, one proceeds in the reverse order and generates both atom and dimer number densities.

Once the number densities have been determined along one chosen line, one simply extends the values of the number densities along equipotential surfaces in the entire space. This overall onion structure is very reminiscent of the structures one has inside the Earth or neutron stars, with a core of one nature, followed by a mantle, which has a structure of a different nature, and eventually by a crust with a different structure as well. In the present case however all these three substructures are fluid. The extension of the present analysis to finite temperatures is straightforward, and one has to consider the free energy density, instead of the energy density, by adding the corresponding contributions arising from the fermionic and bosonic entropies as well as allowing for the bosons to have a non-condensed fraction. I shall limit the analysis here to the zero temperature case only. It is expected that with increasing temperature, first the fermionic superfluidity of the mantle is lost and subsequently the superfluidity of the bosonic core.

Some rather simple relations between these various radii can be obtained by taking into account the fact that the system is dilute, thus $`n_fa^31`$ and $`n_fba^31`$. One then easily obtains from Eqs. (10 \- 12) that the three radii discussed above in a harmonic trap are related (with a 2% error)

$$6R_2^2=R_3^2+5R_1^2.$$
(13)
Even though terms like $`U_{fb}n_f`$ have been dropped when this relation was obtained, the interaction is responsible for the specific numbers appearing here, as well as for the existence of this relation. In particular, this relations shows that the inner core disappears as soon as $`R_2<R_3/\sqrt{6}`$. The atom and dimer number density distribution are then given in a spherical trap by, see Fig. 1,

$`n_f`$ $`=`$ $`{\displaystyle \frac{1}{3\pi ^2}}\left({\displaystyle \frac{m\omega }{\mathrm{}}}\right)^3(0.2R_3^21.2R_2^2+r^2)^{3/2},`$ (15)
$`\mathrm{i}fR_1rR_2,`$
$`n_f`$ $`=`$ $`{\displaystyle \frac{1}{3\pi ^2}}\left({\displaystyle \frac{m\omega }{\mathrm{}}}\right)^3(R_3^2r^2)^{3/2},\mathrm{i}fR_2rR_3,`$ (16)
$`n_b`$ $`=`$ $`{\displaystyle \frac{1}{1.2\pi a}}\left({\displaystyle \frac{m\omega }{\mathrm{}}}\right)^2(R_2^2r^2),\mathrm{i}frR_2.`$ (17)
The magnitude, but not the shape, of the dimer distribution depends on the particular value of the scattering length. And since it has a very simple shape it is not shown in this figure. The specific numbers here are universal in the regime of large positive scattering lengths, when $`ar_0=(C_6m/\mathrm{}^2)^{1/4}`$. Here $`\omega `$ is the trap frequency and $`C_6`$ is the strength of the atom-atom van der Waals interaction. The atom distribution is thus essentially independent of the particular value of the scattering length $`a`$, as long as this is large and the system is dilute. One can easily derive analytical formulas for the number of atoms and dimers as functions of the radii $`R_{1,2,3}`$, even for anisotropic traps. Once the entire system is in the quantum degenerate regime the temperature dependence of the number density profiles are rather weak . The ratio of dimers to atoms scales as a function of $`\omega `$, $`a`$ and cloud radii as follows

$$\frac{N_{dimer}}{N_{atom}}\frac{\mathrm{}R_2^5}{m\omega aR_3^6}.$$
(18)
Since the scattering length $`a`$ is the smallest length scale typically the number of dimers is larger than the number of atoms, but this is not true always.

Besides the rather intricate structure of this new quantum object, one can naturally expect that the spectrum of collective states is going to have a nontrivial character, arguably different from the spectrum of BEC systems studied so far . As noted above, the core and the mantle of this new quantum object are superfluid, the first as a Bose system and the latter as mixed Bose-Fermi superfluids. As the most notable feature of a superfluid system is the possibility to excite a superflow, it would be extremely interesting to study the character of one or more vortex states. Unlike Bose  or Fermi  clouds, since the crust here is a normal Fermi fluid, a vortex will have to be limited to either the core and mantle or the mantle alone, and not extend into the Fermi normal crust (apart from the proximity effects due to the Andreev reflection mentioned above). It is likely that the Bose and Fermi superfluids can sustain largely independent superflows and the presence of a vortex in one component would not necessarily lead to the appearance of a vortex in the other component as well. The rotational properties of this object are clearly nontrivial, as one can expect both rotational and irrotational flows at the same time. The kinetics and thermalization of an object with such a complex structure is going to provide a interesting challenge to unravel, both experimentally and theoretically. However, unlike stars and planets, we should be able to study these processes under a variety of conditions and be able to finely tune various properties of these ”stars in a bottle.”

In summary, it was shown that an atomic cloud consisting of fermionic atoms in a trap, at sufficiently low temperatures and when the atom-atom scattering length is positive and large, will have an onion-like structure, with a superfluid bosonic core made of atomic dimers, a mantle consisting of a mixture of superfluid bosonic dimers and fermionic atoms, all surrounded by an outer layer consisting of fermions alone, forming a normal Fermi gas. The existence of this particular structure is due to a large extent to the fact that various couplings “conspired” to acquire the specific values quoted here, see Eqs. (3 - 5).

This work benefited from partial financial support under DOE contract DE–FG03–97ER41014.


text:
# A Hubble Space Telescope Census of Nuclear Star Clusters in Late-Type Spiral Galaxies: II. Cluster Sizes and Structural Parameter Correlations

## 1 Introduction

In most formation scenarios for spiral galaxies, the central bulge is the “trashbin of violent relaxation” where a dynamically hot stellar component has formed either through external potential perturbations such as early mergers of proto-galaxies (e.g. Carlberg, 1992), or perhaps via internal effects such as violent bar instabilities (Norman et al., 1996). While this “nature or nurture” debate is far from settled, the recent discovery of the so-called $`M`$-$`\sigma `$ relation (Gebhardt et al., 2000; Ferrarese & Merritt, 2000) has made it clear that the structure of the very nucleus, i.e. the black hole and its immediate vicinity ($`1\mathrm{pc}`$), is closely linked to the overall dynamical structure of the stellar bulge. The latest-type spirals, then, must have lived very sheltered and uneventful lives, since their central “trashbin” is virtually empty. Indeed, these disk-dominated, bulge-less galaxies (e.g. Matthews & Gallagher, 1997) often have slowly rising rotation curves that indicate a nearly homogeneous mass distribution on scales $`1\mathrm{kpc}`$. On these scales, gravity therefore hardly provides a vector pointing at the center, and it is not obvious that the nucleus of these galaxies is well-defined and a unique environment.

It therefore came as a surprise that recent observations with the Hubble Space Telescope (HST) have shown that the photocenter of many late-type spiral galaxies is occupied by a compact, luminous stellar cluster (Phillips et al., 1996; Carollo et al., 1998; Matthews et al., 1999). In a recent paper (Böker et al., 2002, hereafter Paper I), we have shown that at least 75% of spiral galaxies with Hubble types later than Sc harbor such a nuclear star cluster. The fraction of nuclear clusters (NCs) appears to be somewhat lower in earlier Hubble types, as indicated by the results of Carollo et al. (1998). However, smaller sample sizes and the observational difficulties due to the more complex nuclear morphologies of early-type spirals and the decreasing contrast between a central star cluster and a luminous bulges make this result somewhat uncertain. Typical luminosities of NCs are in the range $`10^6`$ \- $`10^7\mathrm{L}_{}`$ (Paper I) which makes them much brighter than average stellar clusters in the disks of nearby spiral galaxies (Larsen, 2002), and comparable to young “super star clusters” in luminous merger pairs (e.g. Whitmore et al., 1999) or circumnuclear starforming rings in spiral galaxies (e.g. Maoz et al., 2001).

Despite the recent progress, the formation mechanism of NCs remains largely a mystery. So far, no satisfying explanation has been put forward to explain the high gas densities that must have been present to enable the in-situ formation of such compact and luminous objects in the nuclei of these disk galaxies with shallow potentials. Plausible alternatives to in-situ formation include the infall of evolved clusters, possibly combined with more recent star formation events due to the accretion of gas in the galaxy nucleus. In order to gather evidence for or against any of these scenarios, it is essential to identify the stellar population(s) that comprise the NC, and thus to constrain the cluster formation history. Another crucial piece of information in this undertaking is an accurate estimate for the cluster mass.

The stellar masses of NCs can be derived - under the assumption of spherical symmetry - if one has accurate knowledge of (a) the stellar velocity dispersion and (b) the cluster light distribution. This kind of analysis has been successfully applied to NCs (Böker et al., 1999) as well as young clusters in the disks of nearby galaxies (Smith & Gallagher, 2001; Mengel et al., 2002). Our team has obtained high-resolution spectroscopy of $`15`$ clusters from the sample of Paper I. These data will address point (a) above, i.e. to accurately measure the stellar velocity dispersion of the clusters (Walcher et al., 2003, Walcher et al. 2004, in prep.). The purpose of the present paper is to discuss point (b), namely the size and spatial structure of the NC sample of Paper I. Because these objects are only barely resolved even with HST, one has to carefully correct for the effects of the instrumental point-spread function (PSF). Because standard deconvolution algorithms are not well suited in the realm of “almost point sources”, we have used the specialized software package ISHAPE (Larsen, 1999) to find the best match between the data and PSF-convolved analytical models for the cluster shapes. We describe the details of this exercise in § 2, and compare the resulting structural parameters of NCs to those of other cluster populations in § 3.

Because nuclear star clusters are a relatively new aspect of galaxy morphology, it seems a logical approach to search for correlations between the properties of NCs and those of their host galaxies. As has been shown by Carollo (1999) for spirals of earlier Hubble type, such empirical correlations might reveal clues about the mechanism(s) that regulate the gas supply in the very centers of spiral galaxies, and will certainly improve our understanding of nuclear activity in spirals. Moreover, because the presence of a NC appears to be a common feature in all types of spiral galaxies, the dynamical processes that lead to nuclear starbursts will be an important contribution to the discussion on bulge formation scenarios and a possible evolutionary connection between the various Hubble types. For these reasons, we present in § 4 a number of correlation plots between NC sizes/luminosities and various host galaxy properties such as size, Hubble type, HI mass, star formation rate, or total magnitude. We summarize our results in § 5.

## 2 Structural Analysis of Nuclear Star Clusters

### 2.1 Methodology

Star clusters in all but the closest galaxies are marginally resolved sources, even with the excellent resolution of HST images. In order to measure reliably their intrinsic sizes and structural parameters, it is therefore crucial to have both observations with high signal-to-noise ratio (SNR) and accurate knowledge of the instrumental point-spread function (PSF). Insufficient sampling of the PSF can complicate matters further. For example, the pixel size of the WFPC2 ($`0.1^{\prime \prime }`$ and $`0.0455^{\prime \prime }`$ for the WF and PC chips, respectively) in general does not provide Nyquist sampling of the HST resolution element $`R=\lambda /D`$. While sub-pixel dithering can alleviate this problem to some extent, HST observations of extragalactic star clusters in all but the closest galaxies are almost always confined to the regime of “marginally resolved”.

In this situation, none of the conventional deconvolution methods such as the Maximum Entropy Principle (Burch, Gull, & Skilling, 1983) or the Richardson-Lucy algorithm (Richardson, 1972; Lucy, 1974) work reliably. Instead, differential aperture photometry (e.g. the magnitude difference $`m_2m_4`$ between concentric apertures of 2 and 4 pixels, Holtzman et al., 1992) of stellar clusters and comparison to individual stars broadened by Gaussians has been used to estimate cluster sizes. However, a Gaussian is generally a poor approximation to the surface brightness profile of stellar clusters. In particular, the wings of a Gaussian fall off too rapidly, and cannot adequately describe the outer regions of either old globular clusters (GCs) in the Milky Way or young GCs in external galaxies.

More recently, a number of authors (e.g. Kundu & Whitmore, 1998; Larsen, 1999; Carlson & Holtzman, 2001) have independently developed techniques to determine more reliably the structure of stellar clusters in external galaxies. The various methods are conceptually similar in that they assume a parametric model for the intrinsic light distribution. A particular model is considered to be a good description of reality if - after convolution with the instrumental PSF - it matches the observations. The “best” model is the one that minimizes the $`\chi ^2`$-difference to the data. The most widely used analytic models for the projected density profile of stellar clusters are those of King (1962):

$$f(z)=\{\begin{array}{cc}(\frac{1}{\sqrt{1+z^2}}\frac{1}{\sqrt{1+c^2}})^2\hfill & \mathrm{for}z<c\hfill \\ 0\hfill & \mathrm{for}zc.\hfill \end{array}$$
(1)
Here, $`z`$ is defined by $`z^2=a_1x^2+a_2y^2+a_3xy`$, with $`x`$ and $`y`$ denoting the coordinates relative to the center of the profile and constants $`a_1`$, $`a_2`$, and $`a_3`$ which depend on the major axis, ellipticity and orientation of the model. The so-called concentration index $`c`$ is the ratio between the tidal radius $`r_\mathrm{t}`$ and the core radius $`r_\mathrm{c}`$.

In the regime of marginal resolution, the results of Carlson & Holtzman (2001) have shown that unless the SNR of the observations is exceptionally high ($`500`$), it is very difficult to distinguish between King models of different concentration indices. On the other hand, they find that the effective radius $`r_\mathrm{e}`$, i.e. the radius which contains half the cluster light in projection, can be fairly robustly estimated. The value of $`r_\mathrm{e}`$ does not depend strongly on the exact shape of the assumed model. This is in agreement with the results of Kundu & Whitmore (1998) and has been extensively tested by Larsen (1999). In what follows, we will describe our attempts to measure the effective radii of the NCs found in the survey of Paper I.

For our analysis, we have made use of Sören Larsen’s software package ISHAPE, which is described in detail in Larsen (1999). This software package has been well tested in a number of studies (e.g. Mengel et al., 2002; Larsen et al., 2002b), and has proven to be a robust and user-friendly tool for structural studies of compact star clusters. In brief, the program convolves a two-dimensional analytical model (e.g. Gauss-, Hubble-, Moffat-, or King-profiles) with a user-provided PSF, and finds the particular parameter set that best describes the observations, including the exact location of the cluster centroid on sub-pixel scales. Comparison with the pure PSF (i.e., a PSF-convolved $`\delta `$-function) allows to assess whether the clusters are resolved at HST resolution. Because $`r_\mathrm{e}`$ is not defined for functions with divergent integrals such as Hubble-type profiles, we have only attempted to fit the data with the above-mentioned King models as well as Moffat profiles which are defined as

$$f(z)=\frac{1}{(1+z^2)^i},$$
(2)
with $`z`$ defined as for the King model, and $`i>1`$ denoting the so-called power index.

In order to construct accurate PSFs for our analysis, we have used the TinyTim software (Krist & Hook, 2001), which accounts for the spectral passband used, the plate scale variations across the field of view, and the spectral energy distribution (SED) of the source. Although all of our images were taken through the same filter (F814W), we have constructed a separate PSF, according to the exact position of the NC on the detector array. We note, however, that this only marginally affects the results because the WFPC2 PSF varies only slightly across the PC field. The choice of the correct SED for constructing the PSF is somewhat more critical in estimating the NC sizes. Unfortunately, spectroscopic information is available for only a limited number of NCs. Early results of our spectroscopic follow-up survey (Walcher et al., 2003) show that the light of many NCs is dominated by a stellar population in the age range 50-$`300\mathrm{Myr}`$. This agrees with a number of case studies of individual clusters, e.g. in M 33 (Gordon et al., 1999), NGC 2403 (Davidge & Courteau, 2002), IC342 (Böker et al., 1999), or NGC 4449 (Böker et al., 2001). We therefore chose the SED of an A5 star for constructing the TinyTim PSF for all objects. We have verified that an alternate choice for the cluster SED, e.g. that of a cool giant with spectral type M3, does not significantly change the ISHAPE results.

### 2.2 Results

We have used ISHAPE to fit all 59 NCs of Paper I within the aperture radius that is dominated by emission from the NC (at least 5 $`r_\mathrm{e}`$). Some typical examples that illustrate the ISHAPE results are presented in Figure 1. We used only circularly symmetric models within ISHAPE because visual inspection of the fit results confirmed that the majority of clusters show little evidence for flattening. In a few cases, the clusters are asymmetric or have a complex morphology, most likely due to dust lanes or faint companion clusters. These objects were excluded from the ISHAPE analysis. In addition, we excluded those clusters with insufficient SNR to obtain reliable fits, i.e. with SNR $`<`$ 40 (Larsen, 1999).

The remaining 39 clusters were fit using King profiles with $`cr_\mathrm{t}/r_\mathrm{c}`$ values of 15, 30, and 100, as well as Moffat profiles with power indices of $`i=1.5`$ and $`i=2.5`$. In principle, the galaxy light underlying the NC must be subtracted before fitting the cluster profile. The ISHAPE code simply subtracts a constant background level, measured at the egde of the fitting aperture. While this algorithm certainly is less than perfect, and more sophisticated methods to fit the surface brightness profile of the galaxy disk could be employed, we have verified that the resulting values for $`r_\mathrm{e}`$ are not systematically biased by this shortcoming. In fact, the cluster magnitudes calculated by integrating the best-fit ISHAPE model agree well with those listed in Paper I which were derived with a more careful estimation of the galaxy background (see § 3.1). The median (mean) difference between these two sets of numbers is 0.07 (0.11) magnitudes when integrating the models over the same aperture used for the Paper I photometry.

The best-fitting models and the resulting effective radii, both in angular and physical dimensions, are summarized in Table 1. The majority (29/39) of NCs have $`r_\mathrm{e}`$ $`<5\mathrm{pc}`$, with median (mean) values of $`3.5\mathrm{pc}`$ ($`5.1\mathrm{pc}`$). As discussed further in § 3, these values are comparable to the typical size of old globular clusters in the Milky Way and other nearby galaxies. We emphasize that all clusters are resolved, i.e. models that simply rescale and position the TinyTim PSFs are inconsistent with the data by a large margin. For completeness, we also list in Table 1 the best-fitting analytic model. We caution, however, that in some cases, other models may also be consistent with the data. However, this has little impact on the uncertainties of the $`r_\mathrm{e}`$ measurements, as explained further below. Our goal is not to determine the detailed shape of the cluster profile, but merely to obtain a robust estimate for its size.

The uncertainties listed in Column 7 of Table 1 represent the range of $`r_\mathrm{e}`$ for the best-fitting model profile that yields acceptable fits to the data. Here, “acceptable” is defined by the $`3\sigma `$ confidence limit for one free parameter, i.e. $`\mathrm{\Delta }\chi ^2<9`$ (e.g. Press et al., 1992). Our analysis is conservative in the following sense: before evaluating the quality of the model fits, we rescaled all $`\chi ^2`$ values such that the best ISHAPE model for a given galaxy had a $`\chi ^2`$ value equal to the number of degrees of freedom, i.e. the number of data points (pixels in the aperture) minus the number of free parameters. This rescaling was necessary because ISHAPE appears to underestimate the intrinsic errors in the cluster images. As a result, many of the best $`\chi ^2`$ values do not formally represent “good” models according to $`\chi ^2`$-statistics. As the scale factor was always smaller than one, the $`\chi ^2`$ rescaling is equivalent to increasing the noise level of each data point or pixel, and consequently broadens the range of models that appear consistent with the data. Unfortunately, this rescaling of errors made a more rigorous statistical comparison between the different analytic profiles impossible.

To illustrate the analysis, Fig. 2 shows for the example of NGC 2805 the rescaled $`\chi ^2`$ values as a function of $`r_\mathrm{e}`$ for all models considered. In this case, a Moffat profile with a power index of 1.5 (MOFFAT15) and effective radius $`r_\mathrm{e}`$ $`=0.056^{\prime \prime }\pm 0.003^{\prime \prime }`$ yields the best fit to the data. The $`3\sigma `$ confidence limit and the corresponding range of $`r_\mathrm{e}`$ is indicated by the horizontal line in the inlay panel. As mentioned before, the effective radius is a rather robust quantity, i.e. different models yield similar values for $`r_\mathrm{e}`$, even if they are not formally consistent with the data according to the above criterion for $`\mathrm{\Delta }\chi ^2`$. In the case of NGC 2805, this can be seen in Fig. 2: the “best” values for $`r_\mathrm{e}`$ of all five models are within 20% of each other which is also the typical uncertainty listed in Table 1.

## 3 Discussion of Structural Parameters

### 3.1 Comparison to Milky Way Globular Clusters

A natural question to address with the data we have in hand is that of a possible size-luminosity relation for NCs. Since NCs in late-type disks are potential candidates for the seeds of bulge formation, they might follow the mass-radius correlation established for more massive ellipsoidal stellar populations, i.e. elliptical galaxies and bulges of intermediate- and early-type spirals (e.g. Kormendy, 1985; Bender et al., 1992). On the other hand, the old GCs in the Milky Way – which constitute the best-studied family of stellar clusters – are remarkably free of any such trend: the effective radii of Galactic GCs are independent of cluster luminosity or mass (van den Bergh et al., 1991; Djorgovski & Meylan, 1994; McLaughlin, 2000). The same is true for GCs in the nearby galaxies M 31 (Barmby et al., 2002), M 33 (Larsen et al., 2002a), and NGC 5128 (Harris et al., 2002).

In Fig. 3, we plot the absolute I-band magnitude of the NCs, $`M_I^C`$, against their effective radius, $`r_\mathrm{e}`$. The procedure for measuring $`M_I^C`$ is described in detail in Paper I. In brief, the challenge is to determine the light contribution of the stellar disk underlying the NCs. For this, we use two limiting models for the inward extrapolation of the galaxy disk profile. The $`M_I^C`$ values in Column 5 of 1 represent the average of the resulting two numbers; half their difference determines the size of the error bars in Fig. 3. Aside from a handful of bright clusters with effective radii larger than average, the bulk of the sample—28 of 39 clusters—is found in the narrow range $`1\mathrm{pc}r_\mathrm{e}5\mathrm{pc}`$, while the clusters range over a factor of $`100`$ in luminosity, from $`9`$ $`M_I`$ $`14`$. As we discuss below, relating luminosity to mass is nontrivial for NCs, but taken at face value, the situation illustrated in Fig. 3 is more reminiscent of GCs than of the bulges in early-type spirals.

It is not only the absence of a correlation between size and luminosity that appears to connect the NCs in our sample to GCs. As noted in § 2.2, the distribution of $`r_\mathrm{e}`$ for our NC sample is quantitatively similar to that of Galactic GCs. This is shown explicitly in Fig. 4, where we also compare the distributions of I-band magnitude for our NC sample and for 111 Milky Way GCs without an obvious or suspected core-collapse morphology, i.e., with well-determined $`r_\mathrm{e}`$ values.111The Milky Way GC data are taken from the February 2003 version of the online catalogue compiled by Harris (1996), which contains 141 clusters with known values for angular half-light radius and heliocentric distance, King concentration index, absolute V-band magnitude, and $`(BV)`$ color and excess. For 96 of these, Harris also gives an apparent $`(VI)`$ color, which we correct according to $`E(VI)=1.25E(BV)`$ (Cardelli et al., 1989) in order to derive $`M_I`$ from $`M_V`$. We simply assign the average $`(VI)_0=0.95`$ to the other 45 clusters. Evidently, although the median $`M_I`$ of NCs is fully 3–4 mags brighter than that of the Milky Way GCs, the widths of the two luminosity distributions (which corresponds to a distribution in mass for the GCs) are comparable. The median effective radius of the NCs is $`\overline{r_e}=3.5\mathrm{pc}`$ (interquartile range 2.4–5.0 pc) and that of the GCs is $`\overline{r_e}=3.3\mathrm{pc}`$ (interquartile range 2.4–5.2 pc). A two-sample Kolmogorov-Smirnov test confirms statistically the visual impression that the two $`r_\mathrm{e}`$ distributions are identical, with a probability $`P_{\mathrm{KS}}=0.75`$ that they have been drawn from the same parent distribution. We further note that the GCs in M31, M33, and NGC 5128 have essentially the same median $`r_\mathrm{e}`$ as the clusters in Fig. 4, and the same independence from cluster luminosity (see the references above).

Figure 5 provides an alternate illustration of this basic point, in a representation that allows for direct comparison with larger ellipsoidal stellar systems and comes closer to the “fundamental plane” typically plotted for such systems. Here, we show the average surface brightness, $`I_\mathrm{e}L_\mathrm{I}/(2\pi r_\mathrm{e}^2)`$, vs. effective radius $`r_\mathrm{e}`$ for the NCs, the 111 Galactic GCs from Harris’ (1996) catalogue, and a number of E galaxies, spiral bulges, and dwarf Ellipticals (dE’s) from the compilation of Burstein et al. (1997). The position of Milky Way GCs and NCs to the left of the plot shows again that the latter have essentially the same size as the former, but are substantially brighter and therefore at higher surface brightness.

Of additional interest is the disconnect between the GCs and NCs on the one hand and the elliptical galaxies and spiral bulges on the other, a dichotomy first noted by Kormendy (1985). As we mentioned at the beginning of this discussion, the bulges and galaxies show a size-luminosity relation, while GCs do not. Equivalently, GCs have an enormous spread in surface brightness, while bulges and ellipticals at a given $`r_\mathrm{e}`$ appear much more uniform. Although this fact is still not well understood, the alignment of our NCs with the GCs in Figure 5 would seem to underline a closer link between NCs and GCs, as opposed to one between NCs and classical galaxy bulges.

For completeness, we have included in Figure 5 data for the nuclei of 5 nucleated dE galaxies in Virgo, from the sample of Geha et al. (2002). As those authors have already demonstrated, with a different version of the “fundamental plane”, the dwarf-elliptical nuclei are also more similar to star clusters than to galaxy bulges.

### 3.2 Possible Interpretations

The apparent similarity between NCs and GCs could be an important clue to the nature and origin of the NCs themselves. First, however, their identification as an extension of the GC sequence to brighter magnitudes (Fig. 4) and higher surface brightness (Fig. 5) needs to be understood in terms of mass rather than luminosity alone. Of particular concern in this context is the interpretation of the spread of NC absolute magnitudes in Fig. 3 and Fig. 4. We consider two extreme scenarios—neither of which we mean necessarily to advocate as plausible per se, but which may serve to delimit further discussion:

(1) If the stellar populations of all NCs shared a common age, then the five-magnitude spread in $`M_I^C`$ would translate directly to a range of a factor of $``$100 in NC mass, intriguingly similar to the corresponding variation among old GCs (Fig. 4). Were Fig. 5 then to be re-plotted using mass surface densities, the NC locus would simply shift vertically, by an amount depending on the typical NC age (through the $`I`$-band mass-to-light ratio). This would strengthen the connection of NCs to GCs and sharpen their distinction from galaxy bulges.

(2) If instead all NCs had roughly the same stellar mass, then the spread in $`M_I^C`$ would have to arise from a range in luminosity-weighted age, with fainter cluster magnitudes (corresponding to systematically older ages) being purely the result of stronger population fading. If this were true, then Fig. 5 re-drawn in terms of mass surface densities would have the NC locus contracted into a relatively narrow line of slope $`2`$. The small variation in NC effective radius would then be largely the result of a small mass range in our sample, and it would reveal little about the (non)existence of a mass-radius relation among these clusters.

In order to gain some quantitative feel for the viability of these two alternatives, we have used the PEGASE population-synthesis code (Fioc & Rocca-Volmerange, 1997) to calculate the evolution of $`I`$-band mass-to-light ratio and absolute magnitude for a stellar cluster of mass $`7\times 10^5M_{}`$, formed in a single burst of star formation with the IMF of Kroupa et al. (1993) between lower- and upper-mass cut-offs of $`0.1M_{}`$ and $`120M_{}`$. The results, for a range of assumed metallicities, are presented in Fig. 6. Note that the mass-to-light ratio in the top panel is independent of the assumed total mass of the cluster, while the $`I`$-band magnitude in the bottom panel is of course sensitive to this parameter: for a given metallicity, the absolute magnitude of a single-burst stellar population of any mass $`m`$ is obtained, as a function of age, by adding $`2.5\mathrm{log}(m/7\times 10^5M_{})`$ to the appropriate curve in Fig. 6. The calculations for $`\mathrm{{\rm Y}}_I`$ and $`M_I^C`$ both take into account the mass lost from a cluster due to stellar evolution, i.e., winds and supernova explosions.

If the NCs were in fact a roughly coeval population of clusters with a narrow range in mass-to-light ratio $`\mathrm{{\rm Y}}_IM/L_I`$ (scenario (1) above), then they could in principle have masses similar to the Galactic GCs if $`\mathrm{{\rm Y}}_I`$ were some 15 times smaller for NCs than for GCs (to explain the 3-magnitude difference in the peaks of the $`M_I^C`$ distributions of Fig. 4). The top panel of Fig. 6 shows that a “typical” GC, with an age of 13 Gyr and $`[\mathrm{Fe}/\mathrm{H}]=1.5`$, has $`\mathrm{{\rm Y}}_I1.75M_{}L_{}^1`$ in these models.222The observed number is about $`1.2M_{}L_{}^1`$, which follows from an average $`V`$-band value of 1.45 in solar units (McLaughlin, 2000), and a typical $`VI`$ color of 0.95 (see § 3.1). The same graph then shows that the NCs would have to be of order $`10^7`$ years old (to within a factor of two or so) in order to have $`\mathrm{{\rm Y}}_I0.1M_{}L_{}^1`$, and thus a mass distribution similar to that of the GCs. Interestingly, this would make them potential analogues of the super star clusters forming in many nearby interacting and starburst galaxies (e.g., Whitmore et al., 1999). In this context, it is worth noting that these super star clusters tend to have effective radii $`r_e210`$ pc that do not appear to depend strongly on luminosity—reminiscent of NCs and GCs both (e.g., Barth et al., 1995; Zepf et al., 1999). Nevertheless, it is difficult to imagine how such a large fraction of late-type disks could have conspired simultaneously to form nuclear clusters such a short time ago.

If, on the other hand, scenario (2) were correct, and the spread in NC luminosity reflected an age spread among single-burst stellar populations with essentially a single mass, then the bottom panel of Fig. 6 suggests that this common mass would have to be near $`7\times 10^5M_{}`$, and the age range exceedingly broad: from $`10^7`$$`10^{10}`$ years to just barely reproduce the observed $`14M_I^C9`$ (Fig. 4). (A change from the “fiducial” mass in Fig. 6 by a factor of even 2.5 in either direction would shift all the $`M_I^C`$ curves up or down by a full magnitude, leaving an incomplete overlap with the observed NC magnitude distribution.)

Even as oversimplified as these two scenarios are, it is not possible to discriminate further between them without additional information. We have therefore obtained short-wavelength spectra for a number of the NCs in the current sample, which we will fit with population-synthesis models to constrain the clusters’ ages and star-formation histories, and from which we will extract stellar velocity dispersions to constrain dynamical mass-to-light ratios. Although this analysis is still in progress, the first results (Walcher et al., 2003) suggest a picture that is, not surprisingly, rather more complicated than either of the extremes just discussed.

Walcher et al. (2003) show that the light-dominating population of stars in 10 relatively bright NCs range in age from $`10^7`$$`10^9`$ years old. That is, there is in fact a spread in (luminosity-weighted) age among the NCs — although it does not appear sufficient to explain the full range of $`M_I^C`$ given the single-burst models in Fig. 6, unless there is also a significant (order-of-magnitude) spread in the cluster masses. More importantly, however, Walcher et al. (2003) show that the single-burst hypothesis itself is almost certainly incorrect: the spectra of their sample of clusters show some direct evidence for a mix of stellar ages, and the dynamical mass-to-light ratios they derive are consistently higher than the value of $`\mathrm{{\rm Y}}_I`$ inferred from population-synthesis modeling, suggesting that an underlying old (faint but massive) population of stars may be present.

In view of this, and given the similarity between GC and NC effective radii, an obvious first-order enhancement over the single-burst idea is the hypothesis that NC masses are in fact dominated by old stellar populations – perhaps bona fide GCs that have been dragged by dynamical friction to the centers of their host galaxies – but that their luminosities, and hence the young ages implied by their spectra, are due mainly to more recent star formation triggered by the infall of a smaller mass of circumnuclear gas. A more thorough exploration of this scenario is beyond the scope of this paper and must await the full analysis of our spectroscopic data (Walcher et al. 2004, in preparation). Nevertheless, it is worth noting that recent high-resolution interferometric studies of the nuclear gas dynamics in nearby galaxies have shown that molecular gas flows can indeed reach to within the central few pc of the nucleus (e.g. Schinnerer, Böker, & Meier, 2003) and thus provide a possible mechanism for the fueling of such nuclear starburst events.

## 4 Correlations with Host Galaxy Properties

### 4.1 Luminosity and Mass of the Host Galaxy

The work of Carollo et al. (1998) has shown that the luminosity of nuclear star clusters in early- and intermediate type spirals correlates with the total luminosity of the host galaxy. A similar result has been established by Lotz et al. (2001) for dE galaxies. In Figure 7a, we demonstrate that such a correlation is also found in spirals at the late end of the Hubble sequence. Here, we plot the apparent I-band magnitude of the NCs, $`m_I^C`$ (from Paper I), against the total apparent blue magnitude, $`m_B^G`$, of the host galaxy (Paturel, Bottinelli, & Gouguenheim, 1994) as listed in the Lyon-Meudon Extragalactic Database333http://leda.univ-lyon1.fr/ (LEDA). The relation still holds, albeit with a somewhat larger scatter, when plotting the absolute magnitudes $`M_I^C`$ and $`M_B^G`$ (Fig. 7b). The conversion from apparent to absolute magnitudes and the correction for Galactic reddening was performed using the foreground extinction and galaxy distances listed in Table 1 of Paper I. The solid lines in both panels indicate the best (least-squares) linear fit to the datapoints. The fit coefficients as well as their statistical significance (as measured by the Spearman rank-order coefficient) are listed in Table 2. The fact that despite the large intrinsic variations of the cluster luminosities discussed in §3.2, there still is a rather good correlation with galaxy luminosity suggests that the global properties of the host galaxy play a significant role for the formation processes of NCs.

One possible mechanism that could link the occurrence of nuclear starbursts to the total galaxy luminosity are tidal interactions or minor mergers, because these can both enhance the global star formation rate (SFR) and cause significant infall of matter towards the nucleus via re-distribution of angular momentum. However, as demonstrated in Fig. 8, NC luminosity appears only weakly (if at all) correlated with the IRAS far-infrared (FIR) magnitude $`M_{\mathrm{FIR}}`$444The apparent far-infrared magnitude, taken from LEDA, is calculated according to $`\mathrm{m}_{\mathrm{FIR}}=2.5\mathrm{log}(2.58\times \mathrm{f}_{60}+\mathrm{f}_{100})+14.75`$ where $`\mathrm{f}_{60}`$ and $`\mathrm{f}_{100}`$ are the IRAS fluxes at $`60\mu \mathrm{m}`$ and $`100\mu \mathrm{m}`$ (in Jansky). This expression is identical to the one used in the RC3. The distances listed in Paper I are then used to calculate $`M_{\mathrm{FIR}}`$. which is often used as a proxy for the SFR in galaxies (see e.g. Telesco, 1988, for a review). Assuming that tidal interactions contribute significantly to the SFR and hence the FIR luminosity in late-type spirals, this lack of a strong correlation between SFR and NC luminosity suggests that tidal interactions are not the driving force for the evolution of the very nucleus. This is supported by the WFPC2 images of our galaxy sample (Paper I) which show little evidence for recent star formation activity in their (central) disks. Instead, the smooth, shallow, and undisturbed stellar disks of many late-type spirals indicate a rather uneventful dynamical history, yet most of them have a luminous NC.

An alternative explanation for the strong correlation between $`M_I^C`$ and $`M_B^G`$ could be that the properties of the NC are more governed by galaxy mass than by the global star formation rate. In principle, HI linewidths could be used to address this possibility: the inclination-corrected rotational component of the HI linewidth

$$W_\mathrm{R}^i=(W_{20}W_{\mathrm{rand}})/\mathrm{sin}i,$$
(3)
based on the correction for random motions of Bottinelli et al. (1983), can be used as a proxy for the dynamical mass of the galaxy. Unfortunately, our galaxy sample was selected to have inclinations close to face-on and consequently the $`\mathrm{sin}i`$ corrections are highly uncertain. For the sake of completeness, we have netherless calculated the galaxy masses with the above expression. Following Matthews, van Driel & Gallagher (1998), we adopted $`W_{\mathrm{rand}}=20\mathrm{km}\mathrm{s}^1`$ for our galaxies, which corresponds to a line-of-sight velocity dispersion $`\sigma _z=5.5\mathrm{km}\mathrm{s}^1`$ (Rhee, 1996). However, using the inclination $`i`$ as listed in LEDA, we found that for our sample, the HI width does not even correlate with total galaxy luminosity which is a contradiction of the Tully-Fisher relation. We thus conclude that for our sample of face-on spirals, $`W_\mathrm{R}^i`$ is not a reliable indicator of galaxy mass, and that we therefore cannot draw any meaningful conclusions about possible correlations of NC properties and host galaxy mass.

### 4.2 Central Surface Brightness of the Galaxy Disk

Recent analysis of large galaxy samples (e.g. Blanton et al., 2002) has put on a solid statistical ground the notion that there are strong correlations between a number of physical properties of galaxies. In particular, Blanton et al. (2002) have shown that for disk-dominated, “exponential” galaxies (i.e., with Sérsic (1968) shape parameter $`n<1.5`$), surface brightness correlates with total luminosity. We confirm this finding for our sample, as demonstrated in Figure 9a. Here, we have plotted $`M_B^G`$, the total blue magnitude of the host galaxy, versus $`\mu _I^0`$, the central I-band surface brightness of the stellar disk underlying the NC. The values for $`\mu _I^0`$ are listed in Column 3 of Table 1; they were measured from our WFPC2 images by averaging the two inward extrapolations of the surface brightness profile shown in Fig. 3 of Paper I at a radius of $`0.024^{\prime \prime }`$ (i.e. the border of the central WFPC2 pixel).

Given that NC luminosity $`M_I^C`$ correlates with host galaxy luminosity (Fig. 7b), and that the latter correlates with $`\mu _I^0`$, it is not surprising that $`M_I^C`$ is also correlated with $`\mu _I^0`$, as demonstrated in Figure 9b. Again, the fit coefficients and significance levels of both correlations are listed in Table 2. While it is true in general that higher surface brightness disks make it more difficult to detect faint clusters, we emphasize that this effect does not result in a significant selection bias. In fact, most galaxies without clear identification of a NC have $`\mu _I^0`$ values below average (Figure 4 of Paper I). Note that the values of $`\mu _I^0`$ are as observed, i.e. they have not been corrected for inclination.

So far, we have shown that more luminous galaxies have higher surface brightness disks as well as more luminous NCs. Is this simply a reflection of the popular notion that “bigger galaxies have more of everything”? While this is certainly part of the answer, it is likely not the whole story. This is demonstrated in Fig. 10, which plots $`M_I^C`$, as well as the cluster-to-disk luminosity ratio $`C/D`$, as functions of the physical size of the galaxy disk, i.e., the major axis diameter of the $`\mu _\mathrm{B}=25`$ isophote as listed in LEDA. (The disk luminosity is calculated as total galaxy luminosity minus NC luminosity, but with $`C/D0.01`$ in general the disk light is not significantly different from the total galaxy luminosity.)

Apparently, $`M_I^C`$ depends weakly on galaxy size: the Spearman rank-order correlation coefficient in Fig. 10a is $`s=0.36`$, with slightly less than 3-$`\sigma `$ significance. However, this actually reflects the stronger and more significant correlation between NC luminosity and galaxy luminosity (§4.1), coupled with a similarly significant dependence of disk size on total disk luminosity: $`R_{\mathrm{disk}}L_{\mathrm{disk}}^{0.32}`$ in this sample. From Table 2 we have that NC luminosity scales with galaxy (disk) luminosity as $`L_{\mathrm{cl}}L_{\mathrm{disk}}^{0.78}`$, and thus we can expect $`C/D(L_{\mathrm{cl}}/L_{\mathrm{disk}})L_{\mathrm{disk}}^{0.22}R_{\mathrm{disk}}^{0.69}`$ even if there is no direct connection between NC luminosity and global disk size. This scaling is indeed consistent with the weak trend in Fig. 10b.

To first order, the NC luminosity is a measure of the efficiency with which gas is funneled towards the nucleus and converted into stars, regardless of the exact star formation history of the NC. This is true even if the cluster luminosity is dominated by a minor (in mass) but recent burst, because on average, younger populations imply a higher duty cycle of nuclear “delta”-bursts. It thus appears that, at least to some degree, the mechanisms that regulate the (past) star formation efficiency in the vicinity of the galaxy nucleus – as traced by the central surface brightness of its stellar disk – also play a significant role for the amount of gas that is funneled into the very nucleus and converted into stars. Further support for this notion comes from the fact that the molecular gas content within the central kpc of late-type spirals also appears to be connected to the disk surface brightness in the sense that low surface brightness disks in general have little or no central molecular gas (Fig. 7 of Böker, Lisenfeld & Schinnerer, 2003b).

### 4.3 Hubble-type

Figure 11a compares nuclear cluster absolute magnitudes against the Hubble type of their host galaxies. The very weak correlation apparent here, with a Spearman coefficient $`s=0.26`$, is only marginally significant (at the 2-$`\sigma `$ level, and less if the cluster in NGC 7418, with $`M_I^C=16.33`$, is excluded). It seems in any case to be the result of the strong link between NC and total galaxy luminosities, together with a trend towards slightly fainter galaxy magnitudes, on average, for the later Hubble types. Figure 11b therefore plots the ratio of cluster-to-disk luminosity vs. Hubble type. The correlation coefficient in this case is only $`s=0.13`$, which is not significant at even the 1-$`\sigma `$ level.

Thus, nuclear cluster luminosity is essentially independent of Hubble type. This should not come as a surprise given that the Hubble type is not a strong discriminant for the latest-type spirals, as discussed in more detail by Böker, Stanek, & van der Marel (2003). The main reason is that the NC is quite easily mistaken for a compact stellar bulge in seeing-limited images, and hence a ground-based morphological classification based at least in part on the bulge-to-disk ratio is inaccurate at best.

### 4.4 Nuclear Clusters and Stellar Bars

Using the morphological classifications of our galaxy sample, both from the RC3 (listed in Table 1 of Paper I) and the LEDA database, we have investigated for our complete Paper I sample whether nuclear star clusters are more likely to be found in (or whether they preferentially avoid) galaxies of a certain bar class. Table 3 summarizes the statistics for both subsamples with (59 galaxies) and without (18 galaxies) a NC. The fraction of galaxies with a NC is 76%, 88%, and 67% for the RC3 classifications A, AB, and B, respectively. The LEDA database only lists the types “unbarred” and “barred”; for these, the corresponding fractions are 88% and 73%. The apparently smaller frequency of NCs in barred galaxies is probably not significant, given the small number statistics in some categories, and the rather high number of uncertain classifications (17 galaxies in our sample are classified either as uncertain or peculiar in the RC3). More reliable classifications of the morphology of late-type spirals will require a more accurate isophotal analysis or even kinematic information on their central disks.

This rough analysis suggests that the galaxy bar class is not a significant factor for the likelihood of hosting a NC. In particular, the late-type galaxies in our sample seem to have no especial difficulty in supporting bars and nuclear clusters simultaneously. This might have seemed somewhat paradoxical in the light of some early numerical simulations of barred galaxies, which found that fairly modest central mass concentrations could drive dynamical instabilities that destroyed the bar on rather short timescales (e.g., Friedli, 1994; Norman, Sellwood, & Hasan, 1996). However, more recent work (Shen & Sellwood, 2003) suggests that even an NC as massive as a few percent of the total disk mass of a galaxy (i.e., with $`C/D0.02`$, already much higher than the typical values in Fig. 10b) may not be able to destroy a bar completely within a Hubble time. It is then somewhat unclear what role bars and NCs together might play in the formation of the exponential “pseudo”-bulges found in many late-type spirals (cf. Carollo et al., 2001).

## 5 Summary

We have presented a structural analysis of 39 nuclear star clusters in late-type spiral galaxies. After correction for the instrumental point spread function, we measure a median effective radius for the sample of $`\overline{r_\mathrm{e}}=3.5`$ pc, with 50% of the sample falling between $`2.4\mathrm{pc}r_\mathrm{e}5.0\mathrm{pc}`$. This narrow size distribution is statistically indistinguishable from that of Galactic globular clusters, even though the nuclear clusters are on average 4 magnitudes brighter than the old globulars. The compactness of NCs clearly separates them from larger ellipsoidal stellar systems such as spiral bulges or elliptical galaxies.

In addition, we have compared the luminosity of NCs with various properties of their host galaxies. We find that NC luminosity is correlated with galaxy luminosity which extends a previous result obtained for earlier-type spirals. Because there is a strong correlation between total luminosity and central disk surface brightness for the galaxies in our sample, NC luminosity also correlates with local disk brightness. It is not clear which of these trends is physically more fundamental. Either way, other apparent correlations such as between cluster luminosity and galaxy size can ultimately be explained by the underlying link to galaxy luminosity and/or central disk surface brightness.

We find no systematic trends with the Hubble type of the galaxy, confirming that the latter is not a strong discriminant for the latest-type spirals. We also find no evidence for a connection between the presence/absence of NCs and the presence/absence of large-scale stellar bars in the host galaxy.

Without additional spectral information, it is not possible to put strong constraints on the formation history of NCs. We suggest that the most likely scenario for NC formation consists of multiple nuclear starbursts, possibly caused by gas infall onto a “seed” cluster. However, a more detailed discussion of this hypothesis has to await the analysis of our spectroscopic data which we will present in a future paper.

We are grateful to S. Larsen for help with the implementation and use of his ISHAPE software package. Support for this work was provided by the National Aeronautics and Space Administration (NASA) through grants for project no. GO-8599 awarded by the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract no. NAS5-26555. 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 NASA. It has also benefited greatly from use of the Lyon-Meudon Extragalactic Database (LEDA, http://leda.univ-lyon1.fr).


text:
# Self-sustained activity in a small-world network of excitable neurons

## Abstract

We study the dynamics of excitable integrate-and-fire neurons in a small-world network. At low densities $`p`$ of directed random connections, a localized transient stimulus results in either self-sustained persistent activity or in a brief transient followed by failure. Averages over the quenched ensemble reveal that the probability of failure changes from $`0`$ to $`1`$ over a narrow range in $`p`$; this failure transition can be described analytically through an extension of an existing mean-field result. Exceedingly long transients emerge at higher densities $`p`$; their activity patterns are disordered, in contrast to the mostly periodic persistent patterns observed at low $`p`$. The times at which such patterns die out are consistent with a stretched-exponential distribution, which depends sensitively on the propagation velocity of the excitation.

Recent research in complex networks has provided increasing evidence for their relevance to a variety of physical, biological, and social phenomena Watts and Strogatz (1998); Strogatz (2001); Albert and Barabasi (2002). Two distinct types of topology have been particularly useful in providing insights into the implications of complex connectivity: scale-free networks Albert and Barabasi (2002), characterized by the existence of a small number of hubs with high coordination number, and small-world networks Watts and Strogatz (1998), characterized by the presence of shortcuts that link two randomly chosen sites regardless of the distance between them.

So far, most work on complex networks has focused on their topological and geometrical properties; less attention has been given to the properties of dynamical systems defined on such networks. The interplay between the intrinsic dynamics of the constituent elements and their complex pattern of connectivity strongly affects the collective dynamics of the resulting system. For instance, the addition of shortcuts induces a finite-temperature phase transition even in the one-dimensional Ising model Pekalski (2001) and the introduction of unidirectional shortcuts can change the second-order phase transition in the two-dimensional Ising model into a first-order one Sánchez et al. (2002). In a system of coupled oscillatory elements, the introduction of shortcuts enhances synchronization Barahona and Pecora (2002), while the introduction of hubs eliminates the threshold for epidemic propagation Pastor-Satorras and Vespignani (2001).

The coexistence of shortcuts and regular local connections characteristic of small-world networks (SWNs) mimics a salient feature of the circuitry in the cortex Mountcastle (1997); González-Burgos et al. (2000); McCormick and Contreras (2001); Wu et al. (1999); Demir et al. (2001); Melchitzky et al. (2001), where experimental observations of excitatory traveling waves Wu et al. (1999) provide evidence of some degree of local connectivity, while it is also recognized that long-range excitatory connections are present Demir et al. (2001); McCormick and Contreras (2001). Our goal is to explore the influence of this complex connectivity on the dynamics of neuronal circuits; to this purpose, we choose a minimal model. The underlying network is modeled as a SWN with unidirectional shortcuts that reflect the nonreciprocal character of synaptic connections, and the excitable neurons are modeled as leaky integrate-and-fire units. We find that even this simple model exhibits a rich repertoire of distinct dynamical behaviors as a function of the density $`p`$ of added shortcuts: a low $`p`$ regime characterized by persistent periodic activity that is bistable with the quiescent state, a transition to failure with increasing $`p`$, followed by a reemergence of long-lasting disordered activity. We note that a SWN with unidirectional shortcuts has been considered in a different regime by Lago-Fernández et al Lago-Fernández et al. (2000, 2001) to address the possibility of rapid synchronization among conductance-based neurons of the Hodgkin-Huxley type.

The model considered here consists of a one-dimensional array of $`N`$ integrate-and-fire neurons (IFNs) in which a SWN topology is created through the addition of a density $`p`$ of unidirectional long-range couplings. The membrane potential of the IFNs is determined by

$`\tau _m{\displaystyle \frac{dV_i}{dt}}=V_i+I_{ext}+g_{syn}{\displaystyle \underset{j,m}{}}w_{ij}\delta (tt_j^{(m)}\tau _D).`$ (1)
A neuron fires whenever its voltage exceeds $`1`$; the voltage is then reset to $`0`$. The external current is chosen to satisfy $`I_{ext}<1`$; in this regime the IFNs are not oscillatory, but excitable. The last term models input currents due to presynaptic firing as a delayed impulse: if $`w_{ij}`$ is non-zero, then neuron $`i`$ receives a pulse input of amplitude $`g_{syn}`$ with a delay $`\tau _D`$ after neuron $`j`$ has fired its $`m`$th spike at time $`t_j^{(m)}`$. The synaptic conductance is chosen to satisfy $`I_{ext}+g_{syn}>1`$, so that a single input suffices to sustain firing activity. The local connections are modeled here as nearest-neighbor couplings ($`w_{i,i\pm 1}=1`$) that define an underlying regular lattice. The long-range connections result from randomly adding rather than rerouting Watts and Strogatz (1998) a fixed fraction $`pN`$ of unidirectional couplings $`w_{ij}=1`$ to generate a SWN topology.

At $`p=0`$, any excitation sufficient to cause a neuron to fire will generate two pulses that propagate through the regular lattice in opposite directions with velocity $`v=1/\tau _D`$, and either exit the system or annihilate each other, depending on boundary conditions. No persistent activity results in either case. However, self-sustained activity may arise for nonzero $`p`$, as shown in fig.1. Persistence relies on the reinjection of activity via a shortcut into a previously active domain that has by then recovered; this reinjection can occur only if the shortcuts are unidirectional. For a fixed value of $`p`$, any particular network realization has a different connectivity graph that may or may not sustain persistent activity. We average over typically 2000 realizations to calculate the probability of persistent activity; the complementary probability of failure to sustain activity is shown in fig.2 (inset) as a function of the density $`p`$ of random connections for four different system sizes. In this regime, the probability of failure is an increasing function of $`p`$ that crosses over from $`0`$ to $`1`$ with increasing steepness as the size $`N`$ of the system increases.

Failure to sustain activity is a simple consequence of the intrinsic dynamics of the neurons. Pulses travel outwards from an initial activity seed, and spawn new pulses at a rate that increases with $`p`$. A currently inactive neuron can fire again only if enough time has elapsed from its preceding firing to allow for a recovery to $`V1g_{syn}`$. A single input will be able to elicit a spike only if the elapsed time exceeds $`T_R^{(1)}`$, with

$$T_R^{(n)}=\tau _m\mathrm{ln}\left(\frac{I_{ext}}{I_{ext}+ng_{syn}1}\right).$$
(2)
If activity recurs to a given site too rapidly, the neuron will fail to produce a spike, and the pulse of activity will die out. A critical density $`p_{cr}`$ for the transition from self-sustained activity to failure can be estimated from

$$T_A(p_{cr})=T_R^{(1)},$$
(3)
where $`T_A(p)`$ is the time needed for the activity to spread across the whole network. At a fixed velocity for pulse propagation, this time corresponds to the largest distance across the network. This distance has been calculated for bidirectional shortcuts using a mean-field approach Newman et al. (2000); when extended to the case of unidirectional shortcuts it results in

$$\sqrt{\left(1+\frac{4}{pN}\right)}\mathrm{tanh}\left[\sqrt{\left(1+\frac{4}{pN}\right)}\frac{pT_A(p)}{2\tau _D}\right]=1.$$
(4)
An implicit expression for $`p_{cr}`$ and its dependence on the system size $`N`$ and the propagation velocity $`v`$ follows from combining equations (3) and (4); for large $`N`$, $`p_{cr}(N)\mathrm{ln}N`$. Failure rate curves as a function of $`(pp_{cr}(N))/p_{cr}(N)`$, shown in fig.2, cross at the theoretically predicted value. This observation, together with the increased steepness of these curves with increasing $`N`$, indicates that a well-defined transition to failure occurs in the thermodynamic limit.

This well-defined transition to failure occurs only for sufficiently fast waves, i.e. for short delay $`\tau _D`$. For larger $`\tau _D`$ the dynamics of the system become quite more complex, and the fraction of realizations that fail before a fixed time ($`t_{max}=2000`$ in fig.3) becomes a nonmonotonic function of $`p`$. While at low $`p`$ the firing patterns are highly regular (cf. fig.1) and all failures occur within one or two cycles of the initial activity, for higher $`p`$ the patterns are more disordered (cf. fig.4) and the activity can persist for a very long time before failure. Consequently, the distribution of failure times, shown in fig.5, exhibits an increasingly long tail for longer delay times.

To understand the persistence of activity beyond $`p_{cr}`$, it is important to recognize that the result for the critical density $`p_{cr}(N)`$ hinges on the assumption that each neuron receives but a single excitatory input during each cycle of network activity. Its recovery time is therefore given by $`T_R^{(1)}`$, which sets a lower bound for the interspike interval (ISI). While this assumption is well satisfied for small $`p`$, it does not hold for $`p=O(1)`$. In fact, the likelihood that a neuron has $`n`$ incoming shortcuts follows a multinomial distribution such that the fraction of neurons receiving two incoming connections grows from about 0.05% at $`p=0.1`$ to about 30% at $`p=1`$. Neurons that receive $`n`$ inputs can have a recovery time as short as $`T_R^{(n)}`$. Fig.4 reveals that such neurons, with ISIs lower than $`T_R^{(1)}`$, play a crucial role in maintaining network activity in the regime $`p1`$. While neurons with ISIs greater than $`T_R^{(1)}`$ (marked black in fig.4) can go through silent epochs with near-zero activity, neurons with shorter ISIs (marked red in fig.4) may fire two or three times during a network cycle and carry over the network activity across these silent epochs.

The persistence of activity beyond $`p_{cr}`$ depends sensitively on $`\tau _D`$; this reflects the fact that persistence at higher densities $`p`$ is due to chains or trees of neurons that bridge the silent epochs due to their multiple inputs. Since activity propagates with a fixed speed, a chain of multiple-input neurons of given length can bridge a time interval proportional to $`\tau _D`$. Thus, as $`\tau _D`$ increases, ever shorter chains can contribute to bridging a silent epoch of given duration; the likelihood for failure will decrease accordingly. This picture is of course overly simplistic: whether a topological chain can be utilized as a dynamical bridge over a given silent epoch depends on the amount and timing of the inputs it receives, which in turn depend on the recent history of the entire network. Simulations reveal that the identity of the neurons that form the ‘bridging’ dynamical chains varies from cycle to cycle in an irregular way. This implies that the effective utilization of a dynamical chain over one cycle does not guarantee its availability on the next cycle. Therefore, even systems that persist for long times may still have a finite probability of failing.

The cumulative failure distribution function $`F(t)`$, shown in fig.5 for various values of $`\tau _D`$, exhibits a long tail and is well fit with stretched exponentials: $`F(t)=f_{\mathrm{}}(\tau _D)Ce^{\alpha t^\beta }`$ with $`\beta 0.4`$ (dotted lines). Even though the fits are based on runs up to $`t=300,000`$ for $`\tau _D=1.65`$ and $`\tau _D=1.7`$, they do not provide a value of $`f_{\mathrm{}}(\tau _D)`$ accurate enough for establishing whether true persistent activity exists for a small fraction of the network realizations ($`0.97f_{\mathrm{}}(\tau _D=1.65)1`$). Strikingly, the dependence of $`F(t)`$ on the delay time exhibits a high degree of structure, suggestive of ‘resonances’ at values of $`\tau _D`$ for which some chains and trees can be optimally utilized.

One of the salient features of the emergent dynamics of the model is persistent self-sustained activity at low densities $`p`$ of shortcuts. In this regime, the network is in fact bistable between ‘off’ and ‘on’ states, and can be switched between them with sufficiently large stimuli, as illustrated in fig.6. The synchronous stimulation of a sufficiently large number of neurons while the network is in the ‘on’ state increases the level of activity and effectively pushes the network to the right of the failure transition (cf. fig.2), causing a transition into the ‘off’ state.

To address the fact that real neuronal systems are noisy, the simulation shown in fig.6 includes Gaussian fluctuations in the membrane potential; their amplitude is chosen so as to cause the neurons to spike irregularly at a low rate. To keep the noisy spiking of a single neuron from generating a traveling pulse and initiating the ‘on’ state, we adjust the synaptic conductance so that several adjacent neurons must fire in rapid succession in order to propagate a pulse of excitation. The network topology is modified accordingly: we extend the local coupling to include up to $`2k`$ neighbors ($`w_{i,i\pm j}=1`$ for $`j=1\mathrm{}k`$) and model long-range connections via a population of intermediate excitatory neurons which both receive input from and project to multiple adjacent neurons. Under these conditions, spontaneous activity is highly unlikely to initiate traveling pulses. However, a sufficiently large stimulus, synchronous across several neurons, can again turn the state of elevated activity on and off (cf. fig.6); bistability is thus robust with respect to noise.

Network bistability has been hypothesized to be the neural correlate underlying the type of short-term memory known as working memory in the prefrontal cortex of monkeys and humans. Much more realistic and physiologically plausible models of cortical layers have been studied within the context of working memory (e.g. Compte et al. (2000)). Yet, not much attention has been given to heterogeneities in network topology or to long-range excitatory connections. The work presented here suggests that closer attention be given to the role of connectivity as an additional factor that contributes to the generation of the persistent, active state associated with working memory.

In conclusion, we have investigated the effect of incorporating random unidirectional shortcuts to a one-dimensional network of locally coupled integrate-and-fire neurons. We find that even a very low density of shortcuts suffices to generate persistent activity from a local stimulus through the reinjection of activity into previously excited domains. As the density of shortcuts is increased, the substantial decrease in the effective system size characteristic of small-world networks causes a crossover into a regime characterized by failure to sustain activity for essentially all network configurations. For sufficiently slow propagation velocities of the activity and sufficiently high shortcut densities, an intriguing second crossover occurs into a regime in which the activity still fails but only after often exceedingly long and strongly chaotic transients.

The complex dynamical phenomena we find in this extremely simple model are based on a robust mechanism: propagating pulses of activity that are sustained by branching and reinjection. We therefore expect that more realistic models of neuronal networks, which may include multiple ion channels and continuous synaptic currents as well as inhibitory coupling, will show qualitatively the same behavior upon the addition of shortcuts if they originally support propagating pulses of activity Golomb and Ermentrout (2001) that annihilate upon collision. Preliminary computations show that the failure transition persists if shortcuts are not allowed to exceed a maximal length $`L_m`$, as long as $`L_m\tau _D>T_R^{(1)}`$ Roxin et al. (unpublished). The phenomena should also be accessible in excitable chemical systems Martinez et al. (2002), in which the shortcuts could be implemented using video feedback.

We gratefully acknowledge support from the Department of Energy (DE-FG02-92ER14303), NSF grant DMS-9804673, and the NSF-IGERT program Dynamics of Complex Systems in Science and Engineering (DGE-9987577).