text:
# Fractal clustering of inertial particles in random flows

## Abstract

It is shown that preferential concentrations of inertial (finite-size) particle suspensions in turbulent flows follow from the dissipative nature of their dynamics. In phase space, particle trajectories converge toward a dynamical fractal attractor. Below a critical Stokes number (non-dimensional viscous friction time), the projection on position space is a dynamical fractal cluster; above this number, particles are space filling. Numerical simulations and semi-heuristic theory illustrating such effects are presented for a simple model of inertial particle dynamics.

Dust particles, droplets, bubbles and various impurities advected by turbulent flow have usually a finite size and a mass density differing from that of the carrier fluid. Contrary to passive tracers, whose dynamics is conservative (when the flow is incompressible), the dynamics of such *inertial particles* is rendered dissipative by the Stokes drag, which can lead to strongly inhomogeneous spatial distributions. The full statistical description of such preferential concentrations is still an open question with many natural and industrial applications, such as the growth of rain drops in subtropical clouds,ffs02  the formation of planetesimals in the early Solar system,bcps99  optimization of combustion processes and the coexistence problems between several species of plankton.kpstt00

Maxey and Rileymr83  derived an equation for the motion of a rigid spherical particle embedded in an incompressible flow. They assume that (i) the particle is smaller than the smallest turbulent scale of the carrier flow and that (ii) the Reynolds number associated to its size and its relative velocity with respect to the fluid is sufficiently small to approximate the surrounding flow by a Stokes flow. Then, the forces exerted on the particle are buoyancy, the force due to the undisturbed flow, the Stokes viscous drag, the added mass effect and the Basset–Boussinesq history force. Because of the complexity of the resulting equation of motion, simpler models are generally used. For instance, when the Stokes drag is very strong, the dynamics is close to that of passive tracer particles and the discrepancy can be captured by a spatial Taylor expansion, leading to a model in which the particles are advected by a synthetic flow comprising a small compressible component.ekr96 ; bff01  What singles out the model proposed here is its ability to take into account the full phase-space dynamics of the particles and to capture the essential features of their dissipative motion. We are interested in the “Batchelor régime” of the particles, meaning that we focus on spatial scales smaller than the Kolmogorov dissipation scale $`\eta `$. After rescaling of space, time and velocity respectively by factors $`\eta `$, $`\eta ^2/\nu `$ and $`\nu /\eta `$, and assuming that the particle velocity is sufficiently close to that of the fluid, the Newton equation satisfied by its trajectory $`𝐗(t)`$ reduces to

$$\ddot{𝐗}=\beta \frac{d}{dt}\left[𝒖(𝐗,t)\right]\frac{1}{\mathrm{S}_\eta }\left[\dot{𝐗}𝒖(𝐗,t)\right],$$
(1)
where $`\beta 3\rho _f/(\rho _f+2\rho _p)`$ is the added-mass factor ($`\rho _f`$ and $`\rho _p`$ are the fluid and the particle mass densities) and $`\mathrm{S}_\eta a^2/(3\beta \eta ^2)`$ is the *Stokes number* associated to the dissipative scale of the carrier flow ($`a`$ being the particle radius). Introducing the *co-velocity* $`𝐕\dot{𝐗}\beta 𝒖(𝐗,t)`$, the equation of motion can be interpreted in terms of the $`(2\times d)`$-dimensional dynamical system

$`\dot{𝐗}`$ $`=`$ $`\beta 𝒖(𝐗,t)+𝐕,`$ (2)
$`\dot{𝐕}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{S}_\eta }}\left[(1\beta )𝒖(𝐗,t)𝐕\right].`$ (3)
The motion of the particles is clearly dissipative, even if the carrier flow is itself incompressible : indeed, when $`𝒖=0`$, the contraction rate in phase space reduces to $`d/\mathrm{S}_\eta `$ and is strictly negative, inducing a uniform contraction. As a consequence, the long-time dynamics of the particles is characterized by the presence of an attractor, that is a dynamical fractal set of the phase space toward which the trajectories $`(𝐗(t),𝐕(t))`$ converge. Important information on the dynamical system (2)-(3), regarding stability, Lyapunov exponents, etc is obtained from the linearized equation governing the separation $`𝐑(t)(\delta 𝐗(t),\delta 𝐕(t))`$ between two infinitesimally close trajectories of the phase space. For scales within the viscous scale of turbulence, the velocity field $`𝒖`$ can be considered spatially smooth and the separation $`𝐑(t)`$ obeys the linear differential equation

$$\dot{𝐑}=_t𝐑,_t\left[\begin{array}{cc}\beta 𝝈(t)\hfill & \hfill _d\\ \frac{1\beta }{\mathrm{S}_\eta }𝝈(t)\hfill & \hfill \frac{1}{\mathrm{S}_\eta }_d\end{array}\right],$$
(4)
where $`𝝈`$ is the strain matrix of the carrier flow along the path of a reference particle : $`\sigma _{ij}(t)_ju_i(𝑿(t),t)`$, and $`_d`$ is the $`d`$-dimensional identity matrix. A full stability analysis of the dynamics can easily be done;b02  it relates the eigenvalues of the evolution matrix $`_t`$ to that of the stress tensor $`𝝈`$ representative of the local structure of the carrier flow. In both two and three dimensions, this leads to distinguishing between heavy ($`\beta <1`$) and light particles ($`\beta >1`$): the former are usually ejected from the elliptic regions, while the latter may cluster there in a pointwise manner. We therefore expect the vortex cores to be regions of high concentration of light particles and of low concentration of heavy particles, feature which is generally observed in experiments and simulations (for a review, see Eaton and Fessleref94 ).

We focus here on suspensions of particles with a volume fraction sufficiently small to neglect their interactions and collisions. Typically, the phase-space attractor on which the particles concentrate is a fractal object which may be characterized by various dimensions, in particular a non-random Hausdorff dimension $`d_H`$. As the position of the particles is obtained by projection from the $`2d`$-dimensional phase space onto the $`d`$-dimensional physical space, the convergence to the attractor is responsible for strong inhomogeneities in the large-times distribution of particles. More precisely, a standard result of the geometrical theory of fractal setsf86  states that if $`d_H<d`$, the distribution of particles in the physical space is itself a fractal set with Hausdorff dimension $`d_H`$, whereas if $`d_H>d`$, the particles fill the whole space. Hence, depending on the value of the dimension $`d_H`$ of the attractor, two different régimes are distinguished. Clearly, the dimension of the attractor is a function of the Stokes number $`\mathrm{S}_\eta `$ and of the added mass parameter $`\beta `$, and generally also depends on the statistical properties of the velocity of the carrier fluid. Leaving aside this latter dependence, let us first note what can be easily inferred on the behavior of $`d_H`$ as a function of $`\mathrm{S}_\eta `$ and $`\beta `$. First, in the limit of vanishing Stokes numbers, there is a reduction of dimensionality and the dynamics of simple tracers is recovered. An initially uniform distribution of particles remains uniform and we hence have $`d_Hd`$. Next, for very large Stokes numbers, the particles are less and less influenced by the carrier fluid and their motion becomes ballistic. They thus fill the whole phase space and we have $`d_H2d`$. In between these two asymptotic régimes, although it is not obvious that the dimension $`d_H`$ can fall below $`d`$, we shall actually show that there is a whole range of Stokes numbers for which $`d_H<d`$ and, thus, preferential concentration on fractal clusters occurs in physical space.

Finding theoretically or numerically the Hausdorff dimension $`d_H`$ of the attractor is not, in general, a simple task: its determination demands a full understanding of the global dynamics and its numerical measurement requires very large numbers of particles. To obtain a simple estimate of the attractor dimension, Kaplan and Yorkeky79  proposed to use the *Lyapunov dimension*. It is given by the Lyapunov exponents $`\lambda _1>\mathrm{}>\lambda _{2d}`$ which measure the exponential growth of infinitesimal distances, surfaces and volumes in the phase space and are expressed in terms of the limiting singular values of the Jacobi matrix, i.e.

$$\lambda _j\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}\left|𝒥_t𝒆_j\right|,𝒥_t𝒯\mathrm{exp}_0^t_s𝑑s,$$
(5)
where the $`𝒆_j`$’s are the eigenvectors of the symmetric matrix $`𝒥_t^\mathrm{T}𝒥_t`$ and $`𝒯\mathrm{exp}`$ denotes the time-ordered exponential. The Lyapunov dimension is defined as

$$d_Lj\frac{\lambda _1+\mathrm{}+\lambda _j}{\lambda _{j+1}},$$
(6)
where $`j`$ satisfies $`\lambda _1+\mathrm{}+\lambda _j0`$ and $`\lambda _1+\mathrm{}+\lambda _{j+1}<0`$. Beside being a simple estimate of the dimension of the attractor, it was actually showndo80  that the Lyapunov dimension gives a rigorous upper bound for the Hausdorff dimension $`d_H`$ of the attractor.

We performed numerical experiments in two and three dimensions for a space-periodic carrier flow generated randomly by the superposition of few independent Gaussian Fourier modes with a correlation time of the order of unity (this specific form for the carrier flow was considered by Sigurgeirsson and Stuartss02a  who proved the existence of a random dynamical attractor). The Lyapunov exponents are calculated by the use of the standard technique of Benettin et al.bggs80 , and the resulting Lyapunov dimension is represented both for $`d=2`$ and $`d=3`$ in Fig. 1 for the case of very heavy particles ($`\beta =0`$). Two important observations can be made from these simulations. First, fractal clustering of particles already occurs at very small Stokes number where the Lyapunov dimension behaves as $`d_LdC\mathrm{S}_\eta ^2`$ with $`C>0`$. This quadratic behavior near zero was predicted for the correlation dimension by Balkovsky et al.bff01  using the method of the synthetic compressible flow cited earlier, as an approximation at low Stokes numbers. The second observation is the presence of a critical value for the Stokes number below which the attractor dimension is smaller than $`d`$ and where the particles form fractal clusters in the physical space. The two régimes corresponding to different values of the Stokes number are illustrated for $`d=2`$ in Fig. 2. When the Stokes number is below the critical value (a), the particles concentrate onto a fractal set and both very dense and almost empty regions appear. On the contrary, for a Stokes number above the threshold (b), the particles fill the whole domain, albeit with a non uniform density.

When the particles have a finite mass ($`\beta 0`$), there also exists a critical Stokes number for their concentration onto fractal clusters. Figure 3 shows for $`d=2`$ the phase diagram obtained numerically which divides the parameter space $`(\beta ,\mathrm{S}_\eta )`$ between three different régimes. When the sum of the two largest Lyapunov exponents is negative (light-gray domain), we have $`d_L<d`$ and the particles form fractal clusters in the physical space. When the sum is positive (white domain), we have $`d_L>d`$ and the particles fill the whole domain. The third case occurs only for particles lighter than the fluid and corresponds to a negative largest Lyapunov exponent (dark gray area). The particles form pointwise clusters and, when the domain is bounded, they all converge to a single trajectory.

We now present a heuristic argument, which we already gave in the case $`\beta =0`$ of heavy particles,bd02  and which predicts the threshold in Stokes number. It relies on the use of the *stability exponents*, that are the exponential growth rates of the eigenvalues of the Jacobi matrix $`𝒥_t`$ defined in (5). Using Browne’s theorem which bounds the singular values of a square matrix by its eigenvalues, a necessary condition for the fractal clustering of particles is that the sum of the $`d`$ largest stability exponents is positive. For heavy particles ($`\beta <1`$), this sum can be estimated from the local analysis of the dynamics. First, since such particles tend to cluster within the hyperbolic regions of the flow, they are spending there a fraction of time much larger than in the elliptic regions. Let us assume that the relationship between the eigenvalues of the evolution matrix $`_t`$ and those of the stress tensor $`𝝈(t)`$ can be extended to the stability exponents, at least as an approximation; we can then derive a necessary condition for the presence of fractal clusters in the physical space. For $`d=2`$, this condition can be easily written as

$$\mathrm{S}_\eta \frac{1}{\lambda _f\beta ^2}\left(\beta 2+2\sqrt{1\beta +\beta ^2}\right),$$
(7)
where $`\lambda _f`$ is the (non-dimensional) largest Lyapunov exponent associated to the carrier velocity field and calculated along the trajectory of a simple fluid particle. It was easily verified that this bound is compatible with what is observed numerically in Fig. 3.

It is often stated in the literature that the clustering of inertial particles is essentially due to the presence of coherent structures in turbulence (see., e.g., Squires and Eatonse91 ). It is indeed generally assumed that the structures with long life times appearing in the flow are responsible for a deterministic motion of the particle leading to their concentration inside or outside the vortices. Although this argument, based on the local structure of the carrier flow, allowed us to find an upper bound for the critical value of the Stokes number, it is important to stress that preferential concentrations of particles arise solely as a consequence of the dissipative character of the motion. Indeed, we have also performed simulations with carrier flows which are delta-correlated in time and that are thus completely devoid of structure.b02  These simulations show that the dependence of the Lyapunov dimension on the Stokes number $`\mathrm{S}_\eta `$ is very similar to that obtained in Fig. 1. The main difference is the behavior of $`d_L`$ as $`\mathrm{S}_\eta 0`$: in the delta-correlated case, the Lyapunov dimension tends linearly, and not quadratically, to the space dimension $`d`$. In the delta-correlated case, preferential concentrations are actually stronger than for a finite correlation time, contradicting the mechanism frequently invoked to explain concentrations.

Of course, it is not sufficient to know that the particles are concentrated in fractal objects; a fuller statistical description of their distribution is desirable. In particular, for a quantitative description of their spatial intermittency, one needs their multifractal properties (scaling exponents of the various moments of the mass contained in a sphere of radius $`r`$; see e.g. Chap. 8 of the book by Frischf95 ). Preliminary results, for heavy particles, indicate that strong spatial intermittency can be present, even when the Lyapunov dimension is just slightly below that of the physical space. Let us also mention that more quantitative results are likely to be within reach using multi-time methods in the asymptotics $`S_\eta 1`$. Small Stokes numbers are of considerable interest since in most natural and industrial situations, this is the case. A last remark concerns the extension of this approach to the case of real turbulent flows. To confirm the existence of a threshold for fractal clustering of inertial particles, one needs to resolve scales which are much smaller than the Kolmogorov dissipation scale $`\eta `$. This is quite a challenge, both for laboratory and numerical experiments.

I am deeply grateful to K. Domelevo for initiating my work on inertial particles and to U. Frisch for his continuous support. During this work, I have benefited from stimulating discussions with A. Celani, M. Cencini, P. Horvai, K. Khanin, M. Stepanov and P. Tanga. Part of this research was possible thanks to the support of the *Société de Secours des Amis des Sciences* which is warmly acknowledged. Part of this work was also supported by the National Science Foundation under agreement No. DMS-9729992 and by the European Union under contract HPRN-CT-2000-00162. Numerical simulations were performed in the framework of the SIVAM project at the Observatoire de la Côte d’Azur.


format: text

text:
# A self consistent solution to the Einstein Maxwell Dirac Equations

## I Introduction

Bound solutions of the Einstein Maxwell Dirac equations are of interest as a link between quantum mechanics and general relativity in a regime where both theories are well tested. Finster, Smoller and Yau Finster (1999a, 2000, b) along with other collaborators showed that no time periodic solutions exist for the Dirac equation in Schwarzschild, Reissner Nordstrom and Kerr Newman ($`M^2>a^2+Q^2`$) space-times. Dafermos Dafermos (2000) has confirmed this result by another method. On the other hand Schmid Schmid (2002) showed that such time periodic solutions to the Dirac equation do exist for the extreme Kerr case when $`M^2=a^2+Q^2`$.

In these investigations, it was assumed that the Dirac particle mass $`m_e`$ is very small as compared to the black hole mass $`M`$, that is, it is an external test particle. I show in the self consistent case, when the Dirac particle is the source of the curvature and the electromagnetic field, that a solution is possible for certain values of the parameters ($`m_e,Q,a`$) appearing in the Dirac equation. Further this solution is continuous everywhere and normalisable.

Such a solution corresponds to a generalisation of the free particle Dirac equation in Minkowski space to include the effects of the curvature produced by the particle. The ordinary free particle solutions of the Dirac equation are completely delocalised; the curvature however, now causes the Dirac wave functions to be localised over a region comparable in dimension to the Compton wavelength of the particle. Note that the wave function still has an enormous spread in comparison to the dimensions of the event or Cauchy horizons of the particle.

As is usually the case in single particle Dirac theory we ignore the self interaction between the particle and its own electromagnetic field. Further we assume that we need only to look for solutions with zero azimuthal angular momentum about the rotation axis. In the next section the equations resulting from Chandrasekhar’s separation of the Dirac equation in a Kerr Newman spacetime Chandrasekhar (1976, 1992) are simplified by applying these assumptions.

A necessary requirement for a wave function is that it be single valued at a point. We carry out a bifurcation analysis of the Dirac equations for $`\theta `$ and $`r`$ to find their bifurcation points. The requirement that these lie outside the physical region of interest in both the coordinates $`\theta `$ and $`r`$ gives us conditions on the parameters appearing in the Dirac equation. Here we confine ourselves to the simplest case. The equations are then solved and it is found that the angular part of the solution is continuous and normalisable. The radial equation is also solved in each of the three regions (i) the exterior, (ii) the ergosphere between the event and the Cauchy horizons and (iii) the region inside the Cauchy horizon containing the singularity. The resulting radial solutions are continuous in all three regions and so a simple matching of the wave functions at the horizons gives us a unique normalisable solution, with zero energy ($`\omega =0,\lambda =0`$).

In the last section, the nature of this solution is discussed. In particular, if the fields are centered at $`r=0`$, it is not obvious what it means to have a finite probability of finding the source at $`r0`$. I suggest that this can be explained in a manner similar to the explanation for the Klein paradox of Dirac theory. Namely, in the region of the order of the Compton wavelength within which the wave function peaks, the field strengths are so large that virtual pair creation cannot be neglected. The antiparticle thus produced can recombine with the “original” particle at $`r=0`$ leaving us to detect the particle now at $`r0`$. Consider next the case of the only allowed energy eigenvalue zero. For a single bound particle in a potential, the energy can of course be negative or zero.

## II Solution of the Dirac equation

Following the notation of Chandrasekhar’s book Chandrasekhar (1992), the Kerr Newman line element in Boyer Lindquist coordinates is given by

$`(ds)^2`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Sigma }}}(dta\mathrm{sin}^2\theta d\varphi )^2\mathrm{\Sigma }({\displaystyle \frac{dr^2}{\mathrm{\Delta }}}+d\theta ^2){\displaystyle \frac{\mathrm{sin}^2\theta }{\mathrm{\Sigma }}}(adt(r^2+a^2)d\varphi )^2,`$ (1a)
$`\mathrm{\Sigma }`$ $`=`$ $`r^2+a^2,`$ (1b)
$`\mathrm{\Delta }`$ $`=`$ $`r^22Mr+a^2+Q^2.`$ (1c)
As we are using natural units, the mass appearing in the metric $`(M\frac{GM}{c^2})`$ has a size comparable to the horizon length. Whereas, the mass $`m_e`$ which appears in the Dirac equation is $`(m_e\frac{m_ec}{\mathrm{}})`$ has a size which is the inverse of the Compton wavelength. Thus we use different symbols to distinguish the two length scales in the problem. They both of course refer to the same mass but differ according to which of the universal constants $`G,\mathrm{},c`$ was being suppressed.

The Dirac equation can then be separated by writing the four components of the spinor wave function in the form

$`\mathrm{\Psi }_1`$ $`=`$ $`{\displaystyle \frac{R_{\frac{1}{2}}S_{\frac{1}{2}}}{ria\mathrm{cos}\theta }}e^{i(m\varphi \omega t)},`$ (2a)
$`\mathrm{\Psi }_2`$ $`=`$ $`R_{\frac{1}{2}}S_{\frac{1}{2}}e^{i(m\varphi \omega t)},`$ (2b)
$`\mathrm{\Psi }_3`$ $`=`$ $`R_{\frac{1}{2}}S_{\frac{1}{2}}e^{i(m\varphi \omega t)},`$ (2c)
$`\mathrm{\Psi }_4`$ $`=`$ $`{\displaystyle \frac{R_{\frac{1}{2}}S_{\frac{1}{2}}}{r+ia\mathrm{cos}\theta }}e^{i(m\varphi \omega t)}.`$ (2d)
I have slightly modified the notation to use $`\omega `$ for the frequency rather than $`\sigma `$. The $`R_{\pm \frac{1}{2}}`$ are functions of $`r`$ only while $`S_{\pm \frac{1}{2}}`$ are functions of $`\theta `$ only. If this form is substituted into the Dirac equations, Chandrasekhar (1992) we get the following angular equations.

$`{\displaystyle \frac{d}{d\theta }}S_{+\frac{1}{2}}`$ $`=`$ $`(m\mathrm{csc}\theta +{\displaystyle \frac{1}{2}}\mathrm{cot}\theta )S_{+\frac{1}{2}}(\lambda am_e\mathrm{cos}\theta )S_{\frac{1}{2}},`$ (3a)
$`{\displaystyle \frac{d}{d\theta }}S_{\frac{1}{2}}`$ $`=`$ $`(m\mathrm{csc}\theta {\displaystyle \frac{1}{2}}\mathrm{cot}\theta )S_{\frac{1}{2}}+(\lambda +am_e\mathrm{cos}\theta )S_{+\frac{1}{2}}.`$ (3b)
We now set $`m=0`$ and carry out a bifurcation analysis of these equations. The angle equations (3) will simplify under the substitutions

$$T_{\pm \frac{1}{2}}=\sqrt{\mathrm{sin}\theta }S_{\pm \frac{1}{2}},$$
(4)
to yield

$`{\displaystyle \frac{d}{d\theta }}T_{+\frac{1}{2}}`$ $`=`$ $`(\lambda am_e\mathrm{cos}\theta )T_{\frac{1}{2}},`$ (5a)
$`{\displaystyle \frac{d}{d\theta }}T_{\frac{1}{2}}`$ $`=`$ $`(\lambda +am_e\mathrm{cos}\theta )T_{+\frac{1}{2}}.`$ (5b)
These have bifurcations at

$$(\lambda )^2=(am_e\mathrm{cos}\theta )^2.$$
(6)
Our condition that the wave functions be single valued means that we wish to have no bifurcations in the range $`0\theta <\pi `$. This is possible if and only if $`|\lambda |am_e`$ or $`\lambda =0`$. The second case corresponds to an s wave orbital state in the usual separation for the Dirac equation Greiner (1992) in a central potential. I have only considered this case here. The angle equations then are

$`{\displaystyle \frac{d}{d\theta }}T_{+\frac{1}{2}}`$ $`=`$ $`am_e\mathrm{cos}\theta T_{\frac{1}{2}},`$ (7a)
$`{\displaystyle \frac{d}{d\theta }}T_{\frac{1}{2}}`$ $`=`$ $`am_e\mathrm{cos}\theta T_{+\frac{1}{2}}.`$ (7b)
These are trivial to solve and we obtain

$`S_{\frac{1}{2}}`$ $`=`$ $`A(e^{am_e\mathrm{sin}\theta }\pm e^{am_e\mathrm{sin}\theta })/(\sqrt{2\mathrm{sin}\theta }),`$ (8a)
$`S_{\frac{1}{2}}`$ $`=`$ $`A(e^{am_e\mathrm{sin}\theta }e^{am_e\mathrm{sin}\theta })/(\sqrt{2\mathrm{sin}\theta }).`$ (8b)
These solutions while they have a formal singularity at the poles give a finite probability density every where.

Next consider the radial equations. Again I find it convenient to use a slightly different variable from these used by Chandrasekhar (Chandrasekhar (1992)), we define $`R_{}^{^{}}{}_{\frac{1}{2}}{}^{}=\sqrt{\mathrm{\Delta }}R_{\frac{1}{2}}`$

$`{\displaystyle \frac{d}{dr}}R_{}^{^{}}{}_{+\frac{1}{2}}{}^{}(r)`$ $`=`$ $`i{\displaystyle \frac{\omega (r^2+a^2)+\lambda a}{\mathrm{\Delta }}}R_{}^{^{}}{}_{+\frac{1}{2}}{}^{}(r)+{\displaystyle \frac{\lambda im_er}{\sqrt{\mathrm{\Delta }}}}R_{\frac{1}{2}}(r),`$ (9a)
$`{\displaystyle \frac{d}{dr}}R_{\frac{1}{2}}(r)`$ $`=`$ $`i{\displaystyle \frac{\omega (r^2+a^2)+\lambda a}{\mathrm{\Delta }}}R_{\frac{1}{2}}(r)+{\displaystyle \frac{\lambda +im_er}{\sqrt{\mathrm{\Delta }}}}R_{}^{^{}}{}_{+\frac{1}{2}}{}^{}(r).`$ (9b)
The corresponding bifurcation points lie at those values of $`r`$ which satisfy

$$(\omega (r^2+a^2)+\lambda a)^2=(\lambda ^2+m_er^2)\mathrm{\Delta }.$$
(10)
This fourth order equation permits a wide range in the 5 dimensional parameter space $`(m_e,\lambda ,Q,\omega ,m)`$ of the Dirac equation for which the bifurcations lie outside the region of interest. It is worth remarking that such a solution must always have $`\omega ^2m_{e}^{}{}_{}{}^{2}`$. Further the solutions with $`\omega =\pm m_e`$ are unconditionally stable. As we are considering only the case of no orbital angular momentum, $`\lambda =0`$, equation (10) also permits solutions with $`\omega =0`$ which have bifurcations at $`r=M\pm \sqrt{\mathrm{\Delta }}`$ and $`r=0`$. I will consider only this last case ($`\omega =0`$) in this work. For $`\omega =0`$ and $`\lambda =0`$, the radial equations (9) simplify to yield

$`{\displaystyle \frac{d}{dr}}R_{}^{^{}}{}_{+\frac{1}{2}}{}^{}(r)`$ $`=`$ $`{\displaystyle \frac{im_er}{\sqrt{\mathrm{\Delta }}}}R_{\frac{1}{2}}(r),`$ (11a)
$`{\displaystyle \frac{d}{dr}}R_{\frac{1}{2}}(r)`$ $`=`$ $`{\displaystyle \frac{im_er}{\sqrt{\mathrm{\Delta }}}}R_{}^{^{}}{}_{+\frac{1}{2}}{}^{}(r).`$ (11b)
Again as in the case of the angle equations, these can be solved quite easily. Consider first the exterior region $`r>M+\sqrt{\mathrm{\Delta }}`$ if we change the independent variable to $`u=\sqrt{\mathrm{\Delta }}+M\mathrm{ln}(rM+\sqrt{\mathrm{\Delta }})`$ the radial equations (11) simplify to yield

$`{\displaystyle \frac{d}{du}}R_{}^{^{}}{}_{+\frac{1}{2}}{}^{}(u)`$ $`=`$ $`im_eR_{\frac{1}{2}}(u),`$ (12a)
$`{\displaystyle \frac{d}{du}}R_{\frac{1}{2}}(u)`$ $`=`$ $`im_eR_{}^{^{}}{}_{+\frac{1}{2}}{}^{}(u).`$ (12b)
We can thus write down the solution in the exterior region which is finite as $`r\mathrm{}`$ as

$`\sqrt{\mathrm{\Delta }}R_{\frac{1}{2}}`$ $`=`$ $`B{\displaystyle \frac{e^{m_e\sqrt{\mathrm{\Delta }}}}{[rM+\sqrt{\mathrm{\Delta }}]^{m_eM}}},`$ (13a)
$`R_{\frac{1}{2}}`$ $`=`$ $`iB{\displaystyle \frac{e^{m_e\sqrt{\mathrm{\Delta }}}}{[rM+\sqrt{\mathrm{\Delta }}]^{m_eM}}}.`$ (13b)
Next consider the ergospheric region between the event horizon $`r_+=M+\sqrt{\mathrm{\Delta }}`$ and the Cauchy Horizon $`r_{}=M\sqrt{\mathrm{\Delta }}`$. Here we transform the independent variable to $`v=i\sqrt{|\mathrm{\Delta }|}+M\mathrm{ln}(rM+i\sqrt{|\mathrm{\Delta }|})`$ and obtain the solution

$`\sqrt{\mathrm{\Delta }}R_{\frac{1}{2}}`$ $`=`$ $`Ce^{im_e\sqrt{|\mathrm{\Delta }|}}[rM+i\sqrt{|\mathrm{\Delta }|}]^{m_eM}+D{\displaystyle \frac{e^{im_e\sqrt{|\mathrm{\Delta }|}}}{[rM+i\sqrt{|\mathrm{\Delta }|}]^{m_eM}}},`$ (14a)
$`R_{\frac{1}{2}}`$ $`=`$ $`iCe^{im_e\sqrt{|\mathrm{\Delta }|}}[rM+i\sqrt{|\mathrm{\Delta }|}]^{m_eM}iD{\displaystyle \frac{e^{im_e\sqrt{|\mathrm{\Delta }|}}}{[rM+i\sqrt{|\mathrm{\Delta }|}]^{m_eM}}}.`$ (14b)
Finally in the region interior to the Cauchy horizon $`rr_{}`$ we use the new independent variable $`w=\sqrt{\mathrm{\Delta }}+M\mathrm{ln}(rM+\sqrt{\mathrm{\Delta }})`$ to obtain the solution

$`\sqrt{\mathrm{\Delta }}R_{\frac{1}{2}}`$ $`=`$ $`Ee^{m_e\sqrt{\mathrm{\Delta }}}[rM+\sqrt{\mathrm{\Delta }}]^{m_eM}+F{\displaystyle \frac{e^{m_e\sqrt{\mathrm{\Delta }}}}{[rM+\sqrt{\mathrm{\Delta }}]^{m_eM}}},`$ (15a)
$`R_{\frac{1}{2}}`$ $`=`$ $`iEe^{m_e\sqrt{\mathrm{\Delta }}}[rM+\sqrt{\mathrm{\Delta }}]^{m_eM}+iF{\displaystyle \frac{e^{m_e\sqrt{\mathrm{\Delta }}}}{[rM+\sqrt{\mathrm{\Delta }}]^{m_eM}}}.`$ (15b)
Thus the solutions are continuous in all three regions, the only effect of the bifurcation at $`r_+`$ and $`r_{}`$ being an ambiguity as to whether it is my choice in eq.(14 which is correct or its complex conjugate which is to be used in this region. As the solution is continuous in all three regions, the coefficients $`B,C,D,E,F`$ can be found up to a normalisation constant by requiring that the solutions be continuous at the horizons. So we finally obtain the solution

$`\sqrt{\mathrm{\Delta }}R_{\frac{1}{2}}`$ $`=`$ $`B{\displaystyle \frac{e^{m_e\sqrt{\mathrm{\Delta }}}}{[rM+\sqrt{\mathrm{\Delta }}]^{m_eM}}},(rr_+)`$ (16a)
$`\sqrt{\mathrm{\Delta }}R_{\frac{1}{2}}`$ $`=`$ $`B{\displaystyle \frac{e^{im_e\sqrt{|\mathrm{\Delta }|}}}{[rM+i\sqrt{|\mathrm{\Delta }|}]^{m_eM}}},(r_{}rr_+)`$ (16b)
$`\sqrt{\mathrm{\Delta }}R_{\frac{1}{2}}`$ $`=`$ $`Be^{m_e\sqrt{\mathrm{\Delta }}}[rM+\sqrt{\mathrm{\Delta }}]^{m_eM},(rr_{})`$ (16c)
$`R_{\frac{1}{2}}`$ $`=`$ $`iB{\displaystyle \frac{e^{m_e\sqrt{\mathrm{\Delta }}}}{[rM+\sqrt{\mathrm{\Delta }}]^{m_eM}}},(rr_+)`$ (16d)
$`R_{\frac{1}{2}}`$ $`=`$ $`iB{\displaystyle \frac{e^{im_e\sqrt{|\mathrm{\Delta }|}}}{[rM+i\sqrt{|\mathrm{\Delta }|}]^{m_eM}}}(r_{}rr_+)`$ (16e)
$`R_{\frac{1}{2}}`$ $`=`$ $`iBe^{m_e\sqrt{\mathrm{\Delta }}}[rM+\sqrt{\mathrm{\Delta }}]^{m_eM}(rr_{}).`$ (16f)
## III Discussion

The solution we have just obtained is normalisable and thus represents a proper Dirac spinor wave function. There is of course a second possibility due to the bifurcations at the horizon, which is to take the complex conjugate form inside the ergosphere. This however will give rise to the same Dirac currents etc., outside and is thus indistinguishable form the form presented here.

Next consider the question of the only allowed energy eigenvalue, zero. In any bound state problem, we expect the energy of the state to be less than or equal to the asymptotic value of the binding potential. In this case if we consider the curvature to be producing the attractive force localising the particle close to $`r0`$; then the zero energy eigenvalue should evoke no surprise.

Finally there is the question of what is meant by finding the source of the field at a point other than where the field originates from. Again if this is considered as a potential problem for the Dirac field this can be easily explained. The wave function given by eq.(16) is localised over a region of the order of the Compton wavelength of the particle. In this region the field strengths are very high and so virtual pairs of the same particle type have finite amplitude of appearing. The antiparticle of such a pair can then combine with the “ original” particle leaving the now real particle of the pair to be found in the position $`r0`$.

While the present solution might appear overly simplistic, it is the fact that it exists which is important. This suggests that other such solutions may also be found in the rest of the allowed parameter space which was not explored here.


text:
# IGUANA Architecture, Framework and Toolkit for Interactive Graphics

## I INTRODUCTION

IGUANA - Interactive Graphics for User ANAlysis - is a modular C++ toolkit for interactive visualisation. The project mainly focuses on interactive detector and event visualisation with emphasis on:

* High performance 2D and 3D graphics;
* Graphical user interfaces;
* User access to the experiment services: data access framework, application execution framework, etc.

The goal of the project is to provide easy-to-use coherent interactive graphical application interface for the physicist, where the main application - IGUANA studio - interfaces to other tools and components. Interactive histogram, ntuple analysis is not considered a primary goal. It is assumed that this functionality will be provided by other tools and projects such as JAS Jas , Hippodraw hippo-01 , ROOT Root , or OpenScientist Opensci , which may (or may not) ultimately be integrated with IGUANA.

IGUANA development is mainly driven by the needs of the CMS and D0 experiments, but most of IGUANA is independent of CMS or D0 and can be freely used by others. In CMS, IGUANA releases are available directly via SCRAM scram  which keeps a list of the IGUANA project releases. Instructions versioned by release on how to download and install IGUANA are available for any selected version.

IGUANA is currently supported on Linux, Solaris and Windows. Porting to other unix varieties is reasonably straightforward.

## II ARCHITECTURE AND FRAMEWORK

IGUANA architecture is designed to provide a generic tool that can integrate with a specific experiment or even a single task within it without having to specify a fixed list of abstract interfaces or having to standardise on a common object description. IGUANA has a small kernel, everything else is implemented in demand-loaded extensions which negotiate the specific services they need. This allows new services to be added quickly, to be propagated into existing code gradually, and existing services to be altered to do new things easily.

The full documentation is available on the project website111http://iguana.cern.ch automatically generated for every release from the sources kept in CVS.

### II.1 Core Overview

The main units relating to the IGUANA architecture are:

* A thin portability and utilities layer;
* A small kernel that manages a number of plug-ins:

+ application personalities;
    + session with extensions forming the shared application state;
    + user interface components: sites, browsers;
    + representation methods to map between the experiment objects and the various browsers;
* External software imported into IGUANA for convenience of building and distribution, and external software which remains outside IGUANA.

The portability layer, is mostly implemented by the classlib package which in future will be distributed by the LHC Computing Grid (LCG) SEAL project Seal .

The kernel is implemented by IgPluginManager, a generic plug-in manager package, which is extended by the architecture kernel package IgObjectBrowser that provides the core functionality: session objects, extensions, user interface components, representations, etc. and base interfaces for the various plug-in types.

The plug-ins are in the Ig\_Modules subsystem which includes the event display core, application management services, and thin adaptors to export external software in forms understood by the architecture. For more details see arch-01 .

### II.2 Application Personalities

When an IGUANA program runs, it first creates a session object into which it then attaches an application personality: the main program that determines all subsequent behaviour (IgSession, IgDriver). Typically the personality immediately extends the session object with services pertinent to the purpose of the application (IgExtension). For example, a graphical personality would create a main window and add services that give access to the menus and maintain GUI event loops. The personality then loads a number of extensions into the session; a graphical personality would tell the system to load all “GUI extensions” (IgExtension).

The personality and the extensions form the application. The personality exposes features by installing service interfaces into the session, based on which other extensions can provide further services and make features available to the user.

### II.3 Browsers and Sites

An interactive application personality creates its user interface based on one of the user interface components: browsers and sites (IgBrowser, IgSite). A browser is a way to look at and interact with application objects (see Figure 1). It may also control other things, such as a cut in a histogram display. A browser does not however have to be graphical: it might just sit on the background and respond to other browsers’ queries such as “What actions can be invoked on this object?” Sites host browsers, by providing for example a window for a 3D browser view. More generally sites expose objects such as GUI widgets as hosts for other sites and browsers. Sites need not to be graphical either: a pipe site might host a browser communicating to a subprocess.

### II.4 Object Representations

Application objects can be represented in many different browsers. A browser typically has a model made of a number of object reps (IgModel, IgRep). The use of model-specific reps, such as 3D shape objects, is encouraged rather then the application objects themselves. This permits an object’s representation to be separated from the object itself. It is possible to correlate the object’s reps in different browsers and to create a rep in a right context in each browser.

This implementation requires “representable” application objects to inherit from a common base type (IgRepresentable) or make it an alias for another class. The only requirement is that the type must be polymorphic; IgRepresentable simply defines a virtual destructor and no other methods.

In CMS there is however no base type shared by all objects; for visualisation the proxies are used that inherit IgRepresentable and point to other CMS objects. The primary reason for this choice is the framework’s support for “virtual objects” that can be materialized upon request.

The object – rep mappings are extensions loaded dynamically on demand. IGUANA discovers and chooses the right mapping automatically; the extension code simply does whatever is appropriate for that combination. It is not necessary to have a global list of all the pairings: code for new reps and views can be dropped in without any changes to existing code.

### II.5 Communication

Communication takes place through three channels: session extensions, message services, and the object mapping methods. The first has already been covered. The message services allow browsers to share messages such as “I selected this”: all observers of a service can maintain together a coherent state and to answer queries from each other while still knowing next to nothing about the message sender (see Figure 1).

The final channel is the object–rep–model operations, one of which was already mentioned: the creation of a rep. Another operation commits a changed rep back to the object, yet another lazily expands an object rep (for example to read data only when it is requested). The operations can be extended, both to handle new argument combinations and with new methods with arbitrary parameter types (IgBrowserMethods).

## III VISUALISATION FRAMEWORK FOR GEANT4

IGUANA provides an interactive visualisation sub-framework for GEANT4 Geant4 . The framework implements DCUT, DTREE-like functionality (as in Geant3). It allows to explore and visualise the volume tree, with all the usual IGUANA 3D features: view, correlated picking, slices per object, etc. with the Geant4 command line still accessible. Navigation in the volume tree can be done either by logical or physical volumes, or subsets. There are quick operations for common tasks. Volume property window shows additional information about selected object or volume. Volume tree selectors provide various types of selection:

* by material: ”show all silicon”;
* by sensitive: ”show only sensitive detectors”;
* for a sub-tree, predefined viewpoints/settings;
* forward + reverse: ”show where this is used”.

The framework is used in the CMS OSCAR simulation project where an IGUANA wizard is used to guide through GEANT4/OSCAR settings. The application is integrated with Martin Liendl’s overlap detection: overlaps are found and the result details are shown in a list.

Two generic examples showing how to use an arbitrary GEANT4 geometry with the framework are distributed with IGUANA experiment-independent core: the GEANT4 example number two and the ATLAS calorimeter example.

## IV APPLICATIONS

There are several applications based on IGUANA plug-in architecture. Among those are OSCAR and ORCA visualisation which are well known in CMS and actively used for visualising simulated and reconstructed CMS data. D0Scan (see Figure 2) currently uses IGUANA as a toolkit. It is used at D0 for event and detector visualisation.

### IV.1 Experiment Specific Plug-ins

IGUANA is used by other software projects to define experiment-specific views. Among these projects are the CMS experiment projects222http://cmsdoc.cern.ch/cmsoo/cmsoo.html: COBRA - core framework, OSCAR - simulation project, ORCA orca  \- reconstruction project, FAMOS - fast simulation project, DDD - detector geometry description and the D0 experiment project: D0Scan.

The integration with experiment-specific frameworks is an IGUANA dependent part, usually a sub-system, which is located within the experiment project. Even though IGUANA core knows nothing about these sub-systems, the main driver can still load experiment-specific plug-ins providing that the location of those plug-ins is listed in the IGUANA\_PLUGINS environment variable as path defined by a colon-separated list of directories and the experiment-specific shared libraries are available at run-time.

The projects build plug-ins similar to the described above (section II.1). Usually, each project has a dedicated sub-system called Visualisation for that purpose. The sub-system contains project specific modules, or plug-ins. The kernel loads the modules dynamically as required. Currently a user selects what to load at the interactive setup phase or via a resource file.

The plug-ins are normally demand-loaded shared libraries but IGUANA also supports statically linked applications, therefore no shared libraries need to be involved.

### IV.2 ORCA Visualisation

Visualisation for the CMS reconstruction project (ORCA) (see Figure 1 and Figure 3) provides access via COBRA plug-ins to both simulated (cmsim) and reconstructed event data, and Geant3 based CMS geometry as well as reconstructed geometry. Events are requested interactively one by one from a specified data collection. The IGUANA GUI runs in a separate thread. Communication between the GUI thread and the event thread is handled by the COBRA-based event processor plug-in.

The visualisation content is defined at run-time. Plug-ins for different ORCA sub-systems can be loaded independently. A user can define the content of the application in an ASCII text resource file. Various visualisation cuts can also be defined in the resource file. The on-line configuration editor implemented as an IGUANA service allows to change them interactively or set new configurables during the run.

### IV.3 D0Scan

D0Scan is an event viewer for D0 using IGUANA. The D0Scan design is based upon the IGUANA IgVis module which provides a very flexible set of tools for organizing an OpenInventor scene, creating a GUI to control visibility of items in the scene, and centralizing picking callback. Items such as detector components, hits, tracks, and so forth are implemented as IgQtTwigs and arranged in a hierarchical structure similar to the OpenInventor scene graph. This structure is represented in the Twig browser that controls visibility of graphical objects.

D0Scan makes heavy use of the ability of IgVis to support multiple windows with different characteristics. User demand has led to the development of specialized displays as, for instance, in the $`\eta \varphi `$ projection of the distribution of energy in the calorimetry.

### IV.4 OSCAR Visualisation

Visualisation for CMS simulation project using GEANT4 - OSCAR visualisation - is based on the interactive visualisation sub-framework of IGUANA for GEANT4 described above in section III. CMS detector geometry is described in XML by DDD - Detector Description Database project. This description is converted to the GEANT4 one and visualised.

## V TOOLKIT

IGUANA makes use of suitable existing software where possible and integrates the whole in a coherent fashion. Application developers select from the toolkit only those parts which are relevant for a particular application.

External software used by IGUANA includes the 3D graphics libraries: OpenGL and OpenInventor, the Qt toolkit as a GUI kernel, public-domain Qt and Inventor extensions, various packages for plotting and interactive analysis, an XML parser, etc.

Being accessible via IGUANA services with dedicated GUI extensions, both Oprofile and Jprof are used for debugging. IgNominy Ignominy  \- an IGUANA tool for dependency analysis depends on Graphviz for drawing diagrams. Doxygen is used for auto generated reference documentation. IGUANA studio printing services integrate gl2ps for producing high quality vector postscript. The team has improved culling algorithms of the package which reduce the size of the printed file dramatically and fed them back to the author.

### V.1 LCG Software and Services

IGUANA team actively participates in the LCG project. Part of the software has already migrated to SEAL such as the classlib and plug-in manager. IGUANA uses services provided by SPI such as Savannah bug reporting system. POOL is used indirectly via experiment specific plug-ins.

## VI FUTURE PLANS

IGUANA will provide wider selection of 2D and 3D representations and specialized viewers including the time snapshots of visualized data, algorithm dependent views, and animation.

The project progresses toward coherent physicist desktop emphasizing the role of wizards for experiment specific environment, control center, and other services.

Future possibilities are open for discussion in context of the LHC LCG project, e.g. PI - Physicist Interface project. Among those there are issues of integration and future relationship with ROOT, JAS, etc. and potential front-end GUI for GRID applications.

###### Acknowledgements.

This work is supported by the US National Science Foundation.


text:
# Heterogeneity in oscillator networks: Are smaller worlds easier to synchronize?

Phys. Rev. Lett. 91, 014101 (2003)

Current address: \] Department of Mathematics, 208 Clements Hall, PO Box 750156, Southern Methodist University, Dallas, TX 75275-0156.

Current address: \] Max Planck Institute for the Physics of Complex Systems Nothnitzer Strasse 38, 01187 Dresden, Germany.

## Abstract

Small-world and scale-free networks are known to be more easily synchronized than regular lattices, which is usually attributed to the smaller network distance between oscillators. Surprisingly, we find that networks with a homogeneous distribution of connectivity are more synchronizable than heterogeneous ones, even though the average network distance is larger. We present numerical computations and analytical estimates on synchronizability of the network in terms of its heterogeneity parameters. Our results suggest that some degree of homogeneity is expected in naturally evolved structures, such as neural networks, where synchronizability is desirable.

In their seminal work, Watts and Strogatz Watts and Strogatz (1998); Watts (1999) have shown that real networks of different nature have the small-world (SW) property, characterized by high clustering and the average network distance between two nodes that is as small as pure random networks. Since then, many more examples of real-world networks Albert et al. (1999); Amaral et al. (2000); Strogatz (2001); Newman (2001a); Montoya and Solé (2002); Motter et al. (2002a), including both artificial and natural systems, have been identified to have the SW property. It turns out, however, that there is another seemingly generic feature of networks in the real world. It is called the scale-free (SF) property, which is signified by the power-law connectivity distribution of the network Barabási and Albert (1999); Albert et al. (1999); Faloutsos et al. (1999); Amaral et al. (2000); Newman (2001b); Liljeros et al. (2001); Solé and Montoya (2001); Albert and Barabási (2002); Motter et al. (2002a). Barabási and Albert Barabási and Albert (1999) suggested a model of growing networks, in which preferential attachment of new links to nodes with higher connectivity results in the SF property. There has also been some efforts to incorporate both SF and SW properties in a single model Amaral et al. (2000); Newman et al. (2001); klemm2002 .

So far, much research has been focused on the structural properties of SF and SW network models. Despite the widespread belief that these structural properties must have significant impact on dynamical processes taking place on such networks Strogatz (2001), there has been little work directly addressing this issue. Most work has dealt with synchronization of oscillators whose topology of interaction has either the SF or SW Watts (1999); Lago-Fernández et al. (2000); Gade and Hu (2000); Jost and Joy (2001); Hong et al. (2002); Wang and Chen (2002a); Barahona and Pecora (2002); kwon2002  property, showing that it leads to improved synchronizability when compared to local lattice topology. A general argument underlying this phenomenon is that communication between oscillators are more efficient because of the smaller average network distance. But, does smaller average network distance improve synchronizability?

For many real networks, heterogeneity is a common trait which frequently manifests itself in the form of an SF distribution of connectivities. It is known that such heterogeneity tends to reduce the average network distance CohenBollobas , and this leads naturally to the question of whether heterogeneity improves synchronizability.

The aim of this paper is to demonstrate, by using important classes of SF and SW network models, that heterogeneity of the connectivity distribution causes the opposite to hold; namely, as heterogeneity increases, the average network distance is reduced but synchronization becomes more difficult to achieve. We show that this intriguing behavior can be explained by examining the load distribution on nodes or links, where the load of a node (or a link) quantifies the traffic of communication passing through it. The analytical results we derive suggest that our observations are quite general and are expected to hold for a wide class of complex networks.

Synchronizability of a network of oscillators can be quantified through the eigenvalue spectrum of the Laplacian matrix representing the connection topology of the network. Here we follow the general framework established in Barahona and Pecora (2002); Pecora and Carroll (1998). Consider a network of $`N`$ identical dynamical systems with symmetric coupling between oscillators. The equations of motion for the system are

$$\dot{x_i}=F(x_i)+\sigma \underset{j=1}{\overset{N}{}}L_{ij}H(x_j),i=1,\mathrm{},N,$$
(1)
where $`\dot{x}=F(x)`$ governs the dynamics of each individual node, $`H(x)`$ is the output function, $`\sigma `$ is the overall strength of coupling, and $`L`$ is the Laplacian matrix, defined to be $`L_{ij}=1`$ if nodes $`i`$ and $`j`$ are connected, $`L_{ii}=k_i`$ if node $`i`$ is connected to $`k_i`$ other nodes, and $`L_{ij}=0`$ otherwise. The linear stability of the synchronized state $`\{x_i(t)=x^{}(t),i\}`$ is determined by the corresponding variational equations, which can be diagonalized into $`N`$ blocks of the form $`\dot{y}=[DF(s)+\lambda DH(s)]y`$, where $`y`$ represents different modes of perturbation from the synchronized state. We have $`\lambda =\sigma \lambda _i`$ for the $`i`$th block, $`i=1,2,\mathrm{},N`$, and $`\lambda _1=0\lambda _2\mathrm{}\lambda _N`$ are the eigenvalues of $`L`$ footnote1 . The largest Lyapunov exponent $`\mathrm{\Lambda }(\lambda )`$ for this equation, also called the master stability function Pecora and Carroll (1998), determines the linear stability of the synchronized state for *any* linear coupling scheme. In particular, the synchronized state is stable if $`\mathrm{\Lambda }(\sigma \lambda _i)<0`$, for each $`i=2,\mathrm{},N`$ footnote2 . It was found Pecora and Carroll (1998) that for a large class of chaotic oscillatory systems, there exists a single parameter interval $`(\alpha _1,\alpha _2)`$ on which $`\mathrm{\Lambda }(\lambda )<0`$. In this case, there is a value of the coupling strength $`\sigma `$ for which the synchronized state is linearly stable, if and only if $`\lambda _N/\lambda _2<\alpha _2/\alpha _1\beta `$, where $`\beta `$ is a constant that depends on $`F(x)`$, $`x^{}(t)`$, and $`H(x)`$, but not on $`L`$. The value of $`\beta `$ ranges from 5 to 100 for various chaotic oscillators Barahona and Pecora (2002). Realizing that the ratio $`\lambda _N/\lambda _2`$ depends only on the topology of interactions among oscillators, we see that the impact of having a particular coupling topology on the network’s ability to synchronize is represented in a single quantity $`\lambda _N/\lambda _2`$: the larger the ratio, the more difficult it is to synchronize the oscillators, and vice versa Barahona and Pecora (2002).

Having reduced the problem of synchronizability to finding eigenvalues of the Laplacian matrix $`L`$, we now study the effect of heterogeneity in the connectivity distribution. We first consider the semirandom model of SF networks Newman et al. (2001), where the connectivity $`k_ik_0`$ of each node $`i`$ is chosen at random according to the probability distribution $`P(k)k^\gamma `$ with the scaling exponent $`\gamma `$. An SF network is then generated by randomly connecting nodes, forcing each node $`i`$ to have connectivity $`k_i`$, and prohibiting self- and repeated links Newman et al. (2001); motter2002b .

The average network distance $`\overline{𝒟}`$ of an SF network decreases with increasing heterogeneity of the connectivity distribution, where $`\overline{𝒟}`$ is defined as the minimum number of links that must be followed to go from one node to another, averaged over all pairs of nodes. Figure 1(a) shows the dependence of $`\overline{𝒟}`$ on the scaling exponent $`\gamma `$ for the semirandom model. In order to quantify the heterogeneity of the connectivity distribution, we plot its standard deviation $`s`$ in the inset of Fig. 1(a). We observe that smaller values of $`\gamma `$ result in a longer tail in the connectivity distribution, which in turn makes $`s`$ larger. Note that for completely homogeneous case of $`\gamma =\mathrm{}`$, which corresponds to taking $`k_i=k_0`$ for all $`i`$, we have $`s=0`$. Note also that decreasing $`\gamma `$ increases the average connectivity $`\overline{k}`$ as well, which can be seen in the same inset. This implies that the total number of connections increases. In short, the more heterogeneous the network is, the more connected and “smaller” it is.

Since the SF and SW networks are known to show enhanced synchronizability over lattices, one might expect that an SF network with smaller $`\overline{𝒟}`$ would have improved synchronizability. However, we observe the opposite phenomenon: the ratio $`\lambda _N/\lambda _2`$ increases as $`\gamma `$ is increased, as shown in Fig. 1(b). In other words, as the network becomes more heterogeneous, it becomes less synchronizable, even though $`\overline{k}`$ gets larger and $`\overline{𝒟}`$ gets smaller footnote3 .

The heuristic reason for this surprising behavior lies in the fact that a few “center” oscillators interacting with a large number of other oscillators tend to get overloaded by the traffic of communication passing through them. When too many independent signals with different phases and frequencies are going through a node at the same time, they can cancel one another, resulting in effectively no communication between oscillators. Hence, the more concentrated the traffic is on a few nodes, the more difficult it is to achieve efficient communication between oscillators, leading to reduced synchronizability. The same effect results from overloading links as well.

The amount of communication traffic passing through node $`i`$ can be quantified by the load on $`i`$, defined as the number of shortest paths between two (other) nodes that pass through $`i`$ Newman (2001b); Goh et al. (2001). To quantify the extent to which the load distribution is concentrated on a few nodes, we first normalize the load on each node by the total load of the network, and then focus our attention on its maximum value $`\mathrm{}_{\mathrm{max}}`$ over all nodes. The dependence of $`\mathrm{}_{\mathrm{max}}`$ on $`\gamma `$ is shown in the inset of Fig. 1(b). We see that the behavior of $`\mathrm{}_{\mathrm{max}}`$ is followed closely by the ratio $`\lambda _N/\lambda _2`$. This is consistent with the result in Goh et al. (2001) for another model of SF networks. It is shown there that $`\mathrm{}_{\mathrm{max}}`$ for the largest hub scales as $`N^\delta `$, where the scaling exponent $`\delta 0.2`$ for $`2<\gamma 3`$, but $`\delta `$ gets larger as $`\gamma >3`$ is increased. This means that for a fixed $`N`$, more heterogeneous connectivity distribution implies more heterogeneous load distribution on nodes.

To illustrate the prevalence of this phenomenon, we next consider a growing model of SF networks with aging of nodes dorogovtsev2000b . In this model, the network initially has $`n_0`$ nodes. At each time step, a new node with $`m_a`$ links is added. The probability for an existing node $`i`$ to receive a link from the new node is proportional to $`k_i`$ and $`\tau _i^\alpha `$, where $`\tau _i`$ is the age of node $`i`$, or the number of time steps since it was added to the network. It was shown in dorogovtsev2000b  that the parameter $`\alpha `$ determines the scaling exponent $`\gamma `$ for the resulting SF network, and that smaller $`\alpha `$ $`(1)`$ yields smaller $`\gamma `$ $`(2)`$. Figures 1(c,d) show the result for the growing model, from which one observes the same phenomenon of reduced synchronizability due to heterogeneity.

To show that any type of heterogeneity induces the same behavior, we now introduce a variant of the SW model Watts and Strogatz (1998), which is constructed as follows. First we start with a regular one-dimensional lattice having $`N`$ nodes with connections up to $`\kappa `$th nearest neighbors. We next choose $`n_c`$ nodes at random from all nodes with equal probabilities, and assign them to be centers, or hubs. Finally, we add each of $`m`$ shortcuts by connecting one node chosen at random from all $`N`$ nodes to another node randomly chosen from the $`n_c`$ center nodes. We call this the two-layer SW model since nodes can be naturally divided into two groups: centers, which tend to have higher connectivity, and the others that have lower connectivity. Thus, the heterogeneity is controlled by the parameter $`n_c`$: small $`n_c`$ leads to higher connectivity of the centers, which in turn results in increased heterogeneity.

Two extreme cases in our model exhibit interesting features. When $`n_c=N`$, shortcuts are simply added uniformly at random, effectively reducing the network to the homogeneous SW model of Watts and Strogatz Watts and Strogatz (1998); Newman et al. (2000). If $`n_c=1`$, all shortcuts are connected to a single center, making the configuration close to the one which minimizes the average network distance $`\overline{𝒟}`$ dorogovtsev ; nishikawa2002 .

The network exhibits a stronger SW property when $`n_c`$ is smaller, i.e., when the network is heterogeneous. As can be seen in Fig. 2(a), decreasing $`n_c`$ reduces $`\overline{𝒟}`$ and increases the clustering coefficient $`C`$, where $`C`$ is defined as the probability that two nodes are connected, given that they are connected to a third common node Watts (1999). However, the stronger SW property has a counter-intuitive effect on synchronizability: as the network becomes more heterogeneous (and its SW property stronger), the ratio $`\lambda _N/\lambda _2`$ increases sharply (up to a factor of about 50), indicating that it becomes more difficult to synchronize the oscillators. It is interesting to compare this with the result of Barahona and Pecora (2002), in which the network was never observed to become less synchronizable when $`\overline{𝒟}`$ decreases footnote5 .

In fact, for $`N1`$ and $`m=N2\kappa 1`$, one can show footnote6  that $`\overline{𝒟}`$ scales as $`\mathrm{ln}N`$ for $`n_c=N`$, while it is almost a constant for $`n_c=1`$. However, this is achieved by having almost half of the total load on a single giant center node, and, as a result, $`\lambda _N/\lambda _2`$ scales as $`N`$ for $`n_c=1`$, in sharp contrast to bounded behavior for $`n_c=N`$, making it difficult for the oscillators to synchronize.

The generality of the correlation between homogeneity and synchronizability can be seen in the rigorous bounds on the ratio $`\lambda _N/\lambda _2`$ that are valid for *any connected* network footnote6 :

$$\left(1\frac{1}{N}\right)\frac{k_{\mathrm{max}}}{k_{\mathrm{min}}}\frac{\lambda _N}{\lambda _2}(N1)k_{\mathrm{max}}\mathrm{}_{\mathrm{max}}^e𝒟_{\mathrm{max}}\overline{𝒟},$$
(2)
where $`k_{\mathrm{min}}`$ and $`k_{\mathrm{max}}`$ are the minimum and maximum connectivity, $`𝒟_{\mathrm{max}}`$ is the maximum length of the shortest path between two nodes and $`\mathrm{}_{\mathrm{max}}^e`$ is the maximum of the normalized loads on links (edges), similarly defined as for those on nodes. While $`\lambda _N/\lambda _2`$ is larger for heterogeneous networks because $`k_{\mathrm{max}}/k_{\mathrm{min}}`$ is larger, homogeneous networks have smaller ratios $`\lambda _N/\lambda _2`$ because $`k_{\mathrm{max}}`$ and $`\mathrm{}_{\mathrm{max}}^e`$ are smaller. However, it is the combination of network distances being small (so that $`𝒟_{\mathrm{max}}\overline{𝒟}`$ is small) and the distributions of connectivities and loads being homogeneous that makes the network more synchronizable.

It has been suggested often in neuroscience that smaller $`\overline{𝒟}`$ in the nervous system should imply more efficient communication between neurons, and therefore $`\overline{𝒟}`$ might have been minimized (along with other fitness functions) by evolution Stephan et al. (2000); Karbowski (2001); Latora and Marchiori (2001). The results here, however, suggest that natural selection might have increased homogeneity in the distribution of connections (and hence in the load distribution) at the expense of having larger $`\overline{𝒟}`$, in order to enhance synchronizability. For example, the neurons in some layers of cortical columns are relatively homogeneous in soma size, which has been visualized using the Nissl stain technique kandel . Since the cell body size of a neuron is correlated with the size of its dendritic and axonal expansions braitenberg , which in turn correlate with the number of synaptic connections it receives and sends, it appears that neuronal subnetworks within cortical layers have relatively homogeneous connectivity. The ability of such subnetworks to synchronize facilitates their ability to exhibit bursting behaviors that underlie neural computations.

In conclusion, we have shown that the common belief that smaller networks have better synchronization property can be misleading for a wide class of networks, such as SF and SW networks. Using a general framework for synchronization stability of oscillator networks with arbitrary interaction topology, we have established that SF and SW networks will have reduced ability to synchronize as the heterogeneity of their connectivity distribution increases, even though the average network distance between the oscillators becomes smaller. Increased concentration of load on center nodes, or hubs, appears to be responsible for this behavior. Our results suggest that in order for oscillators in a network to communicate better, and hence to synchronize more effectively izhikevich , a balance between having small communication distance and uniform load distribution is essential.

###### Acknowledgements.

TN is supported by DARPA/ONR grant N00014-01-1-0943. AEM and YCL are supported by AFOSR under Grant No. F49620-98-1-0400. FCH is supported in part by NSF grant DMS-0109001.


text:
# The detection of stellar velocity dispersion drops in the central regions of five isolated Seyfert spirals. Based on observations made with the WHT operated on the island of La Palma by ING in the Spanish Observatorio del Roque de Los Muchachos of the Instituto de Astrofísica de Canarias

## 1 Introduction

The existence of non-axisymmetric components of the galactic potential has been frequently invoked as an efficient way to transport gas from the galaxy scale down to the nucleus to fuel the active galactic nucleus (AGN). In particular, the shocks and gravitational torques induced by a galactic bar provide the non axisymmetric potential invoked in theoretical works (Simkin et al. 1980, Shlosman et al. 1989, Barnes & Hernquist 1991) to make the gas lose angular momentum and therefore facilitate the fueling mechanism. Whether there is no preference for Seyfert nuclei to occur in barred galaxies (Heckman 1980, Simkim et al. 1980, Moles et al. 1995, McLeod & Rieke 1995, Mulchaey & Regan 1997) or whether an excess of bars among Seyferts exists, is still an open issue. Nevertheless, in addition to the large scale bar, another mechanism is needed to drive the gas to the very central regions; one such mechanism is that of nested bars (Shlosman et al. 1989; Friedli & Martinet 1993; Combes 1994; Heller & Shlosman 1994), which has been suggested to fuel molecular gas into an intense central starburst in NGC 2782 (Jogee et al. 1999) and in a Seyfert 2 galaxy, Circinus (Maiolino et al. 2000). However, the analysis of the HST images of Seyferts (Regan & Mulchaey 1999, Martini & Pogge 1999, Pérez et al. 2000, Martini et al. 2003a, 2003b) has shown that such nuclear bars are not ubiquitous, and nuclear spirals are proposed to be one of the channels to feed gas into the central engine.

Stellar kinematical profiles of four double barred galaxy hosts of AGNs obtained with the VLT and ISAAC have recently shown the existence of decoupled kinematics in the nuclear region, as well as a drop of the velocity dispersion in three of them (Emsellem et al. 2001). Although their models could not properly reproduce this drop, Emsellem et al. suggested that it was due to a transient cold nuclear disk fueled by gas inflow along the bar, which has recently formed new stars. Models by Wozniak & Michel-Dansac (2003) have recently succeeded in accounting for this drop in velocity dispersion, which indeed seems to be due to a young stellar population born from dynamically cold gas accreted in a circumnuclear disk. The presence of such disks can be directly related to the circumnuclear structures required to drive the gas to the very central regions, as nuclear (nested) bars or nuclear spirals (see references above).

In most of these studies the environmental effects are in general not considered when studying the properties of Seyfert galaxies and/or comparing them with those of non-active ones (but see Moles et al. 1995). However, the properties of spiral galaxies can be modified even in mild interactions, as shown by Márquez et al. (2002) for a sample of 111 galaxies ranging from truly isolated to mildly interacting spirals. Therefore, we decided to study the properties of isolated Seyfert spirals, and to compare them with a control sample of isolated non-active spirals; the main result is that the active and non-active samples are equivalent in the studied properties: large scale disks, bulges and bars have the same properties, as derived by the analysis of their NIR images (Márquez et al. 1999, 2000). To further delimit eventual differences/similarities and to put some constraints on the possible mechanisms taking place in the circumnuclear regions, kinematical information is needed; therefore, long slit spectra along several position angles were obtained for a number of our sample galaxies, to get both gas (2/3 of the whole sample) and stellar kinematics (for a subsample). The results of the overall kinematical analysis together with the information on metallicities provided by the analysis of the HII regions will be published separately (Márquez et al. 2003). High spectral and spatial resolution spectra in the CaT region have been obtained for a total of eight Seyfert galaxies. In four of them we have found drops in the stellar velocity dispersion of the central regions, together with a hint of such a drop in NGC 6814, comparable to those reported by Emsellem et al. (2001).

Sect. 2 is devoted to the presentation of the observations, data reduction and analysis. The results of the analysis are reported in Sect. 3, and the discussion and conclusions in Sect. 4.

## 2 Observations, Data Reduction and Analysis

The DEGAS sample includes 33 isolated galaxies, 17 with an AGN and 16 without an AGN taken as a control sample. The four galaxies reported here all have an AGN and their infrared properties have already been described (Márquez et al. 1999). As a complement to these infrared images, we have retrieved the HST images for these four galaxies and present them here. All were taken with the F606W filter.

The spectroscopic data for the gas kinematics were collected with the ISIS Spectrograph at the 4m William Herschel Telescope (WHT) in La Palma (Spain). The setup and main characteristics of the observations are given in Table 1. The spectral and spatial samplings were 0.39Å/pixel and 0.36 arcsec/pixel, respectively, and the slit width 1.03 arcsec. The spectral resolution given in Table 1 has been measured as the FWHM of a gaussian fit to several unblended sky lines.

Standard IRAF111IRAF is the Image Analysis and Reduction Facility made available to the astronomical community by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy (AURA), Inc., under contract with the U.S. National Science Foundation. procedures were used for the reduction of the spectroscopic data, following the standard steps of bias subtraction, flat field correction, wavelength calibration with a CuNe lamp observed before and after the target, atmospheric extinction correction, and flux calibration using spectroscopic standards observed through an 8 arcsec wide slit. The sky background level was determined taking median averages over two strips on either side of the galaxy signal. We used cross-correlation techniques as described in Pérez et al. (2000) and references therein, to extract the kinematic information, using as templates for the cross-correlation the observed spectra of several stars. The errors refer to the dispersion in the velocity shifts measured with respect to the different template stars. For clarity, error bars in the stellar velocity distributions are not plotted, but they typically amount to $`\pm `$50 km s-1, with median values of $`\pm `$35 km s-1 for NGC 6951, and up to $`\pm `$70 km s-1 for NGC 2639. The systemic velocity has been chosen as that of the section with the maximum continuum level, except for IC 184, for which the kinematical center (the point that allows the best symmetrization of the rotation curve) is shifted by about 0.8 arcsec to the South with respect to the continuum maximum.

The Calcium triplet equivalent width has been measured following the method described by Terlevich et al. (1990, hereafter, TDT90), which uses two specific continuum bands on either side of the CaT region, after Doppler shift corrections for each section. Due to a shorter spectral range coverage than TDT90, slight modifications of the central wavelength of the “blue” continuum were adopted for all the galaxies except NGC 2639. For the “red” continuum, we had to use much shorter wavelengths than TDT90. The resulting equivalent widths for the inner regions of each galaxy (see Sect. 3) are given in Table 3. In the last two columns the central rest frame wavelength of each continuum band is given (the width used is 15 Å). The differences in the continuum bands used may be the reason why we measure systematically lower equivalent widths (EW) than those reported in TDT90 for the same types of active galaxies. As TDT90 pointed out, to obtain conclusions on the stellar population based on CaII triplet measurements, exactly the same method has to be used. For our data the results are specially tricky due to the fact that the TiO and VO bands contaminate the “red” continuum. In order to quantify the uncertainty due to the choice of the continuum level, we have used three different “red” continua and measured the corresponding EW; the error bars in Fig. 6 are calculated as the dispersion of these three values, and therefore represent lower limits to the true error bars. The minimum uncertainty obtained for the best cases with a flat continuum is 0.5 Å but in most of the galaxies the uncertainty can be larger than 1.5 Å.

CaT equivalent widths have also been measured using the method by Cenarro et al. (2001). The comparison between this method and that of TDT90 has shown that, in the case of dilution due to an important underlying old stellar population, TDT90 seem to be less sensitive to the contamination by TiO and VO lines. Therefore, since our sample galaxies are all early type spirals, we have adopted the TDT90 method.

For our five sample galaxies and for nine of the eleven galaxies with velocity dispersion drops reported in the literature (see Sect. 4) the available archival WFPC2 HST222Based on observations collected with the NASA/ESA HST, obtained at STScI which is operated by AURA Inc., under NASA contract NAS5-26555 optical images (see Table 2 for details) have been used to identify features which, according to numerical models, could possibly be related to the kinematics. To trace the presence of morphological features in the innermost few arcseconds, both raw and sharp–divided images333The result of dividing the original images by median–filtered ones. have been used. Sharp–divided images are very useful to trace asymmetries in the light distribution, such as bars, spiral arms, dust lanes and rings. This technique allows the subtraction of the diffuse background in a very convenient way to look for subtle, small-scale variations and discuss the possible presence of both dust extinguished and more luminous regions (Sofue et al. 1994; Márquez & Moles 1996; Erwin & Sparke 1999; Laine et al. 1999; Márquez et al. 1999; Erwin & Sparke 2002).

## 3 Results

The kinematical data and the images used in the analysis of the data are displayed in Figs. 1 to 4. The points corresponding to gas kinematics are given for comparison (for further details on the gas kinematics see Márquez et al. 2003). The slit positions (measured from North to East) are drawn on the figures, together with the bar sizes and ellipticities as determined in Márquez et al. (1999), for the primary bar in single-barred galaxies (NGC 6814 and NGC 6951) and for primary and secondary (inner) bars in double-barred ones (IC 184, NGC 2639 and UGC 3223). The results for the CaII triplet EW measurements of each galaxy are presented in Table 3. Since we are interested in searching for the properties of stellar populations in the regions where the dips have been detected, the values given in Table 3 refer to the whole region occupied by the drops.

### 3.1 IC 184

The rotation curves along PA=7° (major axis) (Fig. 1) and 97° show that the rotation of the gas and stars are strongly decoupled along PA=7° within $`\pm `$4 arcsec of the center, with a central region of stellar counter-rotation within 1 arcsec, which is also visible along PA=97° (not shown here). In order to know if such kinematics are due to the presence of a counter rotating core, as it is the case for NGC 4621 (Wernli et al. 2002) integral field spectroscopy is required. A small dip of about 20 km s-1 with respect to the maximum value is observed in the velocity dispersion in the innermost $`\pm `$2 arcsec. The stellar velocity distribution along PA=97° shows structure within the innermost $`\pm `$9 arcsec, with a wavy central shape similar to that observed along PA=7°. This extension corresponds to that of the primary bar, which extends along PA=170° (Márquez et al. 2000), i.e. only 17° from the disk major axis. We note that the secondary bar extends $`\pm `$4 arcsec along PA=30° (Márquez et al. 2000) and therefore the velocity difference between the stellar and gas components in this region is most likely due to the non-circular motions caused by the secondary bar. EW(CaT) shows a slight tendency for central values to be higher along both PAs, particularly along PA=97°(Fig. 6a).

The HST image of IC 184 (Fig. 1) shows a structure in the innermost 1 arcsec elongated along PA $``$ 50°, very close to the position angle of the inner bar in the same image (PA $``$ 40°); although this feature is detected in a much smaller region, its PA is consistent with that of the small inner bar already detected in the infrared (Márquez et al. 1999). A structure resembling two straight dust lanes may be seen in the sharp-divided image running along about PA $``$ 45°, below the center to the SW and above the center to the NE, respectively.

### 3.2 UGC 3223

The rotation curve of UGC 3223 along PA=80° (very close to the bar major axis, PA=75°), although rather noisy, shows that the overall rotation of the gas follows, within the error bars, that of the stars (Fig. 2). Positive stellar velocities are systematically lower (down to 80 km s-1) than their negative counterparts for radii between 2 and 20 arcseconds. A small dip ($``$ 10 km s-1) may be observed in the velocity dispersion in the innermost $`\pm `$ 1 arcsec, but needs to be confirmed with data of higher signal to noise. In spite of the low S/N ratio of the spectra, an increase in EW(CaT) is detected in the innermost $`\pm `$1 arcsec, corresponding to the extension of the detected dip (see Fig. 6b).

The HST image of UGC 3223 (Fig. 2) shows an asymmetric structure in the innermost 1 arcsec elongated mainly along PA=80°, and more extended to the West, consistent with that of the small inner bar already detected in the infrared (Márquez et al. 1999). Another extension can be seen to the NE along PA$``$55°, resembling what could be an incomplete ring at r=0.9 arcsec.

### 3.3 NGC 2639

NGC 2639 is a galaxy optically classified as unbarred, but barred at NIR wavelengths (Márquez et al. 1999). The stellar and gas rotation curves of NGC 2639 along PA=135° (disk and main bar major axis, not shown here) show that both components have different velocities (up to 100 km s-1) in the region inside the bar (inner $`\pm `$8 arcsec) (Márquez et al. 2003). The gaseous velocity distribution along PA=0° shows departures of 50 km s-1from pure rotation in the region within r$`<`$4 arcsec (Márquez et al. 2003). Along PA=45° the velocity distribution is quite flat, as expected along a position angle close to the minor axis, whereas clear non-circular motions are traced by the gas kinematics in the central $`23`$ arcsec radius (Fig. 3). A velocity dispersion drop of $``$ 20 km s-1 is detected along the two PAs in the innermost $`\pm `$2 arcsec, being more clearly detected along PA=45° (Fig. 3). A hint of an increase in EW(CaT) is detected along PA=45°and 135°(Fig. 6c).

The HST image of NGC 2639 (Fig. 3) shows a structure elongated along PA=135° in the innermost $`\pm `$1 arcsec (consistent with the thin inner bar detected in the infrared, Márquez et al. 1999), which is somewhat more extended and curved towards the NW.

### 3.4 NGC 6814

NGC 6814 has a small inclination and it is therefore difficult to derive its kinematical parameters accurately. It is a single barred galaxy, with a $`\pm `$12 arcsec bar extending along PA=25° (Márquez et al. 1999). Note that the innermost $`\pm `$3 arcsec show velocity dispersions about 30 km s-1lower than the average values out to about 8 arcsec along PA=120°, an effect which also seems to be detected along PA=30°. Higher S/N data are needed in order to confirm the presence of this drop. Although noiser, this profile resembles that of NGC 6951, both showing a broader central drop than the other three galaxies in our sample. As shown in Fig. 6d, a slight enhancement in EW(CaT) has been detected along the extension of the dip along PA=30°.

The HST sharp-divided image of NGC 6814 shows a small (r $`<`$ 2 arcsec) several-armed spiral structure (in dark grey), with a clear elongation to the West along PA=120° in the innermost 1.5 arcsec and a dust lane along PA$``$ 60° within 1 arcsec (Fig. 4).

### 3.5 NGC 6951

The stellar velocity distributions of NGC 6951 are given in Pérez et al. (2000). We recall that two stellar components have been clearly detected, with amplitude differences up to 50 km s-1along PA=84° and PA=138°. A similar behavior is observed along PA=48° in the region between r=5 and 10 arcsec, which corresponds to stellar velocity dispersion peaks (Fig. 5). Higher S/N data would be required to obtain FWHM measurements for the two separate components, but their physical connection with the FWHM peaks seems to be hinted by our data. The enhancement of EW(CaT) in the central region was already reported and discussed in Pérez et al. (2000).

## 4 Discussion and conclusions

New stellar kinematics for a sample of four active galaxies using the Calcium triplet have been obtained and we have shown that reliable stellar velocity curves can be traced. Three of these galaxies show double bar-like confirmed structure at optical and/or near IR frequencies. The previously reported data for NGC 6951 (a single-barred galaxy with a nuclear spiral, Pérez et al. 2000) have been added to the present study, in particular the velocity dispersion curve.

We have observed a dip in the stellar velocity dispersion in the innermost $`\pm `$(1-3) arcsec of these five active galaxies based on observations made with 4 meter class telescopes. It is worth noticing that the velocity dispersion drops have been detected in three double-barred galaxies and in two single-barred ones, already indicating that the presence of nested bars is not a necessary condition; the drops in the later cases seem to be broader than in the former.

Note that such a dip is until now a rarely observed phenomenon, reported so far only for eleven nearby galaxies (most of them listed in Wozniak et al. 2003b; see references below), three of them from recent VLT observations using the 12CO bands in the near IR (Emsellem et al. 2001). This is most probably due to the difficulty of estimating the FWHM of the stellar component with a sufficient signal to noise ratio and spatial resolution.

Michel-Dansac & Wozniak (2003) have recently modeled the evolution of isolated galaxies over several Gyrs using a self-consistent N-body code including stars, gas and star formation. With such models, they succeeded in accounting for the drop in velocity dispersion observed by Emsellem et al. (2001) in three double-barred and one single-barred active galaxies, with a young stellar population born from dynamically cold gas which is accreted in a circumnuclear disk and creates the dip (Wozniak et al. 2003a, Wozniak & Michel-Dansac 2003). We note that their simulations do not require the presence of embedded bar structures, since such processes can be induced by the primary bar itself (see also Wozniak et al. 2003b).

Since under these circumstances a luminosity dominating young stellar component is expected in the region where the drops are observed, it is worth trying to detect such a component in order to better validate the model. We have attempted to do so by analyzing the equivalent widths of the CaT lines along the dips and by inspecting the inner structure in high resolution HST images. The analysis of HST images for the target galaxies shows that structure is indeed present in the regions where the dips occur in IC 184, UGC 3223 and NGC 2639444With our small sample (5 galaxies) we find no dependence of the drop sizes with respect to those of the outer/inner bars; a larger sample should be needed in order to analyze this point.. The same conclusion is valid when inspecting the sharp-divided frames of those galaxies with available HST archive images and with velocity dispersion drops reported in the literature, namely NGC 1365 and NGC 1808 (Emsellem et al. 2001; for the Sy1 NGC 1097 there are only F218W images, that do not allow to detect such structures), NGC 6503 (Bottema 1993), NGC 3627, NGC 4303, NGC 4579 (Héraudeau & Simien 1998), NGC 4725 (Héraudeau et al. 1999; diffraction spikes impede a reliable determination of the hinted central elongated structure), NGC 4477 (Jarvis et al. 1988), NGC 3623 (de Zeeuw et al. 2002) and NGC 3412 (Aguerri et al. 2003). Note that the images are at a wavelength range that does not include the CaT, so we do not expect to find exactly the same structure. Also the spectroscopic structure is diluted because of the poorer ground-based spatial resolution. Therefore, the presence of structure in HST images at the scales of the velocity dispersion drops is only a consistency check. The existence of a direct correlation between the two observables needs further investigation with similar spatial resolution data.

The velocity dispersion drop in NGC 6814 seems to be as broad as that in NGC 6951, and it is only for NGC 6814 and NGC 6951 (LINER/Sy2) that the drops seem to be present with no associated elongation in the images. This may indicate that such structures either are seen face-on in these cases or are absent. For the other galaxies, if the elongations are interpreted as intrinsically circular disks, these disks would be placed almost perpendicular to the main disk in NGC 2639 (i$`_{innerdisk}`$ $``$ 74°), at intermediate angle for IC 184 (i$`_{innerdisk}`$ $``$ 40°) and coplanar with the main disk in UGC 3223.

The Calcium triplet equivalent width, EW(CaT), is expected to increase for younger populations. Unfortunately, the comparison with stellar population models to test the age of such a population is rather risky, due to the already mentioned problem with the determination of the continuum bands which impedes a quantitative comparison with other authors (for instance, the models developed by García-Vargas et al. 1998 require EWs higher than those we have measured to invoke an important contribution of super-giant stars). The variations on the values of the EWs can still be compared but, due to the accuracy achievable by using EW(CaT) measurements on our data, the analysis is not straightforward since we expect the relative differences in the EW(CaT) to be comparable to our accuracy (from 1 to 1.5 Å), and consequently such variations can, strictly speaking, correspond to a constant value. Nevertheless, a hint of a local enhancement in the drop region is suggested by the data, a behavior that could be due to the presence in these regions of a population different from that of the bulge.

To confirm or discard Wozniak et al.’s models the analysis of higher S/N spectra in a much broader wavelength range is needed to constrain the age of the observed stellar populations detected in the central dips. On this respect, we notice that some of the galaxies show evidence for the presence of an intermediate (NGC 3627 and NGC 6503, González Delgado et al., in preparation) or young (NGC 4303, Colina et al. 2002) population in the regions where the dispersion drops take place and where elongations are seen in the HST sharp-divided images.

The case of NGC 6951 is specially interesting since an alternative explanation for the presence of the apparent drop within r=5 arcsec in this galaxy would be the existence of two stellar kinematic components. This appears as a possibility to be checked in the other four galaxies with higher spectral resolution and S/N data. It would therefore be of particular interest to re-observe these five galaxies with an 8 meter class telescope, to define with better precision the characteristics of the dip, in particular its disk nature (for which several position angles are required), and also to observe a larger sample of galaxies with and without an AGN. Note that among the 16 galaxies with reported drops – 5 of them presented for the first time in this work – 15 host AGNs555Their nuclear classification has been taken from NED and is given in Table 2., either Seyfert 1 (3), Seyfert 2 (7), LINERs (2), or intermediate type LINER/Sy2 (3); the only galaxy with no evident signs of nuclear activity is NGC 3412 (Ho et al. 1995). The DEGAS sample is particularly well adapted for this purpose, since it is well defined and comprises only isolated galaxies, where the complicating effects of galaxy interactions are not present.

###### Acknowledgements.

We thank the referee, Eric Emsellem, for his excellent report, which has largely helped us improve the clarity of the paper. I.M. acknowledges financial support from the IAA and the Spanish Ministerio de Ciencia y Tecnología through a Ramón y Cajal fellowship. This work is financed by DGICyT grants PB93-0139, PB96-0921, PB98-0521, PB98-0684, ESP98-1351, AYA2001-2089, AYA2001-3939-C03-01 and the Junta de Andalucía grant TIC-144. We acknowledge financial support from the Picasso Programme d’Actions Integrées of the French and Spanish Ministries of Foreign Affairs. F.D. acknowledges financial support from CNRS-INSU for several observing trips. E.P. acknowledges financial support from the CONACYT (Máxico) grant 36132-E. J. Maza gratefully acknowledges support from the Chilean Centro de Astrofísica FONDAP 15010003. We acknowledge Ignacio Marrero for his collaboration in some data acquisition and calibration. This research has made use of the NASA/IPAC extragalactic database (NED), which is operated by the Jet Propulsion Laboratory under contract with the National Aeronautics and Space Administration.