# Dynamical models of the elliptical galaxy NGC 4494

S. a. Rodionov, E. Athanassoula

► To cite this version:

S. a. Rodionov, E. Athanassoula. Dynamical models of the elliptical galaxy NGC 4494. Monthly Notices of the Royal Astronomical Society, 2011, 410, pp.111-126. 10.1111/j.1365-2966.2010.17429.x. hal-01428155

HAL Id: hal-01428155

https://hal.science/hal-01428155

Submitted on 30 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamical models of the elliptical galaxy NGC 4494

S. A. Rodionov<fnref n="1" />,<fnref n="2" />,<fnref n="3" />* and E. Athanassoula<fnref n="1" />

Accepted 2010 July 26. Received 2010 July 12; in original form 2009 December 19

## ABSTRACT

We present dynamical models of NGC 4494, which we built using our iterative method presented in a previous paper. These models are live N-body models consisting of equal-mass particles, and they are in steady state as confirmed by a fully self-consistent evolution. Our goals were of twofold. The first one – namely to test whether our iterative method could indeed be used to construct galactic models following given observational constraints, both photometric and kinematic – was fully achieved. Our method allowed us to go beyond a simple spherical model and to make full sets of rotating, axisymmetric models without any limitations to the velocity distribution. Our second goal was to understand the structure of NGC 4494 better, and more specifically to set constraints on its halo mass. For this we tried three families of models: without halo, with a light halo and with a heavy halo, respectively. Our models reproduce well the photometry and the kinematics, the latter except for specific regions where some non-equilibrium or non-axisymmetric structure could be present in the galaxy (e.g. the kinematically decoupled core). However, the lower-order moments of the velocity distribution (up to and including the second order) do not allow us to discriminate between the three haloes. On the other hand, when we extend the comparison to the higher-order moments of the velocity distribution obtained from the long-slit data, we find that our light halo model fits the data better than the no halo, or the heavy halo models. They also reproduce the shape of the angular dependence of the PNe velocity dispersion in the outermost parts of the galaxy, but not the amplitude of its azimuthal variation. This may imply that a more general class of models, such as triaxial, may be necessary for a better fit.

**Key words:** methods: numerical – galaxies: elliptical and lenticular, cD – galaxies: individual: NGC 4494 – galaxies: kinematics and dynamics.

## 1 INTRODUCTION

Dark matter around ordinary elliptical galaxies is one of the hottest topics in dark matter studies today. The main goal is to obtain sufficient constraints on the dark matter mass from observed stellar kinematics. Traditional long-slit absorption line spectroscopy can only very rarely give kinematics outside $2R_e$, where $R_e$ is the effective radius encompassing half the total light of the galaxy (see e.g. Coccato et al. 2009). It is, nevertheless, possible to obtain line-of-sight velocities at larger radii using planetary nebulae (PNe) because their strong emission line at 5007 Å [O III] stands out against the faint galaxy background. It is usually assumed that PNe trace the kinematics of the underlying field

stars.<fnref n="1" /> It is thus possible to obtain from the PNe the stellar kinematic parameters at the periphery of the galaxy, out to $5–7 R_e$ (Goudfrooij et al. 1994; Romanowsky et al. 2003; Douglas et al. 2007; de Lorenzi et al. 2008, 2009; Coccato et al. 2009; Napolitano et al. 2009, hereafter N09).

Romanowsky et al. (2003) studied the kinematics in the outer part of three ordinary elliptical galaxies, namely NGC 821, 3379 and 4494 (out to $4–6 R_e$), and found that their velocity dispersion profiles decline nearly Keplerian-like at radii outside $2R_e$. They modelled the observational data by means of spherical Jeans models and by means of orbit-based models [see Romanowsky et al. (2003) for details] and they showed that their data are consistent only with models with little or no halo. This result, in good

<footnote n="1">¹Laboratoire d'Astrophysique de Marseille (LAM), UMR6110, CNRS/Université de Provence, Technopôle de Marseille-Etoile, 38 rue Frédéric Joliot Curie, 13388 Marseille Cédex 20, France</footnote>

<footnote n="2">²Sobolev Astronomical Institute, St Petersburg State University, Universitetskij pr. 28, 198504 St Petersburg, Stary Peterhof, Russia</footnote>

<footnote n="3">³Center de recherche INRIA Bordeaux – Sub-Ouest, Bordeaux, 351 Course de la Libération, France</footnote>

<footnote n="1">¹Dekel et al. (2005), however, noted that observations of PNe can be biased towards the kinematics of younger stars.</footnote>

<footnote>*E-mail: seger@mail.ru</footnote>
agreement with what was already found for NGC 3379 by Goudfrooij et al. (1994) and for NGC 4697 by Méndez et al. (2001), is very surprising. Indeed, present theoretical and observational data argue that ellipticals are formed from mergings of spirals, which are known to have a considerable amount of dark matter (Bosma 2004, and references therein). So, if the progenitors have a considerable amount of dark matter, how can the merger product not have it? Furthermore, this result is in conflict with the predictions of the standard Λ cold dark matter ΛCDM cosmology.

Dekel et al. (2005) constructed elliptical galaxy models from numerical simulations of mergers of a pair of disc galaxies. Their resulting models have a ‘normal’ massive dark halo and a velocity distribution with a high radial anisotropy in the outer parts. The latter leads to a low observed (from most viewing angles) velocity dispersion in the outer parts, which in turn leads to a low estimated halo mass, contrary to the real dark halo mass of the model, which is normal. In this way, Dekel et al. (2005) explain the results of Méndez et al. (2001) and of Romanowsky et al. (2003) as due to the velocity distribution in the merger remnant.

Since this velocity anisotropy is so crucial to the data interpre-
tation, Athanassoula (2005a) examined whether it was general, or
whether it depended on the specific mergers used. She examined
the velocity anisotropy in multiple mergers, as would occur e.g. in
groups. In such cases, a pair merger is not examined as an isolated
event, but a whole sequence of mergers is considered. This model
would be more realistic in groups, but is also in good agreement
with the standard ΛCDM cosmology. The result of such mergers is
compatible with the observed properties of elliptical galaxies (Weil
& Hernquist 1994, 1996; Athanassoula & Vozikis 1999). Concern-
ing the velocity anisotropy in the outer parts, she found that the
result is more complex than the single pair mergers would predict
and that the anisotropy depends strongly on the time between two
successive mergers. Thus, more work is necessary to establish how
general the result of Dekel et al. (2005) is. A similar conclusion
was reached by Douglas et al. (2007) who modelled the data of
NGC 3379 and argued that there are considerable discrepancies be-
tween the observations and dark-matter-dominated simulations and
re-iterate the question of whether NGC 3379 has the kind of dark
halo that the current ΛCDM paradigm requires.

de Lorenzi et al. (2008, 2009) constructed dynamical models of NGC 4697 and 3379 using the χ²-made-to-measure particle method (Syer & Tremaine 1996; de Lorenzi et al. 2007) implemented in the *NMAGIC* code. Their main result is that the observational data are consistent with a fairly wide range of halo mass profiles, although it was possible to place some limits on the halo mass. For NGC 4697, de Lorenzi et al. (2008) found that models with a low-density halo with $v_c(5R_e) \lesssim 200 \text{ km s}^{-1}$ are not consistent with the data, where $v_c(5R_e)$ is the total circular velocity at $5R_e$. This, however, is a rather weak limit because even their model D with $v_c(5R_e) \approx 210$ (see fig. 15 in de Lorenzi et al. 2008) has a very light halo which contributes only about 35 per cent of the mass within $5R_e$. For the galaxy NGC 3379, de Lorenzi et al. (2009) found that a model without a halo, as well as a model with a heavy halo with $v_c(7R_e) \gtrsim 250 \text{ km s}^{-1}$ would be excluded by the observational data, but only at a $1\sigma$ confident level. So these constraints are not very strong.

The main problem in determining the halo mass profile from the
observed kinematics in elliptical galaxies is the well-known mass–
anisotropy degeneracy. The low-velocity dispersion on the periph-
ery of some elliptical galaxies can be explained either by nearly
isotropic models with a light halo or by radially anisotropic models
with a heavy halo. It is generally accepted that the mass–anisotropy
degeneracy can be broken by means of high-order moments of the

line-of-sight velocity distribution (LOSVD; Gerhard 1993; van der
Marel & Franx 1993). As shown in these papers, isotropic mod-
els have a Gaussian LOSVD, and radially anisotropic models have
centrally peaked LOSVD as well as long tails. So, these models
can be distinguished by means of high-order moments of LOSVD.
However, the highly radial anisotropic models presented in Dekel
et al. (2005) have a LOSVD with relatively weak deviations from
Gaussian (see supplementary information in Dekel et al. 2005).
Unfortunately, this means that, at least in some cases, breaking the
mass–anisotropy degeneracy can be very difficult, if not practically
impossible.

Let us also point out an obvious, but sometimes ignored, prob-
lem. From a mathematical point of view, it is possible to prove the
existence of a given type of model by simply constructing dynami-
cal models, but it is not possible to prove its non-existence. Let us,
for example, construct dynamical models of a real galaxy by means
of the *NMAGIC* method or by means of our iterative method (see
below). More specifically, let us construct an axisymmetric model
with some dark halo. If this satisfies all observational data, then
we have proven the existence of a model with such a halo agree-
ing with the observational data. But if we fail to construct such a
model, then formally we have not proven anything, since we cannot
exclude that our failure is due to the method itself, or to the fact
that we have not searched sufficiently. We would have needed to
prove that, if an equilibrium model with given parameters did exist,
then our method would construct it. This is not straightforwardly
proven for our iterative method, for the *NMAGIC* method or for any
other orbit-based method. Moreover, even if we did prove it, we
would have only proven that an axisymmetric model with this par-
ticular halo is excluded by the observational data. We would not
have proven anything about triaxial models, or about models with
a somewhat different halo mass profile. Consequently, the conclu-
sion of de Lorenzi et al. (2008, 2009) that the observational data of
NGC 4697 and 3379 are consistent with a wide range of halo masses
can be firmly believed. But, on the contrary, if one finds that a set of
models with a heavy dark halo do not agree with the observational
data, then one has at the best only an argument that the galaxies
in question have a light dark halo, or no halo. Unfortunately, this
is no proof, since another, more general, type of heavy halo could
perhaps have fitted the observations. It is thus very useful to try
different approaches to see whether this disagreement persists or
not. If more than one method lead to the same conclusion, then the
argument is considerably strengthened.

Initially, the present work was inspired by an article of N09, where
the authors presented a large amount of new observational data of
the ordinary elliptical galaxy NGC 4494, resulting in positions and
velocities of 255 PNe out to seven effective radii. They also pre-
sented new wide-field surface photometry from MMT/Megacam,
and long-slit stellar kinematics from VLT/FORS2. Using these data
they put constraints on the distribution of dark matter in this galaxy.
They constructed spherical dynamical models of the system using
two different methods, but both of them are based on Jeans equa-
tions. They argue that some dark matter is required by the data and
their best-fitting model has a relatively low halo mass.

In Rodionov, Athanassoula & Sotnikova (2009, hereafter RAS09) we presented an iterative method for constructing equilibrium *N*-body models with given properties. This method has already been widely used to construct initial conditions for *N*-body simulations (e.g. Rodionov & Orlov 2008; Machado & Athanassoula 2010) and can also be directly used for constructing dynamical models of real galaxies from observational data. In contrast with *NMAGIC* models, our models consist of particles with equal masses. They
are in steady state and can be directly used in N-body simulations. For example, we can directly check that the constructed model is indeed in equilibrium. Also, we can easily calculate for this model any parameter which can be directly obtained from the positions and velocities of the particles.

In this article, we will apply our iterative method to the construction of dynamical models of NGC 4494 from the observational data presented in N09. We have two purposes. First, our intention is to demonstrate that our iterative method can indeed be successfully used to construct dynamical models of real galaxies. The second aim is to return to the interesting question of the halo mass of NGC 4494, using a different method from that of N09, to see whether results can in any way depend on the method used. This is particularly important since the discrepancy between the no-halo model of N09 and the observational data is not large (see upper panels of their fig. 12).

We present the observational data in Section 2 and our method in Section 3. In Section 4 we describe our models and compare them to observations. We summarize and conclude in Section 5.

## 2 PREPARATION OF OBSERVATIONAL DATA

We use the same observational data as in N09. These include the surface photometry, the stellar kinematics along the major and minor axes obtained by means of long-slit spectroscopy, and velocities and positions of 255 PNe (see N09). When constructing our models, we use a physical system of units, i.e. kpc and M⊙. We adopt a distance to NGC 4494 of 15.8 Mpc (see N09), so that 1 arcsec is equal to 0.0766 kpc. Our models can be easily rescaled to any other distance. Let us assume we have an equilibrium model constructed for some distance $d_1$, and we want to rescale it to a distance $d_2 = Cd_1$. To keep the surface photometry unchanged we need to rescale all space coordinates of particles as $r_2 = Cr_1$. To keep the model in equilibrium we need to change the mass of the model as $M_2 = CM_1$. Particle velocities need not be changed, so all the observed kinematic parameters are unchanged. The new total luminosity of the galaxy is $L_2 = C^2L_1$, so that the new mass-to-light ratio is $\frac{M_2}{L_2} = \frac{1}{C} \frac{M_1}{L_1}$.

Let us now describe how we prepare the observational data for use in our iterative method.

### 2.1 Surface photometry

We use the combined V-band surface photometry of NGC 4494 presented in table A1 of N09. These data are a combination of HST-based observations of Lauer et al. (2005), ground-based CCD observations of Goudfrooij et al. (1994) and the new observations of N09.

We, furthermore, make the following simplifications. We assume that the ellipticity $\epsilon$ is the same at all radii and equal to 0.162, the mean value found in N09. We also assume that the shape of the isophotes is precisely elliptical, so as not to introduce in the analysis unconstrained high-order isophote shape parameters. The first two rows of table A1 in N09 give the surface brightness as a function of the intermediate axis $R_m$, measured in arcsec. This is related to the ellipticity and to the major axis $R_a$ by $R_m = R_a\sqrt{1-\epsilon}$. We convert surface brightness in units of $L_{\odot,v}/pc^2$ (assuming an absolute magnitude of the Sun in the V band $M_{\odot,v} = 4.8$) and $R_m$ in parsec using the adopted distance of 15.8 Mpc.

Excluding the innermost region ($R_m < 5$ arcsec), this surface brightness profile is fitted very well by the Sérsic law (Sérsic 1968)

$$I(R_m) = I_0 \exp(-(R_m/a_s)^{1/n}) \quad (1)$$

with parameters $I_0 = 41764 L_{\odot,v}/pc^2$, $a_s = 0.008809$ kpc and $n = 3.3$ as shown by N09.

Numerically, the surface brightness profile is described as follows. Inside the region $R_m < 0.46$ kpc ($\approx 6$ arcsec), we interpolate the tabular data linearly. In the region $0.46$ kpc $< R_m < 40$ kpc we adopt a Sérsic profile. In the region $40$ kpc $< R_m < 50$ kpc we truncate the Sérsic profile by means of a fifth-order polynomial (Dehnen 2000b, equation 4). From the profile of the surface brightness and the adopted value of ellipticity we can calculate the two-dimensional distribution of surface brightness. The total luminosity of this model is $L_v = 2.36 \times 10^{10} L_{\odot,v}$, or $M_v = -21.13$. We assume that the stellar mass-to-light ratio (M/L) is constant, in which case the distribution of surface brightness equals the surface mass distribution to within an unknown multiplier M/L.

### 2.2 Kinematical data

#### 2.2.1 Symmetries and system of coordinates

In order to prepare the kinematical data so that they can be used by our iterative method, we first need to define a system of coordinates which will be used for the model and to assume what symmetries the galaxy has.

In all the following, we will assume that the galaxy is axisymmetric. Elliptical galaxies can well be triaxial (Binney & Tremaine 2008). Triaxial models, however, have extra free parameters that add complexity to the modelling and are beyond the scope of this paper.

Let us consider a Cartesian $(X, Y, Z)$ coordinate system such that the sky plane coincides with the XZ plane. We choose the Z-axis so that it coincides with the minor axis of the projected image of the galaxy and then the X-axis will coincide with the major axis.<fnref n="2" /> We assume that the rotation axis is perpendicular to the X-axis and therefore in the YZ plane. This rotation axis is defined by its angle $\alpha$ with the Z-axis. If $\alpha = 0$ then the rotation axis coincides with the Z-axis, so the galaxy is edge-on.

We also assume that the galaxy has a reflection symmetry with respect to the plane of symmetry perpendicular to the rotation axis and centred on the centre of coordinates (0, 0, 0). In that case, the observed image of the galaxy will have a reflection symmetry about both the X and the Z axes. The observed line-of-sight velocity distribution will have a reflection symmetry about the X-axis. Also, the observed line-of-sight velocity distribution will have a reflection symmetry about the Z-axis, except for the velocity sign. It means that the points $(x, z)$ and $(-x, z)$ will have the same line-of-sight velocity dispersion and that their line-of-sight mean velocities will be equal but with opposite sign. The velocity distribution in the points $(x, z)$ and $(x, -z)$ will be fully identical. Consequently, if we know the line-of-sight velocity distribution in any of the four quadrants, we know it for the whole system. Similarly, if we assume such symmetry then all the observed kinematical data can be 'reduced' to the first quadrant of the sky plane $(x > 0, z > 0)$. So in a first stage, we reduce all kinematical data to the first quadrant.

In our iterative method the input kinematical data should be given as mean velocities and velocity dispersions in a set of two-dimensional areas on the sky plane and we need to present the observed velocities in this manner.

<footnote n="2">² These are only the major and minor axes of the projected image, and the real principal axes of the galaxy can, of course, be different.</footnote>
### 2.2.2 Long-slit spectra kinematics

Part of our kinematical data are obtained by means of long-slit spectroscopy (Coccato et al. 2009) from spectra taken along the major and the minor axes of the galaxy. The width of the slit was 1 arcsec. Thus, we have along each axis profiles of the rotation, of the velocity dispersion, as well as two Gauss-Hermite moments ($h_3$ and $h_4$; see Appendix A). When constructing the models we will use the rotational velocity profiles along the major axis and of the velocity dispersion along both principal axes (while the Gauss-Hermite moments will be used for the analysis of the so constructed models; see Section 4). We use these data almost 'as is' without any parametrization. We only bin the data suitably.

Let us describe how we prepare the profile of the mean velocity along the major axis. From Coccato et al. (2009) we get the table containing values of the mean velocity for different points along the major axis. So we have a set of pairs $x_i$, $v_i$, where $x_i$ is the coordinate of the data point on the major axis and $v_i$ is the observed line-of-sight mean velocity at this point. At first we 'reduce' these data to the first quadrant (see the previous section). For the mean velocity this implies multiplying $x_i$ and $v_i$ by $-1$ for each data point with $x_i < 0$. We need to convert these data into a set of two-dimensional areas with known mean velocity. These long-slit data were obtained from rather narrow zones along the major axis. We assign these data to a wider zone. We assign the data along the major axis to the area defined by $x < 0.5$ arcsec $||\varphi < 10^\circ$, where $\varphi$ is the angle between the current radius vector and the x-axis measured on the XZ plane, and $||$ is a logical OR. We divide this area along the x-axis into the pieces so that each long-slit data point corresponds to a single piece. We do this as follows. We sort the data points by $x_i$. For each data point we define two values as $l_i = (x_{i-1} + x_i)/2$ and $u_i = (x_i + x_{i+1})/2$. We set $l_1 = 0$ and $u_n = 2x_n - u_{n-1}$, where $n$ is the number of data points. For each data point we assign a two-dimensional area defined as $x > l_i$ && $x < u_i$ && ($x < 0.5$ arcsec $||\varphi < 10^\circ$), where && is a logical AND. We then join some of these areas in the following way. We create an N-body system with a surface distribution of particles obtained from the surface photometry (see Section 2.1). The number of particles is $N = 500000$. For each area we calculate the number of particles which are situated in this area. If the number of particles in some area is less than 1000 then this area is joined with a neighbouring area. The value of the velocity in each composite area is calculated as the mean value of the velocities of the constituent areas weighed by the number of particles in each area. This way of binning the observational data is rather unusual but it is convenient for us, because we make sure that in each area there is a sufficient number of particles.

The other profiles are prepared in the same way. The data along the minor axis are assigned to the area defined as $z < 0.5$ arcsec $||\varphi > 80^\circ$.

In our algorithm we do not use any information on the errors of the data. But we will use binned long-slit data for figures, so, in order to plot error bars, we need to calculate the errors of the binned data. We calculate the error of a binned datum as $e = \frac{1}{n^{3/2}} \sum e_i$, where $n$ is the number of original data points from which binned datum was calculated, and $e_i$ are the errors of these data points.

### 2.2.3 PNe kinematics

We use the PNe data in the outer part of the galaxy, which long-slit data cannot reach. More precisely, we use all PNe whose distance from the galaxy centre $R_{xz} = \sqrt{x^2 + z^2}$ is larger than 124 arcsec (9.5 kpc).

We 'reduce' the PNe data to the first quadrant (see Section 2.2). We divide the first quadrant into five zones with an equal opening angle $\varphi$ and use the $\varphi$ angle of the individual PNe to place them into the appropriate zone. The parameters of these five PNe groups are shown in Table 1. For each group we define an area ($\varphi^{(\text{min})} < \varphi < \varphi^{(\text{max})}$) && ($R_{xz}^{(\text{min})} < R_{xz} < R_{xz}^{(\text{max})}$). These areas are shown schematically in Fig. 1 and we will refer to these areas as A1, A2, A3, A4, A5.

For each group (area) we calculate the mean velocity and the velocity dispersion with the corresponding errors. We note that here we take the standard deviation as the dispersion. For the calculation of the standard deviation we use an unbiased estimator (Kenney & Keeping 1951, p. 171), although for such relatively large samples ($n > 10$, see Table 1), the use of an unbiased estimator is not essential.

We want to draw attention to two features of the PNe velocity distribution. The first feature is that the velocity dispersion at intermediate angles is noticeably higher than the velocity dispersion along the major or the minor axes. The velocity dispersion in area A3 is twice as high as in areas A1 or A5. This feature can be clearly seen in the right-hand panel of fig. 7 in N09, where, except for the inner areas, the iso-dispersion contours are elongated in a direction intermediate between the two principal axes. Since this figure was constructed under the assumption of triaxial symmetry, this peculiarity remains true independent of the axisymmetric assumption. The second feature is that the value of the mean velocity in A5 is negative and that it is low in the intermediate areas A2, A3 and A4. This feature can be explained by a twisting of the rotation axis in the outer part of the galaxy (see the left-hand panel in fig. 7 in N09 and corresponding discussion), but cannot be reproduced by the

<table><caption>**Table 1.** Parameters of the five PNe groups. The first column gives the name of the corresponding area (see Fig. 1) and the second and third ones give its lower and the upper azimuthal boundaries [$\varphi^{(\text{min})}$ and $\varphi^{(\text{max})}$, respectively]. The remaining columns give parameters relevant to the group of PNe in the corresponding area. Here $n$ is the number of PNe, $R_{xz}^{(\text{min})}$ and $R_{xz}^{(\text{max})}$ are the minimal and the maximal distance to the galaxy centre, $\bar{V}$ is the mean velocity, $\Delta\bar{V}$ is the error of the mean velocity, $\sigma$ is the velocity dispersion (standard deviation) and $\Delta\sigma$ is error of the velocity dispersion.</caption><thead><tr><th>Area</th><th>φ<sup>(min)</sup></th><th>φ<sup>(max)</sup></th><th>n</th><th>R<sub>xz</sub><sup>(min)</sup><br>(arcsec)</th><th>R<sub>xz</sub><sup>(max)</sup><br>(arcsec)</th><th>V̄<br>(km s<sup>-1</sup>)</th><th>ΔV̄<br>(km s<sup>-1</sup>)</th><th>σ<br>(km s<sup>-1</sup>)</th><th>Δσ<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>A1</td><td>0°</td><td>18°</td><td>11</td><td>166.3</td><td>297.1</td><td>72.7</td><td>15</td><td>49.7</td><td>11.2</td></tr><tr><td>A2</td><td>18°</td><td>36°</td><td>18</td><td>126.7</td><td>291.8</td><td>14.1</td><td>19.2</td><td>81.3</td><td>14</td></tr><tr><td>A3</td><td>36°</td><td>54°</td><td>17</td><td>131.1</td><td>378.9</td><td>21.9</td><td>25.5</td><td>105</td><td>18.8</td></tr><tr><td>A4</td><td>54°</td><td>72°</td><td>14</td><td>132</td><td>277.6</td><td>-2.64</td><td>24.5</td><td>91.8</td><td>18.2</td></tr><tr><td>A5</td><td>72°</td><td>90°</td><td>14</td><td>133.4</td><td>269.8</td><td>-24.7</td><td>14.3</td><td>53.4</td><td>10.6</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Five PNe areas. All PNe have been moved into the first quadrant.</figcaption></figure>

models with an axisymmetric velocity distribution. We therefore
cannot take into account in the iterative process the mean velocity
from all regions. Since in axisymmetric systems there is no rotation
along the minor axis and since the value from the A1 area agrees
very well with the long-slit data (Fig. 3d), we will choose to use the
mean velocity from area A1.

3 METHOD

**3.1 General outline**

We want to construct an equilibrium *N*-body model of an elliptical
galaxy from its observational data. As described in the previous
section, for the galaxy under consideration we have surface pho-
tometry, a distance estimate and various line-of-sight kinematics.
Assuming that the stellar *M*/L is constant, we can obtain from the
surface photometry and the distance the surface mass distribution
to within an unknown multiplicative constant *M*/L. In the case of an
*N*-body model, this implies that we have the projected surface dis-
tribution of particles, but the mass of the individual particles is not
known (in our models all particles have the same mass). We note,
however, that the *M*/L cannot be arbitrary, because it is related to
the line-of-sight velocity dispersion in the central part of the galaxy
which we know from observations (see Section 3.4). We also as-
sume that the galaxy is axisymmetric. Observations do not give us
the inclination of the rotation axis, so this is a free parameter.

As result, our task is to construct an equilibrium *N*-body model
with the given projected surface distribution of particles and the
given line-of-sight kinematics. The total mass of the model is un-
known (*M*/L is unknown), but should be found somehow. The model
should be axisymmetric with a given axis of rotation. The last con-
dition is of course optional, but simplifies the modelling.

In RAS09, we presented an iterative method for constructing
equilibrium *N*-body models with given properties. The idea of
the iterative method is very simple. It relies on constrained, or
guided, evolution. We simply evolve the system while constraining
the desired system properties (see RAS09). Conceptually the same
method, with only relatively minor modifications, can be applied to
our present task.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** The scheme of the iterative method which we use for constructing the *N*-body model of an elliptical galaxy from the observed data.</figcaption></figure>

The scheme of the modified iterative method which we use in this
work is outlined schematically in Fig. 2. We will first briefly describe
the whole method and then each part of the algorithm in detail. We
start by building an initial *N*-body model which has a rigid halo
and the desired projected surface distribution of the particles. The
distribution of particles along the line-of-sight can be arbitrary. The
velocities of the particles and the total mass of the system (*M*/L) can
also be arbitrary. This model will be the starting point for the iterative
procedure. At the start of each iteration we calculate the evolution
of the system over a short time. Then we put the model through
a specific procedure which adjusts the total mass of the system
(see Section 3.4) every second iteration. We then, if necessary,
impose the condition of axisymmetry about the given axis. We do
this in the usual way, i.e. by randomizing the particle azimuthal
angles. Next, we fix the surface distribution of the particles (see
Section 3.2). Finally, we fix the line-of-sight kinematics to agree
with the observations (see Section 3.3). We repeat this iteration
cycle a number of times until the velocity distribution and the total
mass of the system stop changing. The model at this stage is already
in equilibrium, or very close to it. To be sure of this, we calculate
the evolution of the system over a long time-scale (3 Gyr in this
article) after which we obtain the final model in stable equilibrium.
Let us stress that we consider this long time-scale evolution only as
a part of our algorithm for constructing equilibrium models.

In the axisymmetric case, we assume that the galaxy rotates about
the axis which lies in the Z-Y plane. This rotation axis is defined by
its angle $\alpha$ with the Z-axis. If $\alpha = 0$ then the rotation axis coincides
with the Z-axis, so the galaxy is edge-on.

Having described the general outline of our method, we will now
describe each individual step.

**3.2 Fixing the projected surface distribution of particles**

Here, we will describe how we fix the projected surface distribution
of particles. We do this using a method very similar to that described
in RAS09 for fixing the particle mass distribution. The idea of that
algorithm is as follows. We start with an *N*-body system, which is
the result of a short evolution (Fig. 2) and to which we will refer
to as the ‘old’ system. We need to fix the mass distribution in this
model according to given constraints. We create a ‘new’ *N*-body
system with the desired mass distribution, and we ‘transfer’ the
velocity distribution from the ‘old’ to the ‘new’ model. The basic
idea of our velocity transfer algorithm is very simple. We assign to
the new-model particles the velocities of those particles from the old model that are nearest to the ones in the new model. This algorithm is described in detail in RAS09 (Section 2.4). We note that this algorithm has a free parameter $n_{nb}$ – the number of neighbours.

It is very easy to modify this algorithm to fix only the projected surface distribution of particles. We need to construct a 'new' N-body system with a given surface distribution of particles. In the new model only the x and z coordinates of the particles are defined. The y-coordinates, the velocities and the mass of the particles should be carried over from the old model. The total mass of the new system is set equal to the total mass of the old system. We note again that in our models all particles have the same mass. In the algorithm described in RAS09 we 'transfer' from the old to the new model only the velocities. In the present case, we need to 'transfer' also the y coordinates of the particles. When we search for the nearest particle we need to do so in the two-dimensional space X-Z. The y-coordinate should not be taken into account because we copy it from the old model particle to the new model particle together with the velocities.

### 3.3 Fixing the line-of-sight kinematics

Let us first describe an algorithm for fixing the line-of-sight mean velocity for the case where we do not assume any symmetry in the system. Our task is as follows. We have an N-body system and some two-dimensional area on the sky plane where we need to fix the line-of-sight mean velocity to the observed value. In our models the sky plane coincides with the XZ plane, so the line-of-sight velocity is the velocity along the Y-axis. The given area is a two-dimensional area in the XZ plane, so when we search the particles which belong to the given area we do it regardless of y-coordinates of the particles. We denote by $\bar{v}_y$ the desired value of the line-of-sight mean velocity in the given area and by $\bar{v}'_y$ the mean value of the y velocities of all particles in the area. We need to change somewhat the particle velocities so that $\bar{v}'_y$ becomes equal to $\bar{v}_y$. This is achieved by setting the new y velocity component of particle i in the given area to

$$v_{yi} = v'_{yi} + (\bar{v}_y - \bar{v}'_y), \quad (2)$$

where $v'_{yi}$ is the current value of the ith particle and $v_{yi}$ is the corrected ith particle y velocity.

The algorithm for fixing the line-of-sight velocity dispersion is very similar. Let us denote by $\sigma_y$ the desired value of the line-of-sight velocity dispersion in the area under consideration. In the given area we calculate the current value of the line-of-sight velocity dispersion $\sigma'_y$ and the current value of the line-of-sight mean velocity $\bar{v}'_y$. The new y velocity component of particle i in the given area is set as

$$v_{yi} = (v'_{yi} - \bar{v}'_y) \frac{\sigma_y}{\sigma'_y} + \bar{v}'_y. \quad (3)$$

Let us now consider a galaxy with the following symmetries (as described in Section 2.2). The galaxy is axisymmetric, with a symmetry axis which lies in the YZ plane, and also has a reflection symmetry with a plane of symmetry perpendicular to the axis of symmetry and containing the centre of coordinates (0, 0, 0). In this case, if we know the line-of-sight velocity distribution in the single quadrant then we know it for the whole system (see Section 2.2 for details). So, if we assume such a symmetry then all observed kinematical data can be 'reduced' to the first quadrant ($x > 0, z > 0$).

In the case of such a symmetry, the task of fixing the line-of-sight kinematics is as follows. We have some two-dimensional area

in the first quadrant of the sky plane for which we have the line-of-sight mean velocity and/or the line-of-sight velocity dispersion. We need to fix these kinematic parameters in the given area taking into account the discussed symmetries. We invert the sign of $v_y$ for each particle with $x < 0$ and then we flip all particles to the first quadrant, i.e. we set $x = |x|$ and $z = |z|$ for each particles. Now we can apply the algorithms for fixing the line-of-sight kinematic parameters which we described in the first part of this section. We then flip all particles back to their original positions and invert the sign of $v_y$ for each particle with $x < 0$.

### 3.4 Adjusting the total stellar mass

Since we have no a priori knowledge of the total stellar mass in the galaxy, we have to construct our models assuming that the stellar $M/L$ is unknown. In our case this implies that the total mass of the system is unknown. This of course raises questions. Does an equilibrium model with properties in agreement with observations exist for a unique value of the mass? Or for a range of values? For example, for a given equilibrium, spherical, isotropic model with known projected surface distribution of particles and line-of-sight velocity dispersion in the centre of the model, the total mass of the model would be uniquely determined. But in our case it is not so obvious. Moreover, as we will show, models with different inclinations of the rotation axis have slightly different total masses (see Table 2, for example models AL0 and AL45). So if we do not fix the inclination of the rotation axis then the total mass may well not be unique.

Nevertheless, we need to find some value of the total mass for which we can construct the equilibrium model. The straightforward way to solve this problem is to construct a series of models with different total masses and choose the one which, according to some definition, is closest to equilibrium, or to find ranges of values for which the resulting equilibrium model is in agreement with the observational constraints. This, however, would be excessively time consuming, due to the large number of models that need to be constructed.

We will, therefore, use a different algorithm, in which we adjust the total mass during the iterative process (Fig. 2). As we described

<table><caption>**Table 2.** Stellar $M/L$ and relative halo mass for the constructed models. The first column gives the name of the model, the second and third columns the stellar mass-to-light ratios in V and B bands, $M/L_V$ and $M/L_B$, respectively. The fourth and fifth columns give the ratio of the halo mass, $M_h(x)$, to the stellar mass, $M_*(x)$, both calculated within a sphere of radius $x = R_e$ and $x = 5R_e$, respectively.</caption><thead><tr><th>Model</th><th><i>M</i>/<i>L</i><sub>V</sub></th><th><i>M</i>/<i>L</i><sub>B</sub></th><th><i>M</i><sub>h</sub>(<i>R</i><sub>e</sub>)<br><i>M</i><sub>*</sub>(<i>R</i><sub>e</sub>)</th><th><i>M</i><sub>h</sub>(5<i>R</i><sub>e</sub>)<br><i>M</i><sub>*</sub>(5<i>R</i><sub>e</sub>)</th></tr></thead><tbody><tr><td>AN0</td><td>4.16</td><td>4.82</td><td>0</td><td>0</td></tr><tr><td>AL0</td><td>3.55</td><td>4.11</td><td>0.13</td><td>0.96</td></tr><tr><td>AH0</td><td>2.76</td><td>3.22</td><td>0.36</td><td>2.52</td></tr><tr><td>AN45</td><td>4.23</td><td>4.90</td><td>0</td><td>0</td></tr><tr><td>AL45</td><td>3.73</td><td>4.32</td><td>0.12</td><td>0.90</td></tr><tr><td>AH45</td><td>3.15</td><td>3.65</td><td>0.30</td><td>2.18</td></tr><tr><td>BN0</td><td>4.21</td><td>4.89</td><td>0</td><td>0</td></tr><tr><td>BL0</td><td>3.55</td><td>4.12</td><td>0.13</td><td>0.96</td></tr><tr><td>BH0</td><td>2.89</td><td>3.35</td><td>0.35</td><td>2.41</td></tr><tr><td>BN45</td><td>4.27</td><td>4.95</td><td>0</td><td>0</td></tr><tr><td>BL45</td><td>3.81</td><td>4.41</td><td>0.11</td><td>0.88</td></tr><tr><td>BH45</td><td>3.19</td><td>3.69</td><td>0.30</td><td>2.16</td></tr></tbody></table>
earlier, in each iterative step we let the system evolve on a short time-scale (Fig. 2). We calculate in the beginning and in the end of this evolution the velocity dispersion along the line-of-sight in some given part of the system, which we denote by $σ_1$ and $σ_2$, respectively. In all experiments described here this given part was a sphere with radius equal to 10 kpc. We note that the value of $σ_1$ is partly defined by the given line-of-sight velocity dispersion which we fix on the previous stage of the iteration (see Fig. 2). On the other hand the value of $σ_2$ is influenced by the total mass of the model. So, if we choose a value of total mass which is not appropriate and we do not change it during the iterative procedure, then $σ_1$ and $σ_2$ will be different after any number of iterations. We want to construct the equilibrium model so that the values of $σ_1$ and $σ_2$ are close, so we adjust the total mass in the system by multiplying all velocities by a factor of $σ_1/σ_2$ and the masses of all the particles by a factor of $σ_1^2/σ_2^2$. By trial and error we found that the iterations converge faster if we adjust the total mass not at every iteration but only every second iteration.

For a model without dark halo this means that, if the system at the end of the evolution was in equilibrium, then the rescaled system will also be in equilibrium, but will have the line-of-sight velocity dispersion in the selected part of the model equal to $σ_1$. For models with halo this explanation is of course not valid, but this is not a serious drawback since there will be further iterations to bring the system to equilibrium. We thus used this algorithm for constructing all our models, including models with dark halo, and it worked well in all cases. That is, all these models constructed by means of the iterative method with this mass-adjusting algorithm were in equilibrium. But, we note that it is possible that for models with a very massive halo dominating the central part of the galaxy, this algorithm may not work and the iterations would not converge. In such cases, one would have to resort to the straightforward but very time-consuming algorithm described above.

# 4 MODELS

## 4.1 Description of the models

We constructed models with three types of halo. The first type is models without halo. The second and third types are rigid NFW haloes (Navarro, Frenk & White 1996, 1997) with density profile

$$ \rho(r) = \frac{\rho_s}{(r/r_s)(1+r/r_s)^2}, \quad (4) $$

where $ρ_s$ and $r_s$ are the characteristic density and scale radius of the halo. The second type was found in N09 to be the best-fitting NFW halo for this galaxy (see section 4.2.5 in N09). This is a relatively light halo with parameters $ρ_s = 0.0019 M_\odot pc^{-2}$ and $r_s = 32$ kpc, i.e. a concentration parameter $c_{vir} \approx 8$ and a virial mass $M_{vir} \approx 10^{12} M_\odot$ (see N09). The third type is a relatively massive halo with parameters $ρ_s = 0.00522 M_\odot pc^{-2}$ and $r_s = 26.5$ kpc. N09 found that a model with such massive halo is a relatively poor fit of the observational data (see their section 4.2.5). Since we use a totally different approach, we include this model to test whether this conclusion is method-dependent or not. This halo has $c_{vir} = 12.3$ and $M_{vir} = 2 \times 10^{12} M_\odot$.

Another parameter of our models is the angle $α$ defining the inclination of the rotation axis (see Section 3). We construct models for two values of the angle $α$, namely 0° and 45°.

As we discussed in Section 2.2.3 the azimuthal variation of the PNe velocity dispersion presents an interesting feature in the outer parts. Namely, the velocity dispersions in areas close to the principal

axes (areas A1 and A5) are considerably smaller than in intermediate areas (areas A2, A3 and A5). It is not clear whether such a feature can be reproduced by an equilibrium axisymmetric system. So we construct two sets of models. In the first ones, which we denote as ‘A’, we do not try to model this feature. Thus, we ask the iteration method to fit the following quantities:

(i) the projected surface distribution of the particles (see Section 2.1).

(ii) the radial profile of the mean velocity along the major axis obtained from long-slit spectroscopy.

(iii) the radial profiles of the velocity dispersion along the major and minor axes obtained from long-slit spectroscopy (see Section 2.2.2).

(iv) the mean velocity in area A1 obtained from PNe kinematics. Note that we do not fix the mean velocity in other areas (see Section 2.2.3).

(v) the velocity dispersion in the two PNe areas A1 and A2, which are close to the major and minor axes, respectively.

For the second set of models, which we denote by ‘B’, we take into account all available information. We therefore try to fit the velocity dispersion in all five PNe areas separately. This introduces a considerable extra difficulty because of the peculiar PNe velocity distribution. As already mentioned, such a fit may, of course, not be possible with models such as ours, i.e. with models that are both axisymmetric and in equilibrium. It is, nevertheless, important to try for at least two reasons. First it is useful in all cases to try and fit all available data, since only that can tell us how far off the attempted fit is from reality. Secondly, since we allow the model total freedom regarding the velocity anisotropy, a fit may be found, in which case an interesting effect of anisotropy would be revealed.

We will, therefore, discuss in total 12 models, which we will denote as follows. The first symbol in the name of the model denotes the set (A or B) and the second denotes the halo. Here ‘N’ is for models without halo, ‘L’ is for models with a relatively light NFW halo ($ρ_s = 0.0019 M_\odot pc^{-2}$, $r_s = 32$ kpc) and ‘H’ is for models with a relatively massive NFW halo ($ρ_s = 0.00522 M_\odot pc^{-2}$, $r_s = 26.5$ kpc). The last number in the model name denotes the angle $α$ of the rotation axis inclination.

Let us now discuss in some detail how we construct these models. There is considerable freedom in choosing the initial model, from which the iterative search will start. Our basic initial model is a model with the given surface distribution of particles (generated using the rejection method), and with zero velocities. The y coordinates of the particles were chosen as random numbers from the interval $[-1, 1]$ and the initial mass was chosen according to $M/L_v = 1$. These choices are of course totally ad hoc, but this does not matter since other values, although starting the iterative procedure (see Fig. 2) from a different initial model, lead to essentially the same final model. In practice, there is only one significant property of this initial model, namely whether it is rotating (or not). Indeed, if we start the iterative search from a rotating, rather than a non-rotating initial model, we will end up with a model that is different, albeit not in all properties. A rotating initial model can be constructed as follows. We first take the basic non-rotating initial model and put it in the iterative procedure (see Fig. 2). After a relatively small number of iterations, the model becomes close to steady state. We then choose the axis of rotation and set, for all particles, the azimuthal velocity with respect to the chosen axis equal to the circular velocity. We have thus obtained a rotating initial model.

Models with $α = 45$ constructed from rotating and non-rotating initial models are practically identical. Non-inclined models
(α = 0), however, are slightly different. In particular, there is a
difference in rotation in areas far from major axis. We note that we
fix the mean line-of-sight velocity only in areas close to major axis
(long-slit data along major axis and mean velocity in A1 area). It is
not surprising that models constructed from initially rotating models
rotate slightly faster in areas far from major axis. But this difference
is not sufficiently significant to warrant further discussion.

Here, we will discuss models constructed from the rotating initial
model. Each model was constructed as follows. At first we construct
a basic non-rotating initial model with $N = 300\,000$ particles. We
put this model into the iterative procedure and make 100 iterations
with relatively low precision. This is possible because each itera-
tion is very short and errors do not accumulate (RAS09). We use the
fast N-body code GYRFALCON (Dehnen 2000a, 2002) with an inte-
gration step and a softening length equal to $dt = 1/2^{12}$ Gyr and $\epsilon =
0.05$ kpc, respectively.<fnref n="3" /> The tolerance parameter for GYRFALCON was
set to $\theta_t = 0.9$ and the duration of each iteration to $t_i = 0.05$ Gyr.
Using this constructed model, we create the rotating initial model
with $N = 500\,000$ particles. Again, we put this model into the iter-
ative procedure and make 500 iterations. The integration step and
the softening length were taken as $dt = 1/2^{14}$ Gyr and $\epsilon = 0.02$ kpc,
respectively. The duration of each iteration and the tolerance pa-
rameter for GYRFALCON were chosen as in the previous stage. The
final stage of our procedure is a free evolution over a time-scale of
3 Gyr (see Section 2), with parameters $dt = 1/2^{15}$ Gyr, $\epsilon = 0.02$ kpc
and $\theta_t = 0.6$. We would also like to mention that before fixing the
parameters to these values, we made a number of tests, such as
increasing the number of particles up to tenfold, and did not find
any significant improvements in the fits.

The stellar mass-to-light ratio and relative mass of the dark
halo inside one and five effective radii, respectively, are shown in
Table 2. According to N09, the value of the effective radius $R_e =
48.2$ arcsec $\approx 3.69$ kpc. We use photometry in the V band, so for us
it is more straightforward to calculate the mass-to-light fraction in
V-band $M/L_V$. To compare our results with N09 we also calculate
mass-to-light fraction in B-band $M/L_B$. We assume that ($B-V$) is
0.65 for the Sun and 0.81 for NGC 4494 (see N09). We note that, if
we change the adopted distance, the mass-to-light ratio will also be
changed (see Section 2). For all models, there is considerably less
halo mass than stellar mass within $1R_e$. For models with a ‘light’
halo, the mass of the dark matter inside $5R_e$ is approximately equal
to the mass of the stellar component, while for model with a ‘heavy’
halo the mass of dark matter inside $5R_e$ is more than twice that of
the stellar component.

**4.2 Discussion of the first set of models**

4.2.1 Basic comparison

When constructing our first set of models, models A, we use only
velocity information from the vicinity of the principal axes. In par-
ticular, we take into account the velocity dispersion in PNe areas
A1 and A5 and ignore observations in areas A2, A3, A4. In so
doing we want to demonstrate that our method can be used to con-
struct galactic models following given observational constraints,
both photometric and kinematic.

Let us first discuss models AN0, AL0 and AH0, which have an
inclination angle $\alpha = 0$, i.e. their rotation axis is parallel to the sky
plane. Our results for these models are summarized in Table 2 and
Figs 3 and 4. In order to increase the resolution and reduce the noise
in our figures we use a trick described by Athanassoula (2005b).
We stack 10 snapshots closely spaced in time, with $\Delta t = 10$ Myr,
and we calculate all values for this combined snapshot. This allows
us to reduce the noise in our plots quite significantly.

Since we include rotation, we compare models to observations
separately for mean velocities and for velocity dispersions, and do
not fold these two quantities into a single root mean square velocity
$v_{rms} = \sqrt{v^2 + \sigma^2}$ (N09).

In general, there is a good agreement between the models and
the observational data which were used by the iterative method for
constructing them (Fig. 3). We get an excellent fit of the density
profile (panels a and b). We also get good agreement with the mean
velocities on the major axis except for radii within the centremost
region (panels d and e) and with the velocity dispersion on major and
minor axes (except for the bump at radii between 60 and 90 arcsec;
see panels f and g). All these agreements hold for all three models,
i.e. models AN0, AL0 and AH0. Also all our models fit very well
mean velocity of PNe in area A1 (panel h). Velocity dispersion in
PNe areas will be discussed later in this section. Let us now discuss
in more detail the parts where the models fail to reproduce the
observations.

(i) As many other ellipticals, NGC 4494 has a kinematically de-
coupled core. This is clear not only from the mean velocity, but also
from the velocity dispersion profiles along the major and minor
axes. Such features are believed to be out of equilibrium, e.g. due
to a merger with a small companion that reached the Galactic Cen-
tre by dynamical friction. Furthermore, the material in this region
may not be axisymmetrically distributed and/or have a different ori-
entation from the remaining galaxy. Arguments in favour of these
possibilities are that the density distribution near the centre does not
follow the Sérsic law and that the profile of the mean velocity does
not have the shape expected for axisymmetric mass distributions in
equilibrium. Departures from axisymmetry and from equilibrium
are beyond the scope of this paper. It is thus normal that neither
N09 – with a spherically symmetric non-rotating equilibrium model
– nor we – with an axisymmetric, rotating or non-rotating equilib-
rium model – reproduce the structure in the innermost part. Models
including kinematically detached cores are beyond the scope of both
studies.

(ii) The velocity dispersion profiles have a bump along both the
major and minor axes in the area between 60 and 90 arcsec (panels
f and g of Fig. 3). N09 considered the long-slit data as kinematical
constraints only up to 60 arcsec, where they considered them as
more accurate. But their plots clearly show that their models have
no such bump in the 60 to 90 arcsec region. We kept all the long-slit
data as constraints, but none of our models reproduced this bump. It
is possible that in NGC 4494 this bump is transient, presumably due
to some collision event and thus cannot be reproduced by the Jeans
method, or by our iterative method, which can only find equilibrium
solutions. This, together with the previous discussion on the centre-
most area, argues that some regions of NGC 4494 are not in exact
steady state, as often observed and as could be expected in the
framework of the ΛCDM paradigm.

(iii) The ellipticity of model AN0 agrees on average with the
mean observed value (panel c). However, models AL0 and AH0
have bigger ellipticity values than observed in their outermost parts.
This is partly due to the fact that these two models have slightly

<footnote n="3">³ We use a system of units where the unit of length is $u_1 = 1$ kpc, the unit of velocity is $u_v = 1 \text{ km s}^{-1}$, the unit of mass is $u_m = 10^{10} \text{ M}_\odot$ and consequently the unit of time is $u_t \approx 0.98$ Gyr. For simplicity, all time values in this papers are presented with the assumption that $u_t = 1$ Gyr.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Comparison of models AN0, AL0 and AH0 with the observational data. (a) and (b) show the dependence of $n_{xz}$ on $R_{xz}$, where $n_{xz}$ is the number of particles in concentric cylindrical shells and $R_{xz} = \sqrt{x^2 + z^2}$. The thick solid line shows the profile calculated for a model with the input data prepared as described in Section 2.1. Panel (c) shows the ellipticity profiles of the models projected on the sky plane (XZ plane), calculated with the IRAF task ELLIPSE. Panels (d) and (e) show the profile of the mean velocity along the major axis and panels (f) and (g) the profiles of the velocity dispersion along the major and minor axes, respectively. In (d) and (f), together with the long-slit data, we show PNe data in area A1 which is close to the major axis (see Fig. 1). The value calculated for this area is assigned to radius $(R_{xz}^{(min)} + R_{xz}^{(max)})/2$ (Table 1). In panel (g), together with the long-slit data, we show the PNe data in area A5 which is close to the minor axis. In panels (d), (f) and (g) the observed data are binned as described in Section 2.2.2. Panel (h) shows the mean line-of-sight velocity calculated for the five PNe areas (see Section 2.2.3). Note that only the mean velocity in area A1 is used for constructing the models we use. Panel (i) shows the line-of-sight velocity dispersion calculated for the five PNe areas. Only dispersions in areas A1 and A5 are used for the construction of A models. The radial profiles for models AN0, AL0 and AH0 are given by solid, dashed and dotted lines, respectively. The open squares show the long-slit data, not binned for panel (e), and binned as discussed in Section 2.2.2 for the remaining panels. The filled circles show the PNe data, which were used as observational constraints for these models, and the X signs the remaining PNe data.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Radial profiles of velocity anisotropy for models AN0, AL0 and AH0. The velocity anisotropy parameter is calculated as $\beta = 1 - \frac{\sigma_r^2}{\sigma_t^2}$, where $\sigma_r$ and $\sigma_t$ are the velocity dispersion in the radial and the $\varphi$ direction, respectively, in a spherical coordinate system.</figcaption></figure>

boxy isophotes at their periphery, which we believe to be connected to the very high velocity anisotropy in their periphery (Fig. 4), i.e. to be linked to the radial orbit instability. As already mentioned, our algorithm for model construction includes a self-consistent N-body evolution on a relatively long time-scale (Fig. 2), so that our models can ‘feel’ such instabilities. Let us, however, underline that this effect is relatively small and the isophote boxyness is rather subtle. As a result the difference with the average ellipticity is of the order of only 0.05, i.e. much smaller than the corresponding difference for spherical models, such as those of N09.

We can conclude that our models reproduce well the projected surface density, the mean apparent ellipticity, the mean velocities and their dispersions (the latter though not over the full radial and angular extent of the galaxy). The corresponding fits are no worse than the N09 models. The discrepancies between our models and observations occur mainly in regions whose data were not considered by N09.

Fig. 4 shows the radial profile of the velocity anisotropy for our $\alpha = 0$ models and clearly illustrates the degeneracy between mass and anisotropy. These three models have different halo masses, but the same stellar density profile and observed kinematics (Fig. 3) and this is possible because they have different velocity anisotropies. As was expected, radial anisotropy is increasing with halo mass (see Fig. 4).

The inclined models AN45, AL45 and AH45 are considerably different from the previously discussed AN0, AL0 and AH0 models, as they have an intrinsic ellipticity of 0.35. The quality of the fits, however, is very similar (Fig. 5).

Seen the size of the error bars and the small effect of the halo mass on the density and velocity radial profiles, it is not possible from the above comparisons alone to rule out any of the models, or even to set a strong preference to one rather than another. The only exception is the ellipticity profiles, which are considerably better fitted by the haloless models AN0 and AN45 (Figs 3 and 5). On the contrary, all models with dark halo have in their outer parts slightly boxy

isophotes (presumably due to their higher velocity anisotropy) and a higher than mean observed ellipticity. In our modelling, however, we use only the mean value of the ellipticity and not the full profile. Also we do not use high-order isophote shape parameters, so we enforce our models to have precise elliptical isophotes. NGC 4494 has a negative $a_4$ isophote-shaped parameter on periphery, i.e. boxy isophotes (see appendix A in N09). So we cannot rule out models with dark haloes relying only on the fact that they have slightly boxy isophotes in their outer parts.

Models without dark halo fit perfectly the velocity dispersion of the PNe that have been used as observational constraints, i.e. in regions close to the major and minor axes [see panel (i) of Figs 3 and 5]. Model AL45, however, with a relatively light dark halo, also fits the observational data along the major and minor axes rather well, while AL0 fits only slightly worst [see panels (i) of Figs 3 and 5]. Seen the error bars, even model AH45 with a relatively massive dark halo cannot be ruled out.

### 4.2.2 Higher-order moments of the velocity distribution

When we constructed the models, we did not use the higher-order moments of the velocity distribution as observational constraints. It is thus interesting to compare the third- and fourth-order moments of the velocity distribution of our models to those of the observational data. Moreover, such moments may help break the mass–anisotropy degeneracy and were used by N09 to argue for the need of a low-mass halo. At the periphery of the galaxy, we have kinematical information only from the PNe. We thus calculate the skewness and the kurtosis in the five PNe areas (Fig. 1) and compare them to the corresponding model values, calculated for our models as described in Appendix B. It should, however, be noted that the number of PNe is very small and that this affects the higher-order moments more. Thus, even assuming a Gaussian velocity distribution, the uncertainties for skewness and kurtosis are very big, while without this assumption the uncertainties are formally infinite (Appendix B).

One of the advantages of our models is that we can easily calculate any parameters for them. This will allow us now to calculate for our models exactly the same high-order moments of the velocity distribution as for observations. We can thus calculate two Gauss–Hermite moments ($h_3$ and $h_4$) along the major and minor axes of the galaxy, as described in Appendix A, and compare them with the corresponding profiles from long-slit spectroscopy. The results are given in Figs 6 and 7 for $\alpha = 0$ and $\alpha = 45$, respectively. These figures clearly show which of the higher-order moments depend on the halo mass and in which way.

The most interesting of the Gauss–Hermit moment profiles is probably the $h_3$ profile along the major axis. This is visibly different for models with different halo mass [see panels (a) of Figs 6 and 7], the most massive halo having the highest $h_3$ values and AN0 the lowest. The best fit seems to be for the intermediate-halo mass. To establish this we calculated the $\chi^2$, excluding the region of the kinematically decoupled core, and normalized it by the number of data points. We find that for $\alpha = 0$, the values are 3.9, 0.4 and 5.7 for models with no halo, light halo and heavy halo, respectively. The corresponding numbers for $\alpha = 45$ are 2.4, 0.8 and 1.8, respectively. These numbers show a preference for the models with light halo, and argue, albeit weakly, for a preference for $\alpha = 0$. The profile of $h_4$ along the major axis also shows some dependence with halo mass, but much less so than the corresponding $h_3$ profile [panel (d) of Figs 6 and 7]. Profiles along the minor axis for both Gauss–Hermite
<figure><img src="image_6.png" /><figcaption>Figure 5. As in Fig. 3, but for models AN45, AL45 and AH45.</figcaption></figure>

moments show no clear dependence on halo mass [panels (b) and (e) of Figs 6 and 7], except for the $h_4$ profiles in the outermost parts, where, however, there are no long-slit data. In the central part of the model, i.e. inside ~30 arcsec, all our models have $h_4$ values which are clearly less than those of the observations for both axes [Panel (a) of Figs 6 and 7] and this is true also for the $h_3$ major axis profiles. This could again be linked to the kinematically decoupled core.

The kurtosis values calculated in each of the PNe areas for all 'A' models are, for practical purposes, almost the same [panels (f) in Figs 6 and 7], the differences being much smaller than the observational errors. The situation with the skewness is the same

[panels (c) in Figs 6 and 7]. We therefore do not believe that the comparisons in panels (c) and (f) of Figs 6 and 7 can be used to distinguish between models of different mass.

To summarize, the $h_3$ profiles along the major axis provide some arguments in favour of models with a light halo (particularly AL0 and AL45), in good agreement with what was found by N09.

## 4.3 Discussion of the second set of models

As discussed in Section 2.2.3, the velocity distribution of the PNe has an interesting feature, namely that the velocity dispersion near
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison of higher moments of the velocity distribution of models AN0, AL0 and AH0 with the observational data. Panels (a) and (b) show profiles of the $h_3$ Gauss-Hermite moment along the major and minor axes, respectively (calculated as described in Appendix A). Panel (c) shows an estimate of the skewness $G_1$ calculated for the five PNe areas (calculated as described in Appendix B). (d) and (e) show profiles of the $h_4$ Gauss-Hermite moment along the major and minor axes, respectively. (f) shows an estimate of the kurtosis $G_2$ calculated for five PNe areas. These estimations of the model skewness and the kurtosis are 'reduced' to a sample size equal to the number of PNe in the corresponding areas (see Appendix B). The error bars in panels (c) and (f) were calculated assuming a Gaussian distribution function (see Appendix B). The symbols and line styles are as in Fig. 3.</figcaption></figure>

the principal axes (regions A1 and A5) has much smaller values than the velocity dispersions in the intermediate areas (regions A2, A3 and A4). When building the models presented in the previous section, we did not use this as an observational constraint and ignored observations in intermediate areas. As results all our 'A' models have almost the same dispersion in all PNe areas [see panels (i) in Figs 3 and 5], so they obviously fail to reproduce observational data in intermediate areas. Here, we consider more realistic models ('B' models) for the construction of which we used all available observational data, to see whether the tangential variation of the velocity dispersion can be reproduced by rotating axisymmetric models.

In general, both sets of models are very similar, with differences only at the periphery of the galaxy. Furthermore, the quality of the fit to the observation data is also similar. Thus, to save space, we show only a very restricted set of figures. Fig. 8 shows the velocity dispersion calculated in the five PNe areas for our B models.

We find that only non-inclined models with a dark halo (i.e. models BL0 and BH0) reproduce this feature, and even only partially. Namely, the velocity dispersions near the principal axes (areas A1 and A5) are too high, while at the middle (region A3) it is marginally too low. Thus, we were able to qualitatively reproduce the form of the azimuthal variation, but not quantitatively, since the model am-

plitude is lower than necessary to fit the observations. For these models, we cannot make specific comparisons with N09, because the latter did not try and reproduce this feature.

We can thus conclude that, although our models reproduce most of the observed features well, they fail to fully reproduce the velocity distribution of PNe. We will discuss possible explanations in the next section.

## 5 CONCLUSIONS

We used our iterative method (RAS09) to construct dynamical models of NGC 4494. One of the advantages of our method is that the models are produced in the form of *N*-body snapshots with equal-mass particles, which can be directly used in *N*-body simulations. This, for example, allows us to directly check whether the constructed model is in steady state and to calculate any quantities, or parameters of these models, which can then be compared to observations. As already discussed in Section 3 (see also Fig. 2), after the main part of the iterative method, we let our models evolve freely on a time-scale of 3 Gyr, and then we calculate the quantities to be compared with the observational data. Since the aim of this evolution is to check whether the models we constructed are indeed in steady state, we used the same halo as when building the model, i.e.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Comparison of higher moments of the velocity distribution of models AN45, AL45 and AH45 with observations. The values shown are the same as in Fig. 6.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of the line-of-sight velocity dispersion calculated for the five PNe areas in the ‘B’ models with observational data. Panel (a) refers to models BN0, BL0 and BH0 and panel (b) to models BN45, BL45 and BH45.</figcaption></figure>

a rigid halo. This evolution clearly showed that all our models are in a steady state. For the sake of completeness, however, we also tried simulations with a live halo and find no considerable differences from the rigid halo simulations, arguing that in these cases there is not much interaction between the dark and the stellar matter.

We used the observational data given by N09, namely surface photometry, stellar kinematics along the major and minor axes as obtained by means of long-slit spectroscopy, and velocities and positions of 255 PNe (see N09). We use PNe data only in the outer part of the galaxy where long-slit data are absent. The surface
photometry gives us the surface distribution of particles but not the total mass of the system because the mass-to-light ratio is unknown. Our algorithm automatically adjusts the total mass of the model (Section 3).

We constructed models with three types of halo, all of which were already discussed in N09. The first type is models without halo. The second type of halo is a relatively light NFW halo, which N09 found to be the best-fitting NFW for NGC 4494. The third type of halo is a relatively heavy NFW halo. We, furthermore, constructed two sets of models: models ‘A’ and models ‘B’. These two sets have only one difference. When constructing models ‘A’ we used as observational constraints only the PNe in areas close to the principal axes, neglecting information on the velocity distribution of the PNe in intermediate areas. On the contrary, in order to construct models ‘B’ we used all available information.

One important goal in making the ‘A’ models was to demonstrate the ability of our method to construct equilibrium models with given observation constraints. This was fully achieved. In general, there is a good agreement between our models and the observational data concerning the projected surface density, the mean apparent ellipticity, the mean velocities and their dispersions (except for specific radial ranges, where some non-equilibrium or non-axisymmetric structure could be present in NGC 4494). We showed that our models reproduce observation data not worse than N09 models. But our models have the added advantage that they are live, non-spherical, rotating *N*-body systems.

A further goal was to see whether it was possible to set constraints on the dark halo mass with our models. If we do not take into account the high-order moments of the velocity distribution, then the best models are the models without dark halo. This, however, is only a slight preference and models with a light halo and even models with a heavy halo cannot be rejected. It is thus necessary to try higher-order moments. We found that the major axis profile of the third-order moments shows the strongest dependence on the halo mass and that the best fit was for the models with the light halo. This is in agreement with what N09 found for their models, though with a different technique.

N09 found that a heavy halo gives worst fits to the data, without referring to higher-order moments, i.e. from the root mean square velocity profile, which is a combination of the mean velocity and the dispersion. This, however, is not the case for our models. For example, model AH45, which has a heavy halo, fits all low-order moments of the velocity distribution very well. These two results taken together show that there are no simple spherical models fitting the lower order models, but that there can well be more general models that do. Thus, our model which is axisymmetric, rather than spherical, rotating and inclined with respect to the line of sight has no difficulty with the lower-order moments. When we place the bar higher, i.e. when we ask our models to fit also the third and fourth-order moments, we find that the models with heavy haloes do worse than the models with lower mass haloes, as we described in Section 4.2.2. This does not necessarily mean that there are no heavy halo models that fit the third- and fourth-order moments. It can simply mean that it is necessary to consider a yet more general model, e.g. a triaxial one, in order to get such a fit. In other words, one has to be careful not to extrapolate a given result further than the class of models for which it was derived. This, of course, makes it very difficult to set any strong constraints to the halo mass in NGC 4494, and to ellipticals in general.

N09 found that some halo is required by comparing the kurtosis of their model with that of the PNe velocity distribution at the outermost parts of NGC 4494. Again, this is not the case for our

models. It is the $h_3$ radial profile along the major axis that proved in our case to be more sensitive to the halo mass, and thus allowed us to set some preference for the light halo model, in agreement with what was found by N09 for the fourth-order moment.

When constructing and analysing the ‘A’ models we ignore the PNe data at angles intermediate between the major and the minor axes. Since this is rather ad hoc, we built the ‘B’ models, in which all PNe data were used as observational constraints. By doing so, we wanted to test whether we could make models which reproduced the azimuthal variation of the velocity dispersion, i.e. which have relatively low values near the principal axes and considerably higher ones in between. Although we were able to build models that qualitatively reproduce this feature, we were unable to reproduce it quantitatively, i.e. our models always display a lower amplitude of this variation than the observations. There are several possible explanation to this.

It should first be noted that the observational data at the periphery of NGC 4494 is very sparse. We have only around 15 PNe in each of our PNe areas (see Table 1). So even the formal uncertainties of the PNe dispersions in these areas are rather big. Moreover, there could still be some contamination in the PNe sample (see section 2.2 in N09), which would further increase the uncertainties. As a result, there is still a considerable probability that the observations do not contradict, e.g. model BL0 (see Fig. 8a, error bars are one $\sigma$).

There is, of course, always the possibility that a model which fully reproduces the observations exists, but that our method failed to construct it. We made very extensive searches, and we do not believe that this possibility is likely, but, strictly speaking, we cannot exclude it. Furthermore, the observations point to other, more likely alternatives, which we discuss below.

(i) In our models the halo is spherical. It is likely that models with an axisymmetric or a triaxial halo would reproduce observations better. Note, however, that this alternative cannot help the models without a dark halo.

(ii) Our models are axisymmetric, while, as discussed in Section 2.2.3, the velocity distribution in the galaxy is not fully axisymmetric. It would thus have been better to consider more general types of models, for example, triaxial models or models with even less symmetry. This, however, would introduce further free parameters and is beyond the scope of this paper.

(iii) It is possible that the outer parts of NGC 4494 are not in a steady state. As we already discussed, there is evidence that NGC 4494 is not exactly in a steady state. If there are transient features in the innermost and perhaps in part some of intermediate part of the galaxy (the kinematically decoupled core and between 60 and 90 arcsec), then such features can also be present at its periphery. Moreover, the dynamical time-scale at the periphery is much longer than in the intermediate part, so that any non-equilibrium structures will tend to equilibrium much slower than in regions further in.

(iv) Seen the low number of PNe in the periphery of the galaxy, the uncertainties are very large and therefore the constraints in these regions very loose. It might therefore not be necessary to dismiss a model if the fits in these regions are poor.

Any of the three first alternatives would necessitate a more general model than what we have considered here. Triaxial models could, at least in principle, be built, but the number of free parameters would increase so that more observational constraints would be necessary. However, neither our techniques, nor other techniques discussed so far can build non-equilibrium models. This could be done only by *N*-body simulations, which, however, would have an extraordinarily unwieldy free parameter space.
To summarize, we were able to reach our first goal, i.e. we showed that our iterative method (RAS09) can be used for building models with given observational constraints. Concerning our second goal, i.e. to set strict constraints on the halo mass of NGC 4494, we can only claim a more limited success. Comparing the third-order moments of our model velocity distribution with the observations, we find that it is the light halo model that gives the best fit. Even this, however, does not fit all available data, and it could be that a more general model (e.g. triaxial) is necessary. It is, however, not possible to extend the preference of light haloes to these more general models, without actually making them and analysing their properties. It could thus well be that it is another type of model (such as with no halo, or with a heavy halo) that gives the best fits in such a case. Thus the present set of models, although showing a preference for a light halo, cannot set a very strict constraint to the halo mass of NGC 4494, particularly seen the sparseness of the PNe in the outermost parts.

## ACKNOWLEDGMENTS

We thank Albert Bosma for useful discussions on the data. This work was partially supported by grant ANR-06-BLAN-0172, by the Russian Foundation for Basic Research (grants 08-02-00361-a and 09-02-00968-a) and by a grant from the President of the Russian Federation for support of Leading Scientific Schools (grant NSh-1318.2008.2).

## REFERENCES

Athanassoula E., Vozikis Ch., 1999, in Barnes J. E., Sanders D.B., eds, Proc. IAU Symp. 186, Galaxy Interactions at Low and High Redshift. Kluwer Academic Publishers, Dordrecht, p. 145

Athanassoula E., 2005a, in Szczerba R., Stasińska G., Górny S. K., eds, AIP Conf. Proc. 804, Planetary Nebulae as Astronomical Tools. Melville, New York, p. 333

Athanassoula E., 2005b, MNRAS, 358, 1477

Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn. Princeton Univ. Press, Princeton

Bosma A., 2004, in Ryder S. D., Pisano D. J., Walker M. A., Freeman K. C., eds, Proc. IAU Symp. 220, Dark Matter in Galaxies, Astron. Soc. Pac., San Francisco, p. 39

Coccato L. et al., 2009, MNRAS, 394, 1249

de Lorenzi F., Debattista V. P., Gerhard O., Sambhus N., 2007, MNRAS, 376, 71

de Lorenzi F., Gerhard O., Saglia R. P., Sambhus N., Debattista V. P., Pannella M., Méndez R. H., 2008, MNRAS, 385, 1729

de Lorenzi F. et al., 2009, MNRAS, 395, 76

Dehnen W., 2000a, ApJ, 536, 39

Dehnen W., 2000b, AJ, 119, 800

Dehnen W., 2002, J. Comput. Phys., 179, 27

Dekel A., Stoehr F., Mamon G. A., Cox T. J., Novak G. S., Primack J. R., 2005, Nat, 437, 707

Douglas N. G. et al., 2007, ApJ, 664, 257

Gerhard O. E., 1993, MNRAS, 265, 213

Goudfrooij P., Hansen L., Jorgensen H. E., Norgaard-Nielsen H. U., de Jong T., van den Hoek L. B., 1994, A&AS, 104, 179

Joanes D. N., Gill C. A., 1998, Statistician, 47, 183

Lauer T. R. et al., 2005, AJ, 129, 2138

Machado R. E. G., Athanassoula E., 2010, MNRAS, 406, 2386

Méndez R. H., Riffeser A., Kudritzki R.-P., Matthias M., Freeman K. C., Arnaboldi M., Capaccioli M., Gerhard O. E., 2001, ApJ, 563, 135

Napolitano N. R. et al., 2009, MNRAS, 393, 329 (N09)

Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493

Kenney J. F., Keeping E. S., 1951, Mathematics of Statistics, Pt. 2, 2nd edn. Van Nostrand, Princeton, NJ

Rodionov S. A., Orlov V. V., 2008, MNRAS, 385, 200

Rodionov S. A., Athanassoula E., Sotnikova N. Ya., 2009, MNRAS, 392, 904 (RAS09)

Romanowsky A. J., Douglas N. G., Arnaboldi M., Kuijken K., Merrifield M. R., Napolitano N. R., Capaccioli M., Freeman K. C., 2003, Sci, 301, 1696

Sérsic J. L., 1968, Atlas de galaxias australes. Observatorio Astronomico, Cordoba, Argentina

Syer D., Tremaine S., 1996, MNRAS, 282, 223

van der Marel R. P., Franx M., 1993, ApJ, 407, 525

Weil M. L., Hernquist L., 1994, ApJ, 431, 79

Weil M. L., Hernquist L., 1996, ApJ, 460, 101

## APPENDIX A: CALCULATION OF GAUSS-HERMITE MOMENTS IN THE CASE OF N-BODY SYSTEMS

Let us now describe how we calculate the Gauss-Hermite moments in the case of N-body systems. Both van der Marel & Franx (1993) and Gerhard (1993), independently, discuss the use of Gauss-Hermite moments to measure the deviations of the observed velocity distribution from a Gaussian, but with different normalization of the Hermite polynomials. Here, we use the normalization of Gerhard (1993).

The Hermite polynomials are defined as

$$H_n(x) = (-1)^n e^{x^2} \frac{d^n}{dx^n} (e^{-x^2}). \quad (A1)$$

The sequence of Hermite polynomials satisfies the recursion

$$H_{n+1} = 2x H_n - 2n H_{n-1} \quad (A2)$$

and the first three Hermite polynomials are

$$H_0 = 1, \quad H_1 = 2x, \quad H_2 = 4x^2 - 2. \quad (A3)$$

The set of functions defined as

$$u_n(x) = (2^{n+1}n!\pi)^{-1/2} H_n(x) \exp(-x^2/2) \quad (A4)$$

obey the orthogonality relation

$$\int_{-\infty}^{+\infty} u_n(x)u_m(x)dx = \frac{\delta_m^n}{2\pi^{1/2}}. \quad (A5)$$

Thus this set of functions is a complete orthogonal system.

The Gauss-Hermite moments for some function $l(v)$ are defined as

$$h_n = 2\pi^{1/2} \gamma_h^{-1} \int_{-\infty}^{+\infty} l(v) u_n(w) dv, \quad w \equiv (v - V_h) / \sigma_h, \quad (A6)$$

where $V_h$, $\sigma_h \neq 0$ and $\gamma_h \neq 0$ are free parameters. The function $l(v)$ can be approximately calculated by means of a truncated Gauss-Hermite series

$$l(v) \approx \frac{\gamma_h}{\sigma_h} \sum_{i=0}^{N_h} h_i u_i(w), \quad (A7)$$

where $N_h$ is the number of terms used (see Gerhard 1993; van der Marel & Franx 1993 for more details).

In the case of an N-body system, we need to solve the following problem. We have the set of values $v_1, v_2, \dots, v_N$, and we need to calculate the Gauss-Hermite moments of the distribution function $l(v)$ of these values. For example, if we need to calculate Gauss-Hermite moments for the line-of-sight velocity distribution in some area of an N-body system then $v_i$ is the line-of-sight velocity of
each particle in this area. In this case, the Gauss-Hermite moments for the function $l(v)$ can approximately be calculated as

$$h_n = \frac{2\pi^{1/2}}{\gamma_h N} \sum_{i=1}^{N} u_n(w_i), \quad w_i \equiv (v_i - V_h)/\sigma_h. \qquad (A8)$$

We use this equation as the definition of Gauss-Hermite moments for our discrete case.

To compare Gauss-Hermite moments with the observations, we need to choose the free parameters $V_h$, $\sigma_h$ and $\gamma_h$ as observer would do. These free parameters should be chosen so as to give $h_0 = 1$, $h_1 = h_2 = 0$ (van der Marel & Franx 1993). From (A3), (A4), (A8) for $n = 1$ and the condition $h_1 = 0$, we have

$$h_1 = 0 \Leftrightarrow V_h = \frac{\sum_{i=1}^{N} v_i \exp(-w_i^2/2)}{\sum_{i=1}^{N} \exp(-w_i^2/2)}. \qquad (A9)$$

From (A3), (A4) (A8) for $n = 2$ and the condition $h_2 = 0$, we have

$$h_2 = 0 \Leftrightarrow \sigma_h^2 = 2 \frac{\sum_{i=1}^{N} (v_i - V_h)^2 \exp(-w_i^2/2)}{\sum_{i=1}^{N} \exp(-w_i^2/2)}. \qquad (A10)$$

We note that $w_i \equiv (v_i - V_h)/\sigma_h$, so (A9) and (A10) cannot be solved directly. We solve (A9) and (A10) together, by means of iterations. Initially, we set $V_h$ equal to the mean value of $v_i$, $\sigma_h$ equal to the standard deviation of $v_i$ and $\gamma_h = 1$. A single iteration is as follows

$$V_h^{(\text{new})} = \frac{\sum_{i=1}^{N} v_i \exp(-w_i^2/2)}{\sum_{i=1}^{N} \exp(-w_i^2/2)} \\ \sigma_h^{(\text{new})} = \left( 2 \frac{\sum_{i=1}^{N} (v_i - V_h)^2 \exp(-w_i^2/2)}{\sum_{i=1}^{N} \exp(-w_i^2/2)} \right)^{1/2}. \qquad (A11)$$

After finding the appropriate $V_h$ and $\sigma_h$ values, we can easily calculate the last free parameter as $\gamma_h^{\text{new}} = h_0$.

## APPENDIX B: CALCULATION OF SKEWNESS AND KURTOSIS FOR A SMALL SAMPLE

In the outer parts of the galaxy, we have information about the line-of-sight velocity distribution only from PNe observations. Our task is to compare the velocity distribution of the PNe in some area on the sky plane with the velocity distribution in a constructed N-body model (see Fig. 1). More precisely, our task can be formulated as follows. We have two random samples which we denote here as ‘A’ and ‘B’. Sample ‘A’ with size $n_a$ consists of the line-of-sight velocities of the PNe in the selected area. Sample ‘B’ with size $n_b$ consists of the line-of-sight velocities of model particles in the same area. We need to assess the probability that these two samples were generated from the same distribution function. This can be achieved by comparing moments calculated for these two samples.

Here, we will discuss the comparison of the skewness and kurtosis calculated for these two samples. In our case the number of PNe in a given area (the size of the sample A) is small, not exceeding 18 (Table 1). On the other hand, the size of sample B, i.e. the number of particles in a given area of constructed model, is rather large, and is about $10^4$. So we need to compare high-order moments calculated for small and for large samples. Of course, for a sample as small as A the uncertainties of the estimators of the higher-order moments are rather large. Moreover, we cannot calculate these uncertainties

without making an assumption about the distribution function. If we assume, for example, a Gaussian distribution, then we can calculate these uncertainties (Joanes & Gill 1998), but if we do not assume a priori any distribution function then the uncertainties are formally infinite. Furthermore, in the general case all commonly used estimators of the sample skewness and kurtosis are biased (Joanes & Gill 1998). For small samples this bias can be very large, especially for the kurtosis (see tables 2 and 3 in Joanes & Gill 1998). Consequently, for sample A the bias can be rather large. The bias depends on the size of the sample and for sample B (which is rather large) the bias is negligibly small. If both of our samples are generated from the same distribution function then the expected values of kurtosis (and skewness) for these two samples can differ significantly. We can solve this problem reducing sample ‘B’ to a smaller sample size.

Let us now describe how we calculate kurtosis. The skewness is calculated in the same way. Let us have some estimator of sample kurtosis. For sample A we simply calculate the value of kurtosis $K_a$ using the chosen estimator. For sample B we calculate the value of the kurtosis ‘reducing’ the sample size to that of sample A, i.e. $n_a$. We randomly get $n_a$ members from sample B and calculate for this sub-sample the kurtosis $k_1$. We repeat this $N_k$ times and calculate $k_1$, $k_2...k_{N_k}$ (in this article $N_k = 100000$). We denote the mean value of $k_1$ as $K_b^{\text{na}}$ and the standard deviation of $k_1$ as $\sigma_b^{\text{na}}$. Let $f_b(x)$ be the distribution function corresponding to sample B. For the chosen kurtosis estimator we can construct the distribution function $g_b(x)$ of the sample kurtosis for samples with size $n_a$ generated from $f_b(x)$. We note that $K_b^{\text{na}}$ is approximately equal to the expected value of $g_b(x)$ and $\sigma_b^{\text{na}}$ is approximately equal to the standard deviation of $g_b(x)$. For example, if the sample A was also generated from the distribution function $f_b(x)$, then $K_b^{\text{na}}$ would be the expected value and $\sigma_b^{\text{na}}$ the standard deviation of the kurtosis of sample A.

We use $K_a$ as measure of the kurtosis of sample A, $K_b^a$ as measure of the kurtosis of sample B and $\sigma_b^a$ as measure of the standard deviation. We note that the main value for us is $\tau = |K_a - K_b^a|/\sigma_b^{\text{na}}$. Joanes & Gill (1998) discussed three different estimators of the kurtosis. For our analysis, however, it does not matter which of the three we use because it can be proven that the value of $\tau$ is the same for all of them. To calculate the sample skewness and the sample kurtosis we use the estimators $G_1$ and $G_2$, discussed in Joanes & Gill (1998). These estimators are unbiased in the case of a normal distribution.

As we noted above, we cannot calculate uncertainties of estimators of kurtosis or skewness without making an assumption about the distribution function. The value $\sigma_b^a$ is the standard deviation of the estimation of the kurtosis or the skewness for sample A (observation) assuming that it is from the same distribution function as sample B (model). So this value depends on the model. We also can calculate the standard deviations $\sigma_n^a$ of the estimation of kurtosis or skewness for sample A assuming that its distribution function is Gaussian (Joanes & Gill 1998). For our models in general, the value of $\sigma_n^a$ is less than that of $\sigma_b^a$. The difference, however, is not so strong. The maximum difference is $\sigma_b^a/\sigma_n^a \approx 1.3$ and $\sigma_b^a/\sigma_n^a \approx 1.4$ for the skewness and the kurtosis, respectively.

This paper has been typeset from a TeX/Latex file prepared by the author.