text:
# Global accretion disk simulations of magneto-rotational instability

## 1 Introduction

Accretion processes in astrophysical disks lead to enormous luminosity and huge changes in disk structure during astronomically short times. Efficient transport mechanisms are necessary to achieve such short time-scales. Anisotropic turbulence appears to be a major physical condition to provide astrophysical disks with strong transport. As these disks generally exhibit increasing specific angular momentum towards larger radii and thus fulfill the stability criterion of Rayleigh, they do not give rise to an instability leading to turbulence by themselves. Searches for instabilities in disks with rotation profiles similar to a Keplerian one unveiled several ways to turbulence being more or less favorable with respect to their prerequisites for the disk configuration. Gravitational instability needs the disk to be either cool or massive. Nonlinear and nonaxisymmetric perturbations require a severe additional perturber near the disk; conditions for instability in a purely hydrodynamical disk were derived by Dubrulle (1993). Hydrodynamic instability essentially comes down to violating the Rayleigh criterion saying that a rotation profile with an increasing specific angular momentum with radius is hydrodynamically stable, that is for $`l^2/r>0`$. Since the typical length of perturbations caused be external forces is supposed to be very large, the required amplitude of the perturbations has to be considerable, too, in order to violate the Rayleigh criterion locally. If the length scale of the perturbation is comparable to the radius, a strong alteration of the Keplerian velocity of several per cent is needed. Convection was shown to deliver either negligible transport (Stone & Balbus 1996) or inward angular momentum transport (Kley et al. 1993).

The requirements for the magnetic shear-flow instability (Balbus & Hawley 1991) do match astrophysical conditions in accretion disks in many configurations. All it needs is a radially decreasing angular velocity and a weak magnetic field threading the rotating object. It can even be shown that the temperature range applicable to the magnetic shear-flow concept is very broad; even very small ionization fractions are sufficient to magnetize a disk in many cases (Balbus & Hawley 1998).

First numerical approaches to the magnetic shear-flow instability tackled the local problem; the linear analysis as described by Balbus & Hawley (1991) were immediately followed by 2D simulations (Hawley & Balbus 1991) of a small box-shaped domain which was cut out of the disk. These computations confirmed the relation between magnetic-field strength and wavenumber derived from the linear analysis earlier, they showed that the maximum growth rate of any wavenumber is independent from the field strength and that the system is capable of transporting angular momentum. Because of being restricted to axisymmetric problems, they could not provide self-sustained turbulence which needs dynamo-action, and the slow decay of the turbulence is an indirect effect of the Cowling theorem. Improved computations dealt with a three-dimensional even though local configuration, and particular care was taken for the radial boundary conditions which are not simply periodic, but account for the shear due to differential rotation (Hawley et al. 1995). These comprehensive computations indeed resulted in magnetically sustained turbulence whose anisotropy causes efficient outward angular momentum transport. This work was followed up by stratified models (Stone et al. 1996) covering more than 50 orbital periods of the box cut-out. As the computational domain covered 2 density scale heights, this work was a first step towards global simulations, followed up by similar approaches such as Ziegler & Rüdiger (2000a, b).

Linear studies of global configurations of disks threaded by magnetic fields in various directions were carried out. Curry & Pudritz (1995) investigated the stability for vertical and azimuthal fields threading the disk. They found in detail that the actual initial field geometry does not strongly depend on field topology as was suggested by the numerical simulations of Hawley & Balbus (1991). Rüdiger et al. (1999) particularly addressed the angular momentum transport in their linear study. These investigations are now followed by nonlinear simulations of a compressible fluid with density stratification in global computational domains. The most recent global approach was presented by Hawley (2000) following the evolution of a thick torus under the influence of an external magnetic field threading parts of the computational domain. The magnetic shear-flow instability was found to set in quickly causing enough turbulent viscosity to soon form a Keplerian velocity profile. All the above mentioned numerical studies integrate the ideal MHD equations omitting magnetic diffusivity. Local simulations including diffusivity were performed by Fleming et al. (2000) which proved the onset of instability even for low conductivity.

The full understanding of accretion disks implies self-excited dynamo action as well as angular momentum transport. Can positive angular momentum transport be coupled with a suitable kinetic helicity providing the expected dynamo action according to the $`\alpha `$-effect principle? As the kinetic helicity in stratified, rotating disks is expected to be negative, a wrong sign (positive) would follow for the dynamo-$`\alpha `$-effect according to the conventional $`\alpha `$-theories. We present global three-dimensional diffusive simulations and study the angular momentum transport as well as dynamo action and the sign of $`\alpha ^{\mathrm{dyn}}`$ as a consequence of the correlation with the flow.

The present computations are not meant as particular simulations of a star formation process. They focus on the applicability of the magnetic shear-flow instability for an astronomically fast process like the formation of a star. Such global simulations are still a challenge for modern computers and fast algorithms.

## 2 The simulation setup

The computations presented here make use of the ZEUS-3D code developed for astrophysical problems of magnetohydrodynamics (see the key papers Stone & Norman 1992a,b; Stone et al. 1992 for numerical, Clarke et al. 1994 for technical details). We use cylindrical coordinates and an extensive computational domain covering radii $`r=5`$ to 6 and $`r=4`$ to 6, a vertical extension of $`h=1`$ to +1, and the full azimuthal range of $`\varphi =0`$ to $`2\pi `$. See Table 1 for a list of models presented in this Paper. In this approach, we assume an isothermal disk to save computation time on the energy equation. The remaining system for integration is

$`{\displaystyle \frac{\rho }{t}}+div(\rho u)`$ $`=`$ $`0`$ (1)
$`{\displaystyle \frac{\rho u}{t}}+div(\rho uu)`$ $`=`$ $`gradp\rho grad\mathrm{\Phi }+J\times B`$ (2)
$`{\displaystyle \frac{B}{t}}`$ $`=`$ $`curl(u\times B)+\eta \mathrm{}B,`$ (3)
where $`\rho `$, $`u`$, and $`B`$ are the density, velocity, and magnetic field resp.; $`p`$ is the pressure ($`p\rho `$ in our model), $`\mathrm{\Phi }`$ is the gravitational potential (solely from a central mass $`M`$), $`J`$ is the current density, and $`\eta `$ is the magnetic diffusivity which is not an original ingredient to the ZEUS code. It is constant in time and space. The computations of the electromotive force in the routines emfs, mocemfs, and hsmoc of ZEUS-3D are extended with the appropriate $`\eta curlB`$ components. We chose $`\eta =0.001`$ to 0.01. The additional time-step criterion resulting from the diffusivity is roughly 0.1 for the finest grid used here. It is therefore irrelevant for the determination of the time-step which is typically $`10^4`$ or, in the case of strong fluctuations of velocity and magnetic field, one or two orders of magnitude lower. In physical units, this diffusivity is large and in fact accounting for a subgrid turbulent diffusivity which cannot be resolved.

The gravitational potential is spherically symmetric, whence the $`z`$-component of the gravitation is retained within the disk. We therefore obtain a density stratification unlike the computations by Armitage (1998) who omits the $`z`$-component of the gravitation and applies periodic boundary conditions for the upper and lower boundaries. The sound speed is always $`c_{\mathrm{ac}}=10`$ which is roughly $`0.07u_\mathrm{K}`$ in terms of the average Keplerian velocity in the simulated ring. The initial density distribution is Gaussian with $`r`$-dependent density scale-height, $`\rho (r,z)=\rho _\mathrm{c}\mathrm{exp}(z^2/H(r)^2)`$.

The magnetic field threads the disk vertically with a mere $`z`$-component. We tried a homogeneous initial field $`B_z(r)=\mathrm{const}`$ for the $`z`$-boundaries as well as a field

$$B_z(r,t=0)=b_0r^1\mathrm{sin}[2\pi (rr_\mathrm{i})/(r_\mathrm{o}r_\mathrm{i})]$$
(4)
– where $`r_\mathrm{i}`$ and $`r_\mathrm{o}`$ are the inner and outer boundary radii resp. – which has zero total magnetic flux through the $`z=\mathrm{const}`$ surfaces. The choice of the initial field will have implications for the topology of the field later on. As the magnetic flux through the surfaces is kept constant, the first approach will always preserve a non-zero average magnetic field through the vertical boundaries, whereas the field penetrating the top and the bottom of the computational domain can completely vanish in the second choice of initial fields. The results presented here were obtained with the second approach of vanishing average initial field $`B_z`$. The parameter $`b_0`$ was chosen between 50 and 100 giving an amplitude of the initial magnetic field of $`\mathrm{max}(B_z)=10`$ to 20 corresponding to maximum Alfvén speeds at the upper and lower boundaries of $`u_\mathrm{A}=13`$ to 26 or $`u_\mathrm{A}=1.3c_{\mathrm{ac}}`$ to $`2.6c_{\mathrm{ac}}`$. The Alfvén velocity in the equatorial plane is subthermal with $`u_\mathrm{A}=0.32`$ to 0.63 or $`0.032c_{\mathrm{ac}}`$ to $`0.063c_{\mathrm{ac}}`$.

The initial velocity field is a merely Keplerian motion following $`u_\varphi =\sqrt{GM/r}`$, where $`G`$ is the gravitational constant and $`M`$ is the central mass which is $`10^5`$ in our computational units. The total disk mass is 35,000 for Model Ia, Ib, and IV and 56,000 for the other models. Times are henceforth measured in orbital periods which convert by $`T_{\mathrm{orb}}=0.159`$ from the arbitrary units of the code.

The number of cells in each coordinate direction was $`31\times 61\times 351`$ for the $`z`$-, $`r`$-, and $`\varphi `$-directions. Tests with up to 621 cells in $`\varphi `$-direction have been carried out, but the increased computation time does not allow reasonable periods to be covered by the simulation. The aspect ratio of the grid cells is not unity. Finest resolution is achieved in radial direction, lowest in azimuthal direction: $`\mathrm{\Delta }z/\mathrm{\Delta }r=2`$, $`r\mathrm{\Delta }\varphi /\mathrm{\Delta }r=2.18`$ for the smallest radius and $`r\mathrm{\Delta }\varphi /\mathrm{\Delta }r=3.28`$ for the outer radial edge. The fastest growing mode of the magneto-rotational instability has a wavelength of $`\lambda _{\mathrm{inst}}=2\pi \sqrt{16/15}u_\mathrm{A}/\mathrm{\Omega }`$. With the angular velocity of $`\mathrm{\Omega }=28.2`$ in the middle of the computation domain, we obtain wavelengths between 0.4 and 0.9 with the above given Alfvén velocities of the initial field strength. These wavelengths are upper limits as they refer to the sites of maximum field and minimum density which need not coincide necessarily. Yet, they are well resolved by the computation domain covering $`z`$-ranges of 2.0 and two different $`r`$-intervals of 1.0 and 2.0.

The integration of the magnetic fields was done with the evolution of the electromotive forces by the Constraint Transport scheme (CT, Evans & Hawley 1988) which preserves a divergence-free configuration.

Two terms adding a source of viscosity to the hydrodynamic equations extend the above Navier-Stokes equation. The von Neumann-Richtmyer artificial viscosity depends on the square of velocity gradients and acts only on decreasing velocity in the direction of propagation, i.e. compression. The strength of this term is denoted by $`q_{\mathrm{con}}`$ in the implementation of the ZEUS-3D code. The second term depends linearly on the velocity gradient and acts on both compression and expansion; the strength is represented by $`q_{\mathrm{lin}}`$. Tests with a stable Keplerian flow found the choice of $`q_{\mathrm{con}}=0.1`$ and $`q_{\mathrm{lin}}=1.0`$ to be suitable. These values were used throughout all the runs described here. The Courant number determining the ’safety’ of a certain maximum allowed time step derived from the velocities, magnetic fields, and the artificial viscosities, is set to 0.5.

The ZEUS-3D code provides three methods of interpolation in the transport step: the first-order donor-cell method, the second-order van-Leer method, and the third-order piecewise-parabolic-advection method. The original ZEUS-3D copy chooses the second-order interpolation for the transport of all quantities, and we have not altered this switch for the computations given here.

Occasionally, pockets of extremely low density may emerge coupled with extraordinarily high speeds exceeding the Keplerian velocity. We suppress the existence of such pockets by a mass replenishment as soon as the density of a certain cell drops below $`10^4`$ the central density. The actual mass being thus added to the total disk mass is found to be negligible. The mass replenishment thus acts like a lower limit for the time step.

### 2.1 Boundary conditions

The vertical boundary condition on the $`z=\mathrm{const}`$-faces of our computational domain are closed for the flow, i.e. $`v_z=0`$. The boundary condition is stress-free in the sense of $`v_r/z=0`$ and $`v_\varphi /z=0`$. The magnetic field has to fulfill the ‘opposite’ boundary condition as it is only allowed to cross the boundary perpendicularly, which is actually implemented affecting the electromotive force $`=u\times B\eta curlB`$ such that $`_z=0`$, $`_r/z=0`$, and $`_\varphi /z=0`$.

No magnetic field lines are allowed to penetrate the radial boundary surfaces. Computations were carried out with two choices of radial boundary conditions regarding the velocity. The first does not allow flow through the top, bottom and outer radial boundaries. The inner boundary was chosen to be ‘outflow’, that is, the physical arrays are constantly extrapolated into the ghost zones of the computational domain. The only exception is an inflow, in our case a $`u_r>0`$ at the inner boundary. This velocity component is reset to zero then. All radial velocities $`u_r<0`$ are copied to the ghost zones beyond the computational domain. We will henceforth refer to these models as Model Ia (only inner boundary has ‘outflow’) and Ib (all both radial and both vertical boundaries have ‘outflow’).

Since the outflow condition is likely to empty the disk on the viscous time-scale, the second choice tries to account for the accreted matter and feeds the disk at the outer radial boundary. We sum up the mass loss rate at the inner boundary $`\dot{M}=r_{\mathrm{min}}\mathrm{\Delta }\varphi \mathrm{\Delta }z_i_k\rho _{ik}(u_r)_{ik}`$ for the total mass loss. A constant, small inflow velocity $`u_{\mathrm{in}}`$ is assumed for the outer boundary. Accordingly, the density at the outer boundary is determined to account for the mass flux at the inner boundary. Since we use the average influx, its $`z`$-dependence (nor even the $`\varphi `$-dependence) is not used for a $`z`$\- (nor even $`\varphi `$-) dependent inflow from far radii. We tested inflow conditions with homogeneous and Gaussian and density distribution. It was supposed that an average influx of mass will establish according to the current viscosity (either numerical or turbulent) in the disk. These models are refered to as Model II-VIII in the following.

The azimuthal velocity, $`u_\varphi `$, is extrapolated into the radial boundary zones by a power law $`r^{1.5}`$ based on the last zone of the computational domain. The velocities are thus not Keplerian, they just follow the same radial power law. Since we cover the full azimuthal range in $`\varphi `$, the boundaries at the $`\varphi =\mathrm{const}`$ surfaces are periodic, naturally.

Even though the ZEUS-3D code is remarkably stable for a broad variety of configurations and physical problems, a comprehensive study of the influence of the computational schemes provided were found to be inevitable for trustworthy results.

## 3 Patterns, transport, spectra

### 3.1 Glancing at the solutions

Typical vertical slices from a run of Model II are shown in Figures 1 and 2. The central mass is located on the left side, outside the computational domain. The shadings represent the azimuthal components of the magnetic and velocity fields, while the arrows are the $`(r,z)`$-projections of the actual vectors. Shaded areas near the equatorial plane of the disk show the typical Keplerian velocity profile. However, at high altitudes above the equator, a strong super-rotation emerges coupled with relatively high infall velocities. (These high radial velocities have been cut in Figure 1 to better show the velocities in the disk.)

The computational domain thus covers part of the corona of the disk where the density falls below $`10^3`$ times the equatorial density. The infall velocity in the corona is of the same order of magnitude as the Keplerian velocity, whereas the non-orbital velocities within the disk are much smaller. The actual flow in the disk as shown in Figure 1 is limited by $`0.35c_{\mathrm{ac}}`$.

After saturation of the turbulence, the temporal behavior of $`\alpha _{\mathrm{SS}}`$ changes. The fairly smooth function turns into an oscillatory behavior. The typical feature of this moment is the zero-crossing of the Reynolds stress which is essentially negative all through the saturation parts of the simulations. It is expected that the infall of matter at the outer radial boundary tends to perform relaxation motions, in particular for Models II–IV where matter falls in homogeneously. The exclusion of the outermost six cell planes from the computation of the spatially average $`\alpha _{\mathrm{SS}}`$ does not alter the result (the actual relaxation motions are observed only in the outermost three cell planes). Relaxation motions will consist mainly of $`u_z`$ components which do not affect radial angular-momentum transport directly.

A problem connected with the homogeneous feeding of the disk from the outer boundary is the appearance of a ‘density front’. This is the reason why the density at the outer boundary is highest which may look odd in the context of an accretion disk. This incoming density enhancement can be partly leveled by the presumed infall velocity at the outer boundary.

### 3.2 Angular momentum transport

According to the standard model of a turbulent disk according to Shakura & Sunyaev (1973), the $`r\varphi `$-element of the stress tensor of velocity and magnetic-field fluctuations,

$$W_{r\varphi }=\frac{\rho u_r^{}u_\varphi ^{}\frac{1}{4\pi }B_r^{}B_\varphi ^{}}{\rho },$$
(5)
scales like the speed of sound such as

$$W_{r\varphi }=\alpha _{\mathrm{SS}}c_{\mathrm{ac}}^2,$$
(6)
parameterized by an unknown $`\alpha _{\mathrm{SS}}`$ which is smaller than unity in the case of subsonic turbulence. We normalize the stress by the average density in the computational domain; other authors used the density averaged in the equatorial plane leading to lower limits for $`\alpha _{\mathrm{SS}}`$. This stress measures the angular momentum transport in the disk; positive values mean outward transport. We average this quantity over the entire computational domain according to

$$\frac{\rho u_r^{}u_\varphi ^{}r𝑑\varphi 𝑑z𝑑r}{r𝑑\varphi 𝑑z𝑑r}$$
(7)
for the velocity and

$$\frac{1}{4\pi }\frac{B_r^{}B_\varphi ^{}r𝑑\varphi 𝑑z𝑑r}{r𝑑\varphi 𝑑z𝑑r}$$
(8)
for the magnetic field fluctuations. As the azimuthal velocity contains a large supersonic axisymmetric mode – the Keplerian flow – we adopt $`u_\varphi ^{}=u_\varphi \sqrt{GM/r}`$ for the azimuthal velocity fluctuations. The averaging is in fact carried out on a discrete computational mesh by

$$\frac{_jr_j_i_k\rho u_r^{}u_\varphi ^{}}{n_in_k_jr_j}$$
(9)
for the velocity part and

$$\frac{1}{4\pi }\frac{_jr_j_i_kB_r^{}B_\varphi ^{}}{n_in_k_jr_j}$$
(10)
for the magnetic part.

The time series of the outflow model Ia is shown in Fig. 3. Model Ia in particular shows the typical short-lived occurrence of an equivalent to the channel solution as was found in local box simulations; a similar but less pronounced behavior is found in Model Ib. This period is the only time in both Models Ia and Ib when Reynolds stresses dominate the Maxwell stresses.

Time series for long runs of Models II and Model III are shown in Figure 5 and 7. Generally, the accretion model shows stronger angular momentum transport than the the simpler inner-outflow Model Ia. Models Ia and Ib run out of matter after about 8 orbital revolutions with enabled magnetic fields, and the $`\alpha _{\mathrm{SS}}`$ values are no longer meaningful.

A striking feature of the accretion Models II and III is the change of sign in the Reynolds stress after about 4–6 orbits. We find that, in a long run, the outward angular momentum transport is a sole consequence of the magnetic stresses. Additionally, the Reynolds part shows strong amplitude variations unlike the Maxwell stress.

Figures 4 and 6 show the respective temporal evolution of the fluctuative kinetic energies. For the $`z`$\- and $`r`$-components, the zero mode is subtracted such as $`u_z^2=u_z^2u_z^2`$; the azimuthal fluctuative energy is found by $`u_\varphi ^2=(u_\varphi \sqrt{GM/r})^2`$. All three models show a clear energy ‘sorting’ of the $`z`$-, $`r`$-, and $`\varphi `$-components in this order.

A particularly long run was performed with a computational domain only half as large as Models II and III and is denoted by IV. The initial and boundary conditions are like in Model II. The long-term behavior is not too similar to that of Model II with very strong oscillations of Reynolds versus Maxwell stresses. A behavior comparable with Model II may be spotted during the first 8 orbital periods. After that time, the energy sorting of the radial and vertical components is given up. All models from Ia to IV appear to be not fully saturated with respect to the ever-growing kinetic energies.

We should note that a homogeneous infall might be suitable for a realistic disk in the sense that the infalling matter does not ‘know’ the structure of the disk. Nevertheless, numerical problems may occur when the incoming homogeneous ’matter front’ meets the near-Gaussian density structure. Models V and VI apply an inflow with pre-defined Gaussian infall profile over $`z`$. The respective $`\alpha _{\mathrm{SS}}`$-profiles are shown in Figures 8 and 10, the first having an infall velocity of $`u_{\mathrm{in}}=0.01`$, the second $`u_{\mathrm{in}}=0.1`$. The initial behavior is similar to the respective Models II and III including the change of sign of the Reynolds stress after a couple of orbits. After as many as 10 orbital periods, the disk turns into a significantly different regime with enhanced outward transport in both models.

The respective kinetic energies are given in Figures 9 and 11. In contrast to the above studied models, energies do not develop in an ever-growing way, but appear to fluctuate about an average. The only exception is the regime change near $`t=42T_{\mathrm{orb}}`$ (Model V) and $`t=29T_{\mathrm{orb}}`$ (Model VI) when the average (essentially of $`u_z^2`$) changes, but the stationary appearance remains.

The distribution of $`\alpha _{\mathrm{SS}}`$ versus time and the $`z`$-direction is shown in Figure 12 for Models V. Outward angular momentum transport is marked with white areas, inward transport by dark areas. Strong transport is observed at large distances from the equator shortly after the onset of the instability. Regions of outward transport migrate towards the equator. This is particularly marked in the Model-V plot for which we found a development into a strong-transport regime (cf. Figure 8). Positive angular momentum transport is not restricted to high altitudes anymore.

### 3.3 Spectral decomposition

Power spectra resulting from Fourier transforms of the magnetic field of Model V over the the azimuthal and radial directions are shown in Figures 13 and 14, respectively. The azimuthal decomposition of the velocity field is shown in Fig. 15. The individual azimuthal spectra for all points in the ($`z,r`$)-plane were averaged, as were the points of the ($`\varphi ,z`$)-plane for the radial spectra. In the azimuthal spectra, we observe slight maxima in power at $`k=5`$ and $`k=8`$, followed by a strong, roughly power-law decline between $`k20`$ and $`k100`$. The magnetic spectrum is steeper than the kinetic one. The underlying power-laws would be $`Ek^3`$ to $`Ek^4`$ which is significantly steeper than a Kolmogorov spectrum, and is similar to what Armitage (1998) found for a global simulation in the far-$`k`$ range. The Kolmogorov spectrum of isotropic turbulence would exhibit a wavenumber exponent of $`5/3`$; here we face MHD turbulence which implies an effect of Alfvén waves impeding the transfer of energy towards smaller scales. The energy spectrum of isotropic MHD turbulence is $`Ek^{3/2}`$, whence more shallow than the Kolmogorov spectrum. Contrasting with the azimuthal decomposition, the radial spectrum matches the theoretical MHD spectrum over the entire range of modes as can be seen in Fig. 14.

### 3.4 Accretion rate

Choosing the rough dimensions of an inner protostellar system of 100 AU radius for scaling the mass flux through our models, we obtain an accretion rate in physical units as shown in Fig. 16. Widely accepted values for the accretion rates of protostellar disks around solar-type stars are $`10^7`$ to $`10^6M_{}`$/yr. The compilation of results for those stars by Strom (1994) rather favors $`10^6M_{}`$/yr. The accretion rate increases with stella mass and may already exceed $`10^5M_{}`$/yr for a 3-solar-mass star. Our results such as given in Fig. 16 are thus not too far from the accretion rates deduced from the observations.

### 3.5 Comparison with other studies

Among stratified local simulations, we find a number of very similar results as regards the magnitude of $`\alpha _{\mathrm{SS}}`$, such as in Stone et al. (1996) who obtained $`0.01`$, Hawley et al. (1996) where an average of 0.016 was found, with typical values of 0.01. The largest box used delivered 0.02, the smallest, 0.003. Upon studying the dependence of the turbulence on the rotation law $`\mathrm{\Omega }r^q`$, Hawley et al. (1999) found $`10^2`$ for $`q=1.5`$. A most recent local-box simulation by Ziegler & Rüdiger (2000b) covering 288 orbital periods resulted in a space-time average of 0.015, where the Reynolds part was $`2.810^3`$ and the Maxwell part was $`1.210^2`$. The stress was normalized with the averaged equatorial-plane pressure and thus provides a lower limit for $`\alpha _{\mathrm{SS}}`$. The early stratified-box simulations by Brandenburg et al. (1995) provided $`\alpha _{\mathrm{SS}}`$ an order of magnitude lower between 0.001 and 0.005. All the models agree upon the dominance of Maxwell stresses over Reynolds stresses.

It is most interesting to note that other global simulation such as performed by Armitage (1998) or Hawley (2000) result in an order of magnitude higher values of 0.2 to 0.3. Contrasting with this, the global spherical, but unstratified simulation by Drecker (2000) gives again an $`\alpha _{\mathrm{SS}}`$ which is an order of magnitude lower than shown here, i.e. $`10^3`$. Those simulations alone included a physical viscosity. The simulations of Hawley (2000) – and with somewhat different objectives, those of Machida et al. (2000) – had open boundaries in radial and vertical directions. This fact limits the run-time of the models which are prone to running out of matter.

## 4 Dynamo effect

The long-term existence of non-zero toroidal fields in each hemisphere as shown for example in Figure 19 for Model V indicates the possible action of a dynamo in the disk.

The local-box simulations by Hawley et al. (1996) impose severe constraints on the growth of an average magnetic field because of the periodicity of the boundary conditions. A non-zero term of the average magnetic-field growth is due to the shear and requires the presence of an average radial component, in the local nomenclature $`B_x`$. We follow their derivation of the averages by replacing the volume curl integral by a vector-produc surface integral of (3). The average

$$\frac{B}{t}=\frac{1}{\pi \mathrm{\Delta }r^2\mathrm{\Delta }z}\underset{V}{}curl(u\times B\eta curlB)dV$$
(11)
is equivalent to

$$\frac{B}{t}=\frac{1}{\pi \mathrm{\Delta }r^2\mathrm{\Delta }z}\underset{S}{}𝑑S\times (u\times B\eta curlB)$$
(12)
and delivers average components as

$`{\displaystyle \frac{B_z}{t}}`$ $`=`$ $`{\displaystyle \frac{}{r}}u_rB_z`$ (13)
$`{\displaystyle \frac{B_r}{t}}`$ $`=`$ $`{\displaystyle \frac{}{z}}u_rB_z+\eta {\displaystyle \frac{}{z}}{\displaystyle \frac{B_z}{r}}`$ (14)
$`{\displaystyle \frac{B_\varphi }{t}}`$ $`=`$ $`{\displaystyle \frac{}{z}}u_\varphi B_z\eta {\displaystyle \frac{1}{r}}{\displaystyle \frac{}{z}}{\displaystyle \frac{B_z}{\varphi }}{\displaystyle \frac{^2B_\varphi }{z^2}}`$ (15)
$`{\displaystyle \frac{}{r}}u_\varphi B_r+\eta {\displaystyle \frac{}{r}}{\displaystyle \frac{1}{r}}{\displaystyle \frac{rB_\varphi }{r}}`$
These considerations show that slight initial fluctuations in $`B`$ are able to create non-zero average magnetic fields at any time, even if their $`B^{}`$ vanishes, as a consequence of shear and radial accretion through $`u_{\mathrm{in}}`$. In particular, the non-constant vertical profile of the orbital velocity produces a $`B_\varphi `$ from the initial perturbation $`B_z`$. Fluctuations in the accretion flow generate $`B_r`$ and $`B_z`$. The actual radial differential rotation then quickly produces $`B_\varphi `$ from even very small $`B_r`$, and toroidal fields are prone to dominate the disk fields.

### 4.1 Sketching mean-field dynamos

Kinematic dynamos have been extensively examined in numerous previous publications concerning various types of geometries. A widely used mean-field approach to circumnavigate the difficulties with resolving the small scales in a model applies the likely amplification of a magnetic field through a fluctuative electromotive force parallel to the actual large-scale magnetic field. This effect is generally expressed by the term ‘$`\alpha `$-effect’ and can be thoroughly studied in e.g. Krause & Rädler (1981). In a model splitting large and small scales, the induction equation extends beyond the large-scale electromotive force such as

$$\frac{B}{t}=curl\left(u\times B+\right).$$
(16)
The $`\alpha `$-effect is part of the development of the small-scale electromotive force as

$$_i=\alpha _{ij}^{\mathrm{dyn}}B_j\eta _{ijk}B_k_{,j}+\mathrm{},$$
(17)
where $`\alpha _{ij}^{\mathrm{dyn}}`$ and $`\eta _{ijk}`$ are tensors in general. If written in components, the induction equation shows that the essential component of the $`\alpha ^{\mathrm{dyn}}`$ tensor is the $`\alpha _{\varphi \varphi }^{\mathrm{dyn}}`$ component which converts azimuthal fields into radial and vertical fields, whereas the differential rotation converts radial and vertical fields back into the toroidal component by the large-scale part of $`u\times B`$. This approximation is usually referred to as an $`\alpha \mathrm{\Omega }`$-dynamo.

Second-order correlation approximation leads to a relation of $`\alpha ^{\mathrm{dyn}}`$ with the helicity of the flow,

$$\alpha ^{\mathrm{dyn}}=\frac{1}{3}\underset{0}{\overset{\mathrm{}}{}}u^{}(t)curlu^{}(t\tau )𝑑\tau ,$$
(18)
which is, when approximated by a typical time-scale $`\tau _{\mathrm{corr}}`$,

$$\alpha ^{\mathrm{dyn}}u^{}curlu^{}\tau _{\mathrm{corr}}.$$
(19)
The same principle is applicable to the current helicity (Keinigs 1983) giving

$$\alpha ^{\mathrm{dyn}}\frac{\eta }{B^2}B^{}curlB^{}.$$
(20)
Among other issues, the evolution of these quantities for the helicity will be subject of the following sections.

### 4.2 Symmetry of magnetic fields

Special attention is payed to the symmetry of the magnetic fields of the solutions. In general, the preference of a certain symmetry is a tool to check the consistency with mean-field dynamo models. Additionally, dipolar (antisymmetric) fields are supposed to support the launch of winds forming jets from the disk better than quadrupolar (symmetric) fields, with the latter rather providing closed field lines within the disk. The mean-field disk dynamo of Rekowski et al. (2000) produced dipolar magnetic fields for negative $`\alpha ^{\mathrm{dyn}}`$ in the northern hemisphere.

A number of mean-field dynamo models in one to three dimensions delivered varying symmetry of the solutions with a major influence of the boundary conditions – vacuum versus perfectly conducting. Three-dimensional investigations of stability by Meinel et al. (1990) and Elstner et al. (1992) found quadrupolar solutions for the condition of low-conductivity surroundings of the disk, whereas perfectly conducting boundaries delivered dipolar solutions. A local model with periodic boundaries will thus not be able to answer the question of the final symmetry of the solutions. An ideal global model would comprise the entire investigated object and is less depending on the surroundings (which will be vacuum in very good approximation).

The parity of a the toroidal magnetic field is measured by

$$P=\frac{E_\mathrm{s}E_\mathrm{a}}{E},$$
(21)
where $`E_\mathrm{s}`$ and $`E_\mathrm{a}`$are the symmetric and antisymmetric energy parts resp., and $`E`$ is the total energy of in toroidal field component. Figures 17 and 18 show the temporal development of parity of the toroidal magnetic fields in the Models II and V along with a repeated plot of $`\alpha _{\mathrm{SS}}`$.

We should emphasize that the initial perturbation of the magnetic field had dipolar symmetry. Since the total flux through the vertical surfaces, however, vanishes, field configurations with either symmetry may arise from the evolution of the model. The dipolar initial perturbation may not have been fully re-organized as the diffusive time-scale

$$\tau _{\mathrm{diff}}=H^2/\eta $$
(22)
is 1600 rotation periods for the low-diffusivity model II and 160 revolutions for the high-diffusivity model V. The latter covered $`0.12\tau _{\mathrm{diff}}`$. An additional experiment used a completely mixed initial magnetic-field perturbation according to a distribution which has zero parity in $`B_z`$, and the toroidal component immediately emerging from the slight vertical has zero parity, too. Again, the magnetic flux through the vertical surfaces is zero. Such initial conditions have been applied to Models II and V. The runs lasted roughly 10 revolution resulting in parity variations with a slight tendency towards quadrupolar (symmetric) solution ($`P`$ from +0.1 to +0.2).

### 4.3 The dynamo of the disk model

Average kinetic helicities were computed based on a layer slightly below one density scale-height in both hemispheres. The average is taken in $`\varphi `$\- and $`r`$-direction. There is a considerable scatter in the time series of the helicity, but temporal averages are all negative on the northern side for all models. The only exception is the northern part of Model II with a near-zero value. Negative sign also holds for the two test models with mixed-parity initial perturbation. Average current helicities are also negative on the northern side (and positive on the southern one) throughout almost all models. The only exceptions are the mixed-parity models which show a current helicity with opposite sign compared with the kinetic helicity in both hemispheres.

Since the angular momentum transport is dominated by magnetic stresses, we also ask about the energy in the flow fluctuations and the magnetic field. The temporal evolutions of these two quantities in Figures 20 and 21 show a clear difference between the Models II and V. Kinetic-energy dominance holds for Model II, whereas magnetic dominance is found for Model V.

In order to evaluate the applicability of the $`\alpha ^{\mathrm{dyn}}`$ approach, we compute the average toroidal field $`B_\varphi `$ and compare it with the actual average electromotive force $`_\varphi =(u\times B)_\varphi `$ derived directly from the simulated vector fields. A correlation between the two quantities may justify the $`\alpha `$-effect for dynamo generation of magnetic fields. As a meaningful $`\alpha ^{\mathrm{dyn}}`$ must change its sign at the equator, we plot the correlation for both hemispheres of the disk separately; the result is given in Figure 23. As the temporal fluctuations are strong, we plot the time-averages of $`B_\varphi `$ and $`(u\times B)_\varphi `$ in Figure 22 only with their resulting sign. The quantities are also averaged in azimuthal direction and plotted in the $`(r,z)`$-plane. White areas denote a negative sign; the dominance of negative averaged electromotive forces indicates a positive $`\alpha ^{\mathrm{dyn}}`$ in the northern hemisphere in accordance with Figure 23.

The classical understanding of the Coriolis force giving preference to left-handed helicity on the northern hemisphere (whence positive $`\alpha ^{\mathrm{dyn}}`$) appears not applicable for the correlation plots of Models II and III. The fluctuations were strong, and no meaningful sign of $`\alpha ^{\mathrm{dyn}}`$ can be derived. However, Models Ia, V, and VI show a mostly negative averaged EMF for both hemispheres, while $`B_\varphi `$ changes its sign. An indication for positive $`\alpha ^{\mathrm{dyn}}`$ is thus found from these simulations.

In fact, Models V and VI are those which show saturated kinetic energies. The other significant difference to the indefinite Models II and III is the 10times higher diffusivity, $`\eta =0.01`$. Notwithstanding, an indication for a negative $`\alpha ^{\mathrm{dyn}}`$ was already found in simulations of a local box cut out of the disk by Brandenburg et al. (1995). The same sign is found by Ziegler & Rüdiger (2000b) from long runs of local simulations, although a clear influence of resistivity on the evolution was found there. Lowest magnetic Reynolds numbers are $`\mathrm{Re}_\mathrm{m}=1450`$ for Model V in this Paper, calculated with the velocity difference between inner and outer radial boundary, just as Ziegler & Rüdiger did for the local box.

These considerations are somewhat limited, since a straight-forward regression line in Fig. 23 will have a significant offset from the origin of the graph; a vanishing $`B_\varphi `$ does apparently not coincide with a vanishing $`(u\times B)_\varphi `$. The effect of the other field components is considered small though; while the average hemispherical $`B_\varphi `$ are of the order of 100, the $`B_r`$ is of order 10 and $`B_z`$ of order unity. A considerable contribution is suspected from the current density through the $`\eta _{ijk}`$-term in (17), and is appears to be quite natural that an offset in the correlation graphs in Figure 23 is found. For this reason, we did not plot such regression lines in the Figure.

Additionally, the sole consideration of the induction equation in the “classical” dynamo theory leads to the omission of the Lorentz force, which is not only essential for the onset of the magneto-rotational instability but for the maintenance of the turbulence and the transport of angular momentum as well. It is thus obvious that the neglect of the back-reaction of magnetic fields on the flow need not lead to representative models of Keplerian disks.

## 5 Conclusions

Three-dimensional global simulations have shown, that the magnetic shear-flow instability is a fast mechanism to generate a turbulent flow in a Keplerian disk. We presented models with a computational domain covering the full azimuthal range, accounting for accretion of matter. The Shakura-Sunyaev parameter $`\alpha _{\mathrm{SS}}`$ increases rapidly during the first revolutions after switching on the magnetic field and reaches $`10^2`$$`10^1`$ at this stage of the computations, though it appears not yet fully saturated in several of the models. Maxwell stresses exceed Reynolds stresses almost entirely. The latter undergo strong fluctuations and are partly negative.

Indications for a dynamo action in the disk are found, the corresponding dynamo-$`\alpha ^{\mathrm{dyn}}`$ tending to be negative north of the equatorial plane and positive south of the equatorial plane. The generated magnetic fields may maintain the turbulence even if the external field ceases. The strong excitation of low-order azimuthal modes in the magnetic field is another promising fact for dynamo action with respect to the Cowling theorem.

### 5.1 The dynamo of Model V

For various reasons, a representative example for an accretion disk dynamo is provided by Model V (and nearly as well by the similar Model VI with higher inflow velocity). The time series is sufficiently long covering more than 18 revolution periods. The angular momentum transport reaches a saturated level of high efficiency during the last four orbital periods. The magnetic diffusivity is roughly two orders of magnitude higher than the numerical diffusivity. The striking dominance of one sign of the averaged EMF in Figure 22 indicates a physically evolved, relevant state of the system.

The results connected with the outcomes of mean-field dynamo theory include include the following facts: (i) a negative average toroidal magnetic field is found for the northern hemisphere, a positive on the southern one. The averaged EMF is negative in both hemisphere indicating a positive $`\alpha ^{\mathrm{dyn}}`$-effect on the northern side (negative on the southern side). (ii) The correlation plot of average toroidal field and EMF gives the same picture. (iii) negative kinetic and magnetic helicities on the northern hemisphere (positive on the southern one) are also consistent with a positive $`\alpha ^{\mathrm{dyn}}`$. (iv) The dipolar structure of the solution contrasts with the results from mean-field dynamo models of disks which yield quadrupolar solutions for positive $`\alpha ^{\mathrm{dyn}}`$. We have to add though that the disks of our Paper are not thin. Mean-field simulations by Covas et al. (1999) show a transition from quadrupolar to dipolar fields when the opening angle of the disk is enlarged, rather representing a torus than a disk.

The extension of the run-times of such simulations exceeding the order of diffusion times promises further results about dynamo action in accretion disks. The connection with mean-field concepts are most interesting as well as the implications for jet-launch models.

###### Acknowledgements.

The invaluable help by D. Elstner, Potsdam, in adapting the ZEUS-3D code for our purposes is greatly acknowledged. We are grateful to the John v. Neumann-Institut for Computing at the Forschungszentrum Jülich, Germany, for the opportunity to use the Cray T90 computer. R.A. thanks for the kind support by the Deutsche Forschungsgemeinschaft.


format: text

text:
# Finding Galaxy Clusters using Voronoi Tessellations

## 1 Introduction

Wide field imaging is becoming increasingly common since new large format CCD cameras are, or soon will be available at several telescope (see e.g. MEGACAM, Boulade 1998, Boulade et al. 1998, WFI, Baade 1999, etc). The possibility to perform wide field imaging of the extragalactic sky allows a systematic search of medium-high redshift galaxy clusters in two-dimensional photometric catalogs of galaxies. These candidate clusters are of cosmological interest and are primary targets for subsequent follow-up spectroscopical observations (see, e.g., Holden et al. 1999, Ramella et al. 2000).

Several automated algorithms have already been developed for the detection of clusters within two dimensional galaxy catalogs. The classical techniques used for this task are the “box count” (Lidman & Peterson 1996) and the “matched filter” algorithm proposed by Postman et al. (1996, hereafter P96) and its recent refinements (see Kepner et al., 1999, Kawasaki et al. 1998, Lobo et al 1999).

The box-counting method uses sliding windows (usually squares) which are moved across the point distribution marking the positions where the count rate in the central part of the window exceeds the value expected from the background determined in the outermost regions of the window. The main drawbacks of the method are the introduction of a binning to determine the local background, which improves count statistics at the expense of spatial accuracy, and the dependence on artificial parameters like bin sizes and positions or window size and geometry.

The “matched filter” is a maximum-likelihood (ML) algorithm which analyzes the galaxy distribution with the assumption of some model profiles to fit the data (e.g. a density distribution profile and a luminosity function). This last technique has been used to build the Palomar Distant Cluster Survey (PDCS, see P96) catalog and the EIS cluster catalog (Olsen et al. 1999, Scodeggio et al. 1999). These two catalogs are two of the largest presently available sets of distant clusters, with 79 and 302 candidate clusters respectively.

However, the main drawback of the matched filter method is that it can miss clusters that are not symmetric or that differ significantly from the assumed profile. This can be a serious problem since we know that a large fraction of clusters have a pronounced ellipticity (see e.g. Plionis et al. 1991, Struble & Ftaclas 1994, Wang & Ulmer 1997, Basilakos et al. 2000). Furthermore, the matched filter technique is sensitive to border effects because the ML is computed on a circular area. This is a problem because real galaxy catalogs are finite and usually present holes in regions corresponding to camera defects or bright stars.

Other interesting methods to detect overdensities in a galaxy catalog are the LRCF method (Cocco & Scaramella 1999), kernel based techniques (Silvermann 1986, Pisani 1996), and wavelet transforms (Bijaoui 1993, Fadda et al. 1997). In general, all these techniques are very sensitive to symmetric structures and suffer from border effects.

Ebeling & Wiedenmann (1993) develop a method based on Voronoi tessellation (VTP) in order to identify overdensities in a Poissonian distribution of photon events. Their goal is to detect sources in ROSAT x-ray images. The VTP does not sort points into artificial bins and does not assume any particular source geometry for the detection process. These interesting properties are also very attractive for other applications. For example, Meurs and Wilkinson (1999) present an application of VTP to a slice of the CfA redshift survey with the purpose of isolating the bubbly and filamentary structure of the galaxy distribution. Several other applications of the Voronoi tessellation to astronomical problems have been published (see, e.g., van de Weygaert 2000, El-Ad & Piran 1996, 1997, Ryden 1995, Goldwirth et al. 1995, Zaninetti 1995, Doroshkevich et al. 1997, Ikeuchi & Turner 1991, Icke & van de Weygaert 1991). Clearly, Voronoi tessellation is also interesting for the search of galaxy clusters (Ramella et al. 1999, Kim et al. 2000).

In this paper we present an automatic procedure for the identification of clusters within two-dimensional photometric catalogs of galaxies. Our technique, announced in Ramella et al. 1999, is based on the VTP of Ebeling & Wiedmann (1993). The procedure, hereafter Voronoi Galaxy Cluster Finder (VGCF), is completely non-parametric, and therefore sensitive to both symmetric and elongated and/or irregular clusters. Moreover, the VGCF is not significantly affected by border effects, and automatically assigns members to structures.

In Sect. 2 we describe the method and in Sect. 3 we show an application of the VGCF to the $`\alpha =13^h26^m,\delta =+30^o52^{^{}}`$ field of the Palomar Distant Cluster Survey (see P96). In particular, we illustrate the performances of the algorithm and discuss some aspects of the its application that are specific to the construction of a catalog of galaxy clusters. Finally, in Sect. 4 we summarize our results.

In what follows we assume an Hubble constant $`H_0=100kms^1Mpc^1`$ and a deceleration parameter $`q_0=1/2`$.

## 2 The method

A Voronoi tessellation on a two-dimensional distribution of points (called nuclei) is a unique plane partition into convex cells, each of them containing one, and only one, nucleus and the set of points which are closer to that nucleus than to any other (see Fig. 1). For us, any point is a galaxy of a given galaxy catalog. The algorithm used here for the construction of the tessellation is the “triangle” C code by Shewchuk (1996). In order to avoid border effects, we restrict our study to Voronoi cells which do not intersect the convex hull of the set of points. In the case of galaxy catalogs with holes, we do not consider Voronoi cells that intersect the hole borders. To do that, we model the hole with a combination of rectangles and we reject cells which have at least one vertex inside one or more of these rectangles.

Since we want to find overdensities in a catalog of galaxies, the interesting quantity for us is the density, i.e. the number of galaxies for area unit, rather than the cell area associated by the Voronoi tessellation to each galaxy. An empirical distribution for randomly positioned points following Poissonian statistics, has been proposed by Kiang (1966):

$$dp(\stackrel{~}{a})=\frac{4^4}{\mathrm{\Gamma }(4)}\stackrel{~}{a}^3e^{4\stackrel{~}{a}}d\stackrel{~}{a},$$
(1)
where $`\stackrel{~}{a}a/<a>`$ is the cell area in units of the average cell area $`<a>`$. As each cell contains exactly one galaxy the corresponding density is the inverse of the cell area $`f1/a`$. The idea of Ebeling & Wiedenmann (1993) is to estimate the background by fitting the Kiang cumulative distribution to the empirical cumulative distribution resulting from the real catalog in the region of low density which is not affected by the presence of structures ($`\stackrel{~}{f}f/<f>0.8`$). We can thus establish a density threshold (Fig. 2) and define as overdensity regions those composed by adjacent Voronoi cells with a density higher than the chosen threshold.

In this paper, the density threshold is to separate background regions and fluctuations which are significant overdensities at the 80% c.l.. Since some of the selected overdensities can still be random fluctuations rather than physical clusters, it is necessary to further suppress random overdensities in order to save only significant clusters. We proceed as follows: we compute the probability that an overdensity corresponds to a background fluctuation by:

$$N(\stackrel{~}{f}_{min},n)=N_{bkg}\frac{n_{src,fluc}(\stackrel{~}{f}_{min},0)e^{b\stackrel{~}{f}_{min}n}}{b(\stackrel{~}{f}_{min})},$$
(2)
where $`n`$ is the number of points associated to the cluster, $`N_{bkg}`$ is the number of background galaxies, and $`\stackrel{~}{f}_{min}`$ is the minimum normalised density of the cells belonging to the cluster. Both $`n_{src,fluc}`$ and $`b`$ are functions of $`\stackrel{~}{f}_{min}`$, which can be computed from simulations (Ebeling & Wiedemann, 1993). We reject overdensities whose probability to be a random fluctuation of the background is greater than 5%.

Once we detect significant clusters, we regularize their shape. We assume that the points inside the convex hull defined by the set of points belong to the cluster itself. Then we fit a circle to each cluster and we expand it until the mean density inside the circle is lower than the density of the original cluster. We perform this expansion of the original cluster because the external points of a cluster have usually large area cells, typical of the low density background. These cells can not be associated from the beginning to the cluster by the Voronoi tessellation algorithm.

We stress that fitting a circle to the Voronoi cluster has absolutely no part in the detection of the cluster. It only provides a convenient way to catalog the cluster with a center and a radius. In principle, we could fit an ellipse rather than a circle, since our cluster finder is sensitive to elongated clusters. In practice this is not convenient because, as Ebeling (1993) shows, fitted ellipse eccentricities become insignificant for values lower than about 0.7.

At this point we have a list of clusters with position, number of members associated (we statistically subtract the background), contrast against the background, and projected size on the sky. We also have a file of galaxies with position, magnitude, area of the Voronoi tessel, and a label for appartenence to one of the detected clusters. As a final (optional) step, we can assume a density model and fit it on the galaxy distribution at the position where the cluster has been found. In this way we can have a better estimate of the number of members and of the sizes of the cluster, in case clusters are well described by the model.

## 3 Building a cluster catalog: the case of a PDCS field

In this section we test the VGCF on part of the galaxy catalog used by P96 to produce the Palomar Distant Cluster Survey (PDCS).

In particular, we run the VGCF on the galaxy catalog derived from the $`V_4`$ band image of the PDCS field centered at $`\alpha =13^h`$ 26m and $`\delta `$ = +$`29^o`$ 52’ (J2000). The galaxy catalog contains 25432 galaxies. For each galaxy the catalog lists right ascension and declination, and total magnitude in the $`V_4`$ band (see P96, Sect. 2). The effective sky area covered by the field is 1.062 square degrees (see Fig. 3).

We also have a catalog of bright stars areas to be avoided by the VGCF. As already pointed out, the VGCF can very easily deal with “holes” in the galaxy catalog.

### 3.1 Magnitude bins

The PDCS galaxy catalog is deep, i.e. the number density of galaxies is high. It is therefore a priori possible that only the highest number density contrasts survive the “dilution” caused by foreground/background galaxies. For this reason we decide to run the VGCF in magnitude bins. As we will show, to run the algorithm in magnitude bins also helps in further suppressing spurious clusters.

The choice of the bin size is somewhat arbitrary and cannot be equally optimal for the detection of all clusters. In fact, the magnitude range of the excess counts produced by clusters varies with the cluster distance and richness. We choose to optimize the bin size for the detection of PDCS-like clusters at intermediate redshift ($`z0.3\mathrm{\hspace{0.17em}0.6}`$), i.e around the peak of the galaxy selection function.

In Fig. 4 we plot the magnitude distribution (normalized to unit area) of galaxies within the projected area of P51 ($`z_{PDCS}0.4`$ and richness $`N_R80`$, see P96 for the definition of $`N_R`$), superimposed on the average counts per square degree of the whole PDCS field. From inspection of Fig. 4 it appears that a reasonable choice for the bin size is two magnitudes. With such a bin most of the counts of the cluster peak fall within one bin and the cluster detection is optimal. In fact, a wider bin would decrease the contrast of the cluster counts, while a narrower bin would increase the statistical uncertainty of the detection because of the lower number of cluster counts.

In order to further justify our bin size, we build simulated clusters at different redshifts ($`z=`$0.3, 0.5 and 0.8 respectively). Our simulated clusters have a projected density profile with a slope $`r^1`$ outside a core radius of 0.1 Mpc $`h^1`$ and out to a cut-off radius of 1 $`h^1`$ Mpc. This choice is motivated by the fact that PDCS clusters have, on average, a slope $`r^{1.4}`$ (Lubin & Postman 1996) and that steeper slopes produce clusters that are more easily detected by the VGCF. We simulate the magnitudes of cluster galaxies by sampling a Schechter luminosity function with a slope $`\alpha =1.2`$ and an absolute characteristic magnitude $`M_V^{}=20`$. We apply both evolutionary and K-corrections to cluster galaxies according to, respectively, Pence (1976) and Poggianti (1996). We also take into account selection effects in the visibility of galaxies (Phillips et al. 1990). We set the richness of our simulated clusters according to the P96 definition of the richness parameter $`N_R`$ (P96 clusters typically vary from $`N_R=10`$ to $`N_R=90`$). As an example, in Fig. 5, we plot both the counts within the area of a simulated cluster at z=0.3 (radius of 1 Mpc $`h^1`$ and $`N_R=60`$) and the average counts of the PDCS field expected within the same area. From the count distribution of our simulated clusters, we find that a bin size of two magnitudes is indeed adequate to contain the peak of counts corresponding to clusters within the redshift range $`z0.3\mathrm{\hspace{0.17em}0.6}`$.

In the case of real clusters, the position of the peak of counts in the magnitude distribution is not known a priori. Therefore, we “slide” at small steps the magnitude bin over the whole magnitude range of the catalog and run the VGCF at each step. We set the step to 0.1 magnitudes, the minimum possible step considering the typical uncertainties of the PDCS magnitudes.

For simplicity, we start the VGCF with the bin $`18.0V_4<20.0`$ and stop at the faint end of the survey, i.e. with the bin $`21.8V_4<23.8`$. We neglect 76 galaxies brighter than $`V_4=18.0`$ without consequences for our identification of clusters within the redshift range $`z0.3\mathrm{\hspace{0.17em}0.6}`$.

### 3.2 Clusters and their detections in bins

Before producing our final cluster catalog, we need to define a criterium to identify clusters by associating fluctuations detected in different magnitude bins. The criterium will consists of a maximum projected distance between the centers of fluctuations to be associated and of a minimum number of coincident fluctuations required for a positive detection, $`N_{min}`$. In fact, because of the statistical noise of the foreground/background galaxy distribution, the centers of the fluctuations produced by a cluster in different magnitude bins will be slightly different. As far as the number of coincident fluctuations produced by a cluster is concerned, it will depend on the cluster distance, richness and luminosity function. Clearly the choice of the criterium has to be made having in mind the goal of the detection algorithm and the characteristics of the galaxy catalog.

In order to define the criterium of association for our test application, we peform extended tests both on the PDCS field and on simulations of Poissonian fields with embedded simulated clusters. Our simulated fields have the same general properties of the PDCS field. In particular, we run the VGCF on 100 catalogs each containing 18 simulated clusters (a number similar to the number of clusters in the PDCS field) embedded within a Poissonian galaxy field. Clusters are simulated as in the previous subsection. All simulated catalogs contain 25432 galaxies, i.e. the same number of PDCS galaxies.

As a result of our tests, we consider coincident two fluctuations with centers separated on the sky by a projected distance $`d_{12}0.3min(R_1,R_2)`$, where $`R_1`$ and $`R_2`$ are the radii of the two fluctuations. A tighter criterium would break the sequence of fluctuations corresponding to a real cluster, a looser criterium would incorrectly associate fluctuations produced by adjacent clusters/fluctuations to the same cluster.

We now set the minimum number of fluctuations, $`N_{min}`$, required for the detection of a cluster. In Fig. 6 we plot the average number of spurious fluctuations mis-identified as clusters as a function of $`N_{min}`$. The number of spurious clusters drops dramatically as $`N_{min}`$ increases from 1 to 5. For $`N_{min}=5`$ there are on average 1.5 spurious clusters per field. This number decreases slowly as $`N_{min}`$ increases further. This result indicates that $`N_{min}=5`$ is a a good choice to keep the number of Poisson fluctuations low while still being sensitive to poor or distant clusters.

In Fig. 7 we show the variation with $`N_{min}`$ of the detection efficiency of simulated clusters at redshifts z=0.3 (panel a), 0.5 (panel b), and 0.8 (panel c). At each redshift the different curves correspond to richnesses ranging from $`N_R=10`$ to $`N_R=60`$. Clearly the fraction of detected clusters decreases as $`N_{min}`$ increases.

We note that the curve in Fig. 6 corresponds to an evaluation of the False Positive Rate (FPR), i.e. the probability of detecting as a cluster a random fluctuation of the galaxy distribution. For $`N_{min}=5`$, our FPR is very low. By comparison, P96 give 4.2 spurious detections per square degree with peak signal greater than $`3\sigma `$. The curves in Fig. 7 correspond to a measure of the False Negative Rate (FNR).

### 3.3 The VGCF run and comparison with P96

We present here the results of the run of the VGCF on the $`V_4`$ catalog of the PDCS field at $`\alpha =13^h`$ 26m and $`\delta `$ = +$`29^o`$ 52’ (J2000). As discussed in the previous section, we run the VGCF in bins two magnitude wide, “sliding” with 0.1 magnitude steps within the magnitude range $`18.00V_4<23.8`$. In total we run the VGCF in 39 magnitude bins and identify as clusters at least five fluctuations that, according to our criterium, are angularly coincident (see previous subsection).

On a Sun ULTRASparc 30 workstation, the time required to run the whole procedure is about 40 minutes.

Our output cluster list consists of 37 objects. We characterize each cluster with the properties of the fluctuation with the highest signal-to-noise ratio as estimated by the ratio of the number of cluster galaxies to the square root of the number of background galaxies expected within the cluster area.

In Table 1, for each cluster we list: 1) identification number, 2) J2000 right ascension and 3) J2000 declination, 4) radius, 5) the cluster signal-to-noise ratio, 6) the estimated number of cluster galaxies and 7) the number of expected background galaxies, 8) the central magnitude of the bin where we detect the cluster with the highest signal-to-noise ratio, 9) the cross-identification with the PDCS catalog.

In Fig. 8 we plot circles on the sky corresponding to our clusters (solid lines). We label our clusters with a “V” followed by their order number in Table 1. In Fig. 9, we give a graphic summary of all the fluctuations of each cluster. The abscissa is the order number of the magnitude bin of the fluctuation and each cluster is represented by a row of circles parallel to the magnitude bin axis. The radii of the circles are scaled with the signal-to-noise ratio of the detection.

We remind here that we fit a circle to a fluctuation after its detection, and that the only use of the circle is to provide a convenient way to catalog the cluster with a center and a radius.

The richest and more reliable clusters exhibit, in Fig. 9, a sequence of fluctuations. The signal-to-noise ratio of these fluctuations regularly increases up to a maximum and then decreases. Several fainter clusters show the same behavior, although at a generally lower S/N level. Clearly, the position of the fluctuations along the magnitude axis is related to the cluster distance.

Some clusters, for example V7 and V12, display substantial gaps along the sequence. The suspicion is that these clusters may actually consist of unrelated overdensities projected along the same line-of-sight.

We now proceed to a detailed comparison with P96 clusters. The original PDCS catalog consists of 19 clusters detected by the matched filter algorithm both in the $`V_4`$ band and in the $`I_4`$ band (see P96, Table 4). Furthermore, the PDCS catalog lists 7 clusters detected only in the $`V_4`$ band and added to the catalog after visual inspection of the CCD images (see P96, Table 5). In total there are 26 P96 $`V_4`$-band clusters, 18 of which are detected at the 3$`\sigma `$ level. Of these last clusters, 13 have estimated redshift within the range where our VGCF is optimal ($`0.2z0.6`$). In Fig. 8 we plot circles corresponding to all the P96 clusters detected by the matched filter technique (dotted line).

We identify 12 of the 13 clusters detected by P96 at the 3$`\sigma `$ level and with estimated redshifts from $`z=0.2`$ to $`z=0.6`$. The only cluster we do not identify, P54 ($`z=0.5`$), falls just under the limit of our selection criterium with only four coincident fluctuations.

Beside the above 12 clusters, we detect other 25 galaxy systems. Of these 25 systems, 2 are P96 clusters with confidence level less than 3$`\sigma `$ and redshift $`z0.6`$. The remaining 23 do not have a counterpart in the PDCS catalog.

Most of the 23 new clusters have properties comparable to those of P96 clusters as detected by the VGCF. At least part of the new systems have been detected by the VGCF thanks to its specific properties. For example, because the VGCF does not smooth the data, it is able to separate two large overdensities (V2 and V18) that are unlikely to be one system (P49). In fact V2 (which we identify as the counterpart of P49) and V18 are well separated both on the sky and in magnitude.

The property of the VGCF of detecting overdensities irrespective of their shape is probably the key for the detection of V26, a very elongated system. In Fig. 10 we show the galaxies of V26 as detected by the VGCF before the regularization processes (black triangles). We also show as grey tessels those of galaxies added at the end of the regularization processes and as white tessels those of background galaxies. The plot refers to the tessellation of the PDCS field in the magnitude bin of maximum contrast for V26.

As far as the richness, $`N_R`$, of VGCF clusters is concerned, we compute it almost in the same way as P96, the only difference being that we use a fixed angular radius of 0.1 deg within which to count cluster galaxies. For systems within the redshift range 0.3-0.6, our angular radius roughly corresponds to a linear radius of 1 $`h^1`$ Mpc. As P96 shows, $`N_R`$ is related to the Abell richness of the system.

In Fig. 11 we plot the distribution of $`N_R`$ of our new systems (solid histogram) and of P96 clusters as detected by the VGCF (dotted histogram). In general the richness of the new systems is similar to that of P96 clusters (as detected by the VGCF), but the VGCF detects a larger fraction of poor systems than the matched filter algorithm. In fact, we have 5 systems with $`N_R<20`$ while only 1 P96 cluster falls within this range. We probably find these poorer systems because we do not smooth the data.

That the general properties of new VGCF systems are not very different from P96 systems is also evident from Fig. 11 where we label with a “P” and the PDCS identification number the VGCF identifications of P96 clusters.

A more detailed comparison between our and P96 clusters is not useful since, in the end, only spectroscopic observations will allow a real comparison between the performances of the two algorithms. What we would like to stress here is that the VGCF is an algorithm that can usefully complement the matched filter (and probably other parametric cluster searches) for the definition of a complete and reliable 2D cluster catalog.

In many problems we think that the VGCF could be preferred over parametric searches and, in these cases, the example we discuss in this section shows that the VGCF retrieves most of the reliable identifications of the parametric searches.

## 4 Summary

We present an objective and automated procedure for detecting clusters of galaxies in imaging galaxy surveys. Our Voronoi Galaxy Cluster Finder (VGCF) uses galaxy positions and magnitudes to find clusters and determine their main features: size, richness and contrast above the background. The VGCF uses the Voronoi tessellation to evaluate the local density and to identify clusters as significative density fluctuations above the background. The significance threshold needs to be set by the user, but experimenting with different choices is very easy since it only requires a new selection from the output of the Voronoi tessellation and not a whole new run of the algorithm. The VGCF is non-parametric and does not require a smoothing of the data. As a consequence, clusters are identified irrispective of their shape and their identification is only slightly affected by border effects and by holes in the galaxy distribution on the sky.

The algorithm is fast, and automatically assigns members to structures. For example, a run on about 25000 galaxies only requires 7 minutes on a Sun ULTRASparc 30 workstation.

We perform test runs of the VGCF both on simulated galaxy fields and on the 25432 galaxies identified in the $`V_4`$ band image of the PDCS field centered at $`\alpha =13^h`$ 26m and $`\delta `$ = +$`29^o`$ 52’ (J2000, see P96). Given the depth of the PDCS survey, we run the VGCF in 39 overlapping 2 magnitude wide bins, from $`V_4=18.00`$ to $`V_4=23.75`$. For the detection of a cluster we require at least five angularly coincident fluctuations. Based on 100 simulations of the PDCS field, we find that there are, on average, 1.5 spurious detections in a simulated PDCS field with Poisson background. Our choice of bin size and number of bins, together with the minimum number of fluctuations required for a cluster detection, optimizes the VGCF for the redshift range $`0.3<z<0.6`$.

We find 37 clusters, 12 of which are VGCF counterparts of the 13 PDCS clusters detected at the 3$`\sigma `$ level and with estimated redshifts from $`z=0.2`$ to $`z=0.6`$. Of the remaining 25 systems, 2 are P96 clusters with confidence level $`<3\sigma `$ and redshift $`z0.6`$.

23 VGCF clusters have no counterpart in the PDCS catalog. An inspection of their properties indicates that some of these clusters may have been missed by the matched filter algorithm for one or more of the following reasons: a) they are very poor, b) they are extremely elongated, c) they lie too close to a rich and/or low redshift cluster.

In conclusion, the VGCF can usefully complement the matched filter algorithm (and probably other parametric cluster searches) for the definition of a complete and reliable 2D cluster catalog. In fact, we think that the VGCF could be preferred over parametric searches in many cases since the example we discuss shows that the VGCF anyway retrieves the reliable identifications of parametric searches.

###### Acknowledgements.

We wish to thank Marc Postman for the galaxy catalog of the PDCS field we used to prepare this paper.