text:
# Superfluid Hydrodynamics of an Electron Gas in a Superstrong Magnetic Field

\[

## Abstract

We derive the equations of hydrodynamics of a fully polarized electron gas placed in a strong magnetic field. These equations reveal the existence of solitons – immobile or propagating domain wall-like defects whose plane is perpendicular to the field direction. The solitons are used to construct weakly excited stationary states, and novel nonuniform persistent current states of the system.

\]

Understanding the properties of matter subject to an extremely strong magnetic field  has fundamental importance. At fields of order $`10^9G`$ the interaction of atomic electrons with the external magnetic field begins to dominate over their Coulomb interaction which results in a novel kind of atomic physics . In semiconductors and semimetals due to the smallness of the effective mass of the charge carriers and dielectric screening of the Coulomb interaction , the corresponding field is substantially smaller, and accessible in the laboratory .

Since a uniform magnetic field localizes the motion of a charged particle in the plane perpendicular to the field direction, and does not affect the motion along the field , strong fields can be used to study essentially one-dimensional physics using three-dimensional samples. The one-dimensionality is most pronounced if the field is strong enough to pull all the carriers into the lowest Landau level. In this limit in an interacting system a variety of phases of one-dimensional nature such as charge-density waves , Wigner crystals , and excitonic insulators  have been predicted, and a uniform-density state (adiabatically evolving from the noninteracting limit) occupies only a portion of the phase diagram . Recent work  has demonstrated that this phase has Luttinger-liquid-like properties.

In this Letter, by studying elementary excitations, we show that the uniform-density phase behaves as a unidirectional superfluid. The governing equations of hydrodynamics have a rich mathematical structure which owes its existence to the one-dimensional physics of the problem. A detailed treatment will be presented elsewhere . Below we restrict ourselves only to the salient physics effects.

We start with reminding the reader of the method of constructing the ground state of a noninteracting electron gas of density $`n`$ placed in a fully polarizing magnetic field $`H`$ pointing along the $`x`$ direction of a coordinate system. The motion of electrons is confined within infinitely long Landau tubes parallel to the field direction: the motion perpendicular to the tube axis is localized within the magnetic length $`l=(\mathrm{}c/eH)^{1/2}`$ while the motion along the tube is free. The degeneracy of the lowest Landau level, i.e. the number of Landau tubes in the system of cross-sectional area $`S`$ perpendicular to the field direction is given by  $`eHS/2\pi \mathrm{}c=S/2\pi l^2`$ which implies that the cross-sectional area of the individual tube is $`2\pi l^2`$. Therefore the magnetic flux confined within the tube is the flux quantum $`hc/e`$, and thus each Landau tube can be thought of as a magnetic field line. Hence the quantum state of every electron can be characterized by specifying the Landau tube it belongs to, and the wavevector $`k_x`$ describing the free motion along the field. Because of the Pauli principle, for a given Landau tube all the states with $`|k_x|k_F`$ are occupied while all the states having $`|k_x|>k_F`$ are empty. The Fermi wave vector $`k_F`$ can be found as $`k_F=\pi \rho `$ where $`\rho `$ is the linear electron density in the tube, which can be computed by dividing the total number of electrons $`N`$ by the number of Landau tubes $`S/2\pi l^2`$, and by the tube length $`L`$. Thus $`\rho =2\pi l^2n`$, $`k_F=2\pi ^2l^2n`$ and the Fermi energy can be found as $`ϵ_F=\mathrm{}^2k_F^2/2m=(2\pi ^2\mathrm{}l^2n)^2/2m`$ where $`m`$ is the electron mass. The expression for the ground-state energy

$$E_0=\frac{V}{(2\pi l)^2}\underset{k_F}{\overset{k_F}{}}\frac{\mathrm{}^2k_x^2}{2m}𝑑k_x=\frac{V(2\pi ^2\mathrm{}l^2)^2n^3}{6m}$$
(1)
(where $`V=SL`$ is the system volume) can be rewritten as $`(L\pi ^2\mathrm{}^2/6m)(2\pi l^2n)^3`$ times $`S/2\pi l^2`$: the ground-state energy of a one-dimensional spinless Fermi gas of linear density $`\rho =2\pi l^2n`$, multiplied by the total number of Landau tubes. The condition of being at the lowest Landau level is satisfied if the Fermi energy lies below the bottom of the second Landau level $`\mathrm{}eH/mc`$ . This leads to the condition $`nl^31/\pi ^2\sqrt{2}`$ which gives the range of applicability of Eq.(1). In what follows we will adopt the stronger condition $`nl^31`$ which also allows us to neglect transitions to the higher Landau levels.

The presence of the Coulomb interaction introduces a new length scale into the problem, the “exciton” Bohr radius $`a=ϵ\mathrm{}^2/me^2`$, where $`ϵ`$ is the dielectric constant of the medium. In astrophysical applications when $`ϵ=1`$, and $`m=m_e`$, the electron mass in vacuum, $`a`$ reduces to Bohr’s radius. In the limit $`nl^31`$ the average distance between the electrons belonging to the same Landau tube, $`1/2\pi l^2n`$, is much bigger than the average distance between the nearest tubes $`l`$. Thus the electrons belonging to the same Landau tube interact much more weakly among themselves than the electrons belonging to different tubes. Therefore the physical picture of independent Landau tubes survives if the Fermi energy is much bigger than the Coulomb interaction of electrons confined to different tubes. The latter (in an electrically neutral system) is estimated  as the plasmon energy $`\mathrm{}\omega _p`$ where $`\omega _p^2=4\pi ne^2/mϵ_0`$ is the plasma frequency. Combining the requirement $`ϵ_F\mathrm{}\omega _p`$ with the condition that all electrons are at the bottom of the lowest Landau level, we arrive at the concentration interval

$$(a/l)^{8/3}na^3(a/l)^3$$
(2)
which is the range of applicability of our theory. The interval (2) exists if the exciton Bohr radius $`a`$ is much bigger than the magnetic length $`l`$; the magnetic fields satisfying this condition are refered to as superstrong . Both length scales coincide at the field $`H_0=m^2e^3c/ϵ^2\mathrm{}^3`$. In astrophysical applications ($`m=m_e`$, $`ϵ=1`$) $`H_0=2.3510^9G`$, while in condensed matter applications the corresponding field is $`(ϵm_e/m)^210^410^6`$ times smaller, and for $`a/l510`$ the density range (2) exists in doped semiconductors and semimetals.

We postulate that the long-wavelength physics of a polarized electron gas placed in a slowly varying external potential $`V(𝐫)`$ is described by the energy functional of the complex order parameter $`\mathrm{\Psi }(𝐫,t)`$:

$$F=𝑑V\left[\frac{\mathrm{}^2}{2m}\left|\frac{d\mathrm{\Psi }}{dx}\right|^2+V(𝐫)|\mathrm{\Psi }|^2+\frac{(2\pi ^2\mathrm{}l^2)^2}{6m}|\mathrm{\Psi }|^6\right]$$
(3)
The physical meaning of $`\mathrm{\Psi }`$ is that $`|\mathrm{\Psi }|^2=n(𝐫,t)`$ gives the electron density which is a function of position $`𝐫`$ and time $`t`$. If $`\mathrm{\Psi }`$ is restricted to be real, then (3) can be recognized as a starting point of a generalized Thomas-Fermi density functional theory : for the translationally-invariant case the density is uniform, and then (3) reduces to the ground-state energy (1). The second term of (3) is the energy in the external field. The derivative term, known as “Weizsacker inhomogeneity correction” , which is a lowest-order term of a gradient expansion, takes into account the fact that deviations from uniformity in the field direction are costly. Density functionals similar to (3) are used to describe static properties of multielectron atoms placed in very strong magnetic fields .

We further postulate that the order parameter $`\mathrm{\Psi }`$ evolves according to the equation of motion

$$i\mathrm{}_t\mathrm{\Psi }=\frac{\delta F}{\delta \mathrm{\Psi }_{}}=\left[\frac{\mathrm{}^2_x^2}{2m}+V(𝐫)+\frac{(2\pi ^2\mathrm{}l^2)^2|\mathrm{\Psi }|^4}{2m}\right]\mathrm{\Psi }$$
(4)
In the representation $`\mathrm{\Psi }=n^{1/2}\mathrm{exp}i\theta `$ and $`v=(\mathrm{}/m)_x\theta `$, Eq.(4) reduces to the system of two coupled equations:

$$_tn+_x(nv)=0$$
(5)
$$m(_tv+v_xv)=_x[V+\frac{(2\pi ^2\mathrm{}l^2)^2}{2m}n^2\frac{\mathrm{}^2_x^2\sqrt{n}}{2m\sqrt{n}}]$$
(6)
In the long-wavelength limit when the last term in the right-hand side of (6) can be dropped, the system of equations (5) and (6) can be identified with those of hydrodynamics of an ideal fluid  whose pressure and density are related by the equation of state $`p=(2\pi ^2\mathrm{}l^2)^2n^3/3m`$. The variable $`v=(\mathrm{}/m)_x\theta `$ has a meaning of the superfluid velocity which only has a component along the field direction. We note that the equations of hydrodynamics of an electron gas in a superstrong magnetic field could be written down from the outset and the density dependence of pressure can be computed from the ground-state energy (1) as $`p=dE/dV`$. This provides an important test of the underlying equations (3) and (4); specifically, selection of the gradient term of (3) in the Weizsacker form plays a central role in getting the correct hydrodynamical limit. We stress that hydrodynamical description requires two fields, density and velocity; that is why the order parameter $`\mathrm{\Psi }`$ must be complex.

Spectrum of small density oscillations: \- Let us look for solutions to (4) (with $`V=0`$) of the form $`\mathrm{\Psi }(x,t)=\psi (x,t)\mathrm{exp}(i\mu t/\mathrm{})`$. The function $`\psi (x,t)`$ then obeys the equation

$$i\mathrm{}_t\psi =\frac{\mathrm{}^2}{2m}\left[_x^2+(2\pi ^2l^2)^2(|\psi |^4\psi _0^4)\right]\psi $$
(7)
where the chemical potential $`\mu =\pi ^2\mathrm{}^2(2\pi l^2)^2\psi _0^4/2m`$ is selected so that the ground state corresponds to the uniform electron density $`n_0=\psi _0^2`$. Consider a small oscillation of $`\psi `$ around its constant value $`\psi _0`$: $`\psi =\psi _0+A\mathrm{exp}[i(qx\omega t)]+B^{}\mathrm{exp}[i(qx\omega t)]`$, where $`A`$ and $`B^{}`$ are small complex amplitudes, and $`q`$ is the wavevector along the field direction. Substituting the expression for $`\psi `$ into (7), linearizing, and solving the resulting system of two equations, we find for the dispersion law

$$(\mathrm{}\omega )^2=(\mathrm{}^2q^2/2m)^2+(\mathrm{}cq)^2$$
(8)
The spectrum takes the Bogoliubov form : for large momenta the dispersion law coincides with the free-particle energy $`\mathrm{}^2q^2/2m`$, while at low momenta we get instead the phonon dispersion $`\omega =cq`$, where $`c=2\pi ^2l^2\mathrm{}n_0/m`$ is the sound velocity. We note that this velocity coincides with the hydrodynamical expression $`c=(dp/mdn)^{1/2}`$ for a gas with equation of state $`p=(2\pi ^2\mathrm{}l^2)^2n^3/3m`$. The Landau criterion applied to (8) then implies  that the critical velocity of superfluid motion along the direction of magnetic field equals the sound velocity.

Solitons: \- Without the external potential ($`V=0`$) Eq.(7) has soliton solutions . These describe moving (gray soliton) or standing (dark soliton) density depressions accompanied by a corresponding profile of phase . Because of the system uniformity in the plane perpendicular to the magnetic field, the solitons are domain wall-like defects whose plane is perpendicular to the direction of magnetic field; the solitons can only propagate along (or opposite) to the field direction (Figure 1).

In dimensionless variables $`\nu =n/n_0`$, $`y=2\pi ^2l^2n_0x`$, $`\tau =(2\pi ^2l^2n_0)^2\mathrm{}t/m`$, $`\zeta =y\beta \tau `$, the density ($`\nu `$) and phase ($`\theta `$) profiles for the soliton are given by

$$\nu =1\frac{3(1\beta ^2)}{2+(1+3\beta ^2)^{1/2}\mathrm{cosh}[2(1\beta ^2)^{1/2}\zeta ]}$$
(9)
$$2\theta =\mathrm{cos}^1\left[\frac{(3\beta ^2/\nu )1}{(1+3\beta ^2)^{1/2}}\right]$$
(10)
where the dimensionless velocity $`\beta `$ is measured in units of the sound velocity. These formulas imply that as $`\beta `$ varies between zero and unity, the soliton becomes more shallow, and the phase shift across the soliton $`\mathrm{\Delta }\theta =\mathrm{cos}^1[(3\beta ^21)/(1+3\beta ^2)^{1/2}]`$ changes from $`\pi `$ ($`\beta =0`$) to zero ($`\beta =1`$); the soliton solution ceases to exist for velocities exceeding the sound velocity. Microscopically the domain wall is composed of aligned point-like solitons belonging to different Landau tubes, and there is a positive energy per unit area associated with the planar defect. It can be calculated as the energy of the soliton belonging to the Landau tube (computed in ), divided by the cross-sectional area of the Landau tube $`2\pi l^2`$:

$$\sigma (\beta )=\frac{\sqrt{3}(\pi l\mathrm{}n_0)^2}{m}(1\beta ^2)\mathrm{ln}\left[\frac{2+[3(1\beta ^2]^{1/2}}{(1+3\beta ^2)^{1/2}}\right]$$
(11)
The positivity of the surface energy (11) guarantees the stability of the domain wall (9), (10) against transverse fluctuations.

Stationary states: \- Solitons form a natural language to describe weakly excited states of the system. For example, let us look at electron gas placed inside a layer of macroscopic thickness $`L`$; the plane of the layer is perpendicular to the direction of the magnetic field. If the electrons cannot leave the layer, then in the ground state the density profile will look as shown schematically in Figure 2 (left): the density stays constant almost everywhere except for the narrow vicinity of the walls where it sharply plunges to zero.

It is convenient to think of these two dips as having one dark soliton present in the system. It is also convenient to take the convention that what is shown in Figure 2 (left) is a first excited state with respect to the strictly uniform state which we will call the ground state. Thus the first excited state has just one soliton present in the system. Then the second excited state will have two solitons \[Figure 2 (right)\]. The argument can be extended to any number of solitons. For example, if there are Q dark solitons present, then as long as soliton tails do not overlap, $`Ql^2n_0L`$, the system energy is given by $`E(Q)=E_0+S\sigma (0)Q`$ where $`E_0`$ is the ground-state energy (1). A similar conclusion is applicable to the case of a ring-shaped sample with the only constraint that now $`Q`$ must be even because the phase shift due to the dark soliton is $`\pm \pi `$. Having an even number of dark solitons will guarantee that the overall change of phase of the order parameter upon going around the ring is multiple of $`2\pi `$, and thus the order parameter is unique.

Because translational symmetry is not broken for the ring geometry, dark soliton states do not exhaust all possible excited states of the system: gray solitons can also be involved. However it will be easier to anticipate their role after discussing the persistent current states.

Persistent currents: \- Persistent currents are only possible in the ring geometry. Experimentally the magnetic field directed along the circumference of a ring-shaped sample can be generated by an electric current flowing along the symmetry axis perpendicular to the plane of the ring. In the ring geometry the phase of the order parameter can only change by a multiple of $`2\pi `$. The simplest way to accommodate a nonzero phase change is to have a constant phase gradient along the ring . This will lead to a spatially uniform persistent current state in which all the electrons are moving around the ring with the same velocity. This state can be also viewed as a stationary state.

There is also a novel nonuniform way of accommodating the phase difference. The corresponding persistent current states can be constructed from gray solitons, and they have a series of remarkable properties. For example, imagine we want to insert into the ring $`Q`$ solitons to accommodate an overall phase shift of $`2\pi `$. Then the phase shift due to the individual soliton must be $`2\pi /Q`$. This will in turn determine the soliton velocity $`\beta `$ which will then fix the surface energy (11). An example of four gray solitons moving along the ring is sketched in Figure 3.

The energy of the train consisting of $`Q`$ solitons running around the ring of circumference $`L`$ is given by $`E(Q)=E_0+S\sigma [\beta (Q)]Q`$ as long as soliton tails do not overlap, $`Q(1\beta ^2)^{1/2}l^2n_0L`$. In the reference frame moving with the soliton train the solitons will be standing but the background electron gas will be moving with a constant velocity. This gives an example of a stationary state involving gray solitons.

For the uniform moving state the total angular momentum of the system is quantized in units of $`N\mathrm{}`$ where $`N`$ is the total particle number . We speculate that for the soliton train moving state the total angular momentum might be a fraction of $`N\mathrm{}`$. Our contention is based on the observation that while the soliton is moving, the electron motion is localized only in the vicinity of the soliton core - merely a tiny fraction of particles are involved in the motion.

In applying these results to the real system one has to keep in mind that since a soliton is a density depletion, for solitons which do not travel fast enough, the density at the soliton core can drop below the left bound of (2). Therefore our results are strictly applicable only to solitons whose velocity is larger than some critical velocity. Future work is necessary to understand how interactions between the electrons belonging to different Landau tubes affect the slow solitons. One possibility is that they convert the depleted soliton core into a two-dimensional Wigner crystal. Another possibility is a modulational instability of the soliton in which Coulomb repulsion of the electrons decreases at the expense of increase of the surface energy (11).

The unique mathematical structure of the underlying equations (3)-(6) which so far did not enter the discussion deserves a special commentary. First we note that in the hydrodynamical limit when the last term of (6) can be neglected (and $`V=0`$) the problem of arbitrary motion of a one-dimensional gas with the equation of state $`pn^3`$ can be solved in a closed form . Second, (even for $`Vx^2`$) equations (4)-(6) have self-similar solutions  which only exist for the equation of state $`pn^3`$. Physically some of them describe free expansion of the electron gas initially confined along the field direction .

The author thanks T. J. Newman for valuable assistance, and X. Qi for useful discussions. This work was supported by the Thomas F. Jeffress and Kate Miller Jeffress Memorial Trust.


format: text

text:
# DYNAMICS OF STARS AND GLOBULAR CLUSTERS IN M87

## 1 INTRODUCTION

Elliptical galaxies lack the dynamical simplicity of spiral galaxies, posing well-known challenges for determining their intrinsic properties (e.g., mass distribution, shape, orbit structure). The projected stellar velocity dispersion radial profile, $`\widehat{\sigma }_\mathrm{p}(R)`$, is easily measured, but is subject to a degeneracy between the mass distribution and the orbital anisotropy, where radial variations in the stellar orbit types mask or mimic changes in the mass-to-light ratio ($`M/L`$) (Binney & Mamon 1982; Tonry 1983). Crucial for resolving this degeneracy is the measurement of higher-order moments of the stellar line-of-sight velocity distributions (LOSVDs). Such LOSVD information has been used to constrain the masses of central supermassive black holes (van der Marel et al. 1998; Emsellem, Dejonghe, & Bacon 1999; Cretton & van den Bosch 1999; Gebhardt et al. 2000) and dark halos (Carollo et al. 1995; Rix et al. 1997; Gerhard et al. 1998; Saglia et al. 2000; Kronawitter et al. 2000). Unfortunately, the surface brightness of an elliptical declines rapidly with radius, making these measurements increasingly difficult at the large radii where the dark matter problem becomes most interesting ($`>`$ 2 effective radii $`R_{\mathrm{eff}}`$). Other techniques are necessary to probe this outer mass distribution.

In nearby galaxies, the two main candidates for such probes are X-ray emission from extended hot gas (see Fabbiano 1989 §4.4 for a review), and the discrete velocities of either planetary nebulae (e.g., Ciardullo, Jacoby, & Dejonghe 1993; Tremblay, Merritt, & Williams 1995; Hui et al. 1995; Arnaboldi et al. 1998) or globular clusters (e.g., Cohen & Ryzhov 1997, hereafter CR; Kissler-Patig et al. 1999; Zepf et al. 2000). However, the analysis of discrete kinematical measurements calls for more refined dynamical models than either simple virial estimators (e.g., Heisler, Tremaine, & Bahcall 1985; Kent 1990; Haller & Melia 1996) or binning the data and using the Jeans equations (e.g., Federici et al. 1993; Grillmair et al. 1994; Tremblay et al. 1995). These procedures do not guarantee a physical solution, make unwarranted assumptions about the orbital anisotropy, and do not make full use of the constraints provided by the data. The first two deficiencies are alleviated by techniques employing distribution function (DF) basis sets that permit the full freedom of orbit types (e.g., Merritt, Meylan, & Mayor 1997; Mathieu & Dejonghe 1999; Saglia et al. 2000). But these techniques typically bin the data in radius and velocity, destroying potentially useful information, such as the maximum observed velocity as a direct constraint on the escape velocity. More general methods are needed to make better use of the information contained in the discrete velocity data; Merritt & Saha (1993) and Merritt (1993) have demonstrated one such approach using DF basis sets and maximum likelihood methods. We develop a general dynamical method based on orbit modeling, and apply it to a real galaxy, M87.

The giant central Virgo Cluster galaxy M87 (= NGC 4486), although not a typical elliptical, is nevertheless a prime candidate for modeling because of the abundance of data available on it. In addition to high-precision measurements of its higher-order stellar velocity moments and of the radially-extended velocity dispersion, it has the largest available sample of globular cluster (GC) velocity measurements. The stellar velocity measurements extend out to $`1.5R_{\mathrm{eff}}`$, and the GC measurements to $`5R_{\mathrm{eff}}`$, probing well into the area where the dynamics are presumably dominated by the dark halo of the galaxy or of the Virgo Cluster. The combination of these kinematical constraints should provide robust limits on the dark matter distribution, and on the orbital structures of both the stars and the GCs.

While a variety of modeling techniques have been used to study the stellar dynamics of M87 (Sargent et al. 1978; Duncan & Wheeler 1980; Binney & Mamon 1982; Newton & Binney 1984; Richstone & Tremaine 1985; Dressler & Richstone 1990; Tenjes, Einasto, & Haud 1991; van der Marel 1994; Merritt & Oh 1997), they concentrated on its central regions in order to constrain the mass of a central black hole. The stellar dynamics of its halo have not been modeled in detail, nor have the higher-order velocity moments been used as strict constraints. But the best available tracers of the halo mass distribution are ultimately the GCs and the hot gas, which can be easily observed at large radii. Some simple dynamical methods have been used to estimate the halo mass using the GCs (Merritt & Tremblay 1993; CR). They claim that the radially-rising velocity dispersion profile rules out a constant $`M/L`$ model, and that it suggests that the GCs are in fact tracing the Virgo Cluster potential rather than M87’s own potential (as suggested also by analyses of the X-ray emission and the kinematics of the cluster galaxies; e.g., see Nulsen & Böhringer 1995; McLaughlin 1999). But a more rigorous dynamical model, allowing for the systematic orbital uncertainties, is necessary to confirm these findings. We will apply our general orbit modeling method, using a large sample ($`200`$) of discrete GC velocities, to see how the GCs — alone and jointly with the stellar dynamics — constrain the mass distribution. In §2 we detail the observational constraints for M87. We describe our modeling methods in §3, §4 presents the results, and the conclusions are in §5.

## 2 OBSERVATIONAL CONSTRAINTS

Here we summarize the observational constraints on M87. These include the stellar surface brightness (§2.1); velocity measurements of the stars (§2.2); the globular cluster surface density (§2.3); and globular cluster discrete velocity measurements (§2.4). We assume a distance to M87 of 15 Mpc (Pierce et al. 1994), so that $`1{}_{}{}^{\prime \prime }=73`$ pc and $`1{}_{}{}^{}=4.4`$ kpc. The stellar effective radius is $`R_{\mathrm{eff}}100^{\prime \prime }`$ (Burstein et al. 1987; Zeilinger, Møller, & Stiavelli 1993).

### 2.1 Stellar Surface Brightness

For the stellar surface brightness radial profile $`\mu (R)`$ in the core region ($`R`$=0.′′01–8′′) we use the HST $`I`$-band data of Lauer et al. (1992). We use their seeing-deconvolved profile, and approximate the photometric uncertainties by $`\mathrm{\Delta }\mu =\mathrm{max}[0`$.′′$`011/R,0.02]`$. For the outer regions ($`R`$=8′′–745′′), we use the $`B`$-band data of Capaccioli, Caon, & Rampazzo (1990), with the photometric uncertainties estimated from their Figure 3, and with an offset of 2.78 magnitudes to match the profile to the Lauer et al. (1992) data. The color gradients are quite small over the entire range of the galaxy ($`|\mathrm{\Delta }(BI)|_{\mathrm{max}}<0.2`$ mag; e.g., Boroson et al. 1983; Cohen 1986; Zeilinger et al. 1993), so there should be little systematic problem in combining the $`I`$-band and $`B`$-band data. The isophote ellipticity varies from 0.02 at $`2^{\prime \prime }`$ to 0.1 at $`75^{\prime \prime }`$ to 0.35 at $`500^{\prime \prime }`$. The dynamics of the galaxy depend on the ellipticity of its gravitational potential, which is much rounder than its density distribution ($`ϵ_\mathrm{\Phi }ϵ_\rho /3`$; Binney & Tremaine 1987 §2.3.1), so it will be a reasonable approximation to treat the galaxy as spherical. We map all of the data from the major axis $`a`$ to the intermediate axis $`m=a\sqrt{1ϵ}`$.

For the purposes of our mass modeling, we project and fit a parameterized luminosity density model to the data of the form $`\nu _{}(r)=\nu _0(r/s_1)^{\alpha _1}(1+r^3/s_2^3)^{\alpha _2}(1+r/s_3)^{\alpha _3}`$, and find $`\nu _0=9.9\times 10^6L_{B,}`$ arcsec-3, $`\alpha _1=1.22`$, $`\alpha _2=0.361`$, $`\alpha _3=1.26`$, $`s_1=1^{\prime \prime }`$, $`s_2=9`$.′′$`38`$, and $`s_3=170^{\prime \prime }`$, with a $`\chi ^2`$ statistic for the fit of 35.1 for 142 degrees of freedom. This functional fit will only be used to generate a constant $`M/L`$ gravitational potential, while we will directly fit the surface brightness data in our dynamical models (§3). Note that in M87, $`\mu (R)`$ does not taper off steeply in the outer regions as in normal $`R^{1/4}`$ ellipticals — a typical $`R^2`$ cD envelope (Schombert 1986) begins at $`R>300^{\prime \prime }`$. An old dispute over the brightness of this envelope (see e.g., de Vaucouleurs & Nieto 1978; Carter & Dixon 1978) has not been resolved, presumably because of sky-subtraction problems for such a large, low-surface-brightness structure. Our asymptotic radial slope lies between the reported extremes.

### 2.2 Stellar Velocities

For the core region ($`R<29^{\prime \prime }`$), we use the G-band data from van der Marel (1994, hereafter vdM) for the projected stellar velocity dispersion $`\widehat{\sigma }_\mathrm{p}(R)`$ and the Gauss-Hermite velocity moments $`h_4(R)`$ and $`h_6(R)`$ (see §3). For the outer regions ($`R`$=28′′–168′′), we use $`\widehat{\sigma }_\mathrm{p}(R)`$ from Sembach & Tonry (1996, hereafter ST). We combine the data from positive and negative radii into one radial profile. There appears to be a systematic velocity offset between the ST data and most other data (Sargent et al. 1978; Davies & Birkinshaw 1988; Jarvis & Peletier 1991; Winsall & Freeman 1993), presumably due to the large slit width; we estimate that this corresponds to an additional instrumental dispersion of $`183\pm 11`$ km s-1, which we remove from the ST data in order to match them to the vdM data. The final velocity dispersion profile is nearly constant for $`R<1^{\prime \prime }`$, outside of which it falls off slowly with radius ($`\widehat{\sigma }_pR^{0.1}`$). Comparing $`h_4`$ and $`h_6`$ from positive and negative radii, we find differences larger than are consistent with the stated uncertainties. We assume the errors were underestimated by 12% and 22% to bring the profiles into statistical agreement. The departures from Gaussianity of the LOSVD are everywhere small (typically, $`|h_4|,|h_6|<0.02`$).

### 2.3 Globular Cluster Surface Density

Our GC surface density radial profile $`N(R)`$ is taken from the number counts of Kundu et al. (1999) for $`R=0^{\prime \prime }`$$`96^{\prime \prime }`$ (with their quoted uncertainties). For $`R=84^{\prime \prime }`$$`472^{\prime \prime }`$, we take the data from McLaughlin, Harris, & Hanes (1993), adding in their 0.4 arcmin-2 background uncertainty, and multiplying by 1.45 to match the normalization of the Kundu et al. (1999) data. For $`R=419^{\prime \prime }`$$`1351^{\prime \prime }`$, we use the data from Harris (1986) after subtracting their background count of $`5.8\pm 0.3`$ arcmin-2 and normalizing by a factor of 2.19. We do not convert any of these profiles to an intermediate axis since this is already effectively accomplished by the derivation method (counting in circular annuli). Like the stellar surface brightness, the GC surface density decreases slowly with radius in the central regions ($`R^{0.3}`$), and changes over to a steeper power law in the outer parts ($`R^{1.6}`$); however, the radius where this break occurs is much larger for the GCs ($`60^{\prime \prime }`$ vs. $`10^{\prime \prime }`$), demonstrating that these systems are dynamically distinct.

### 2.4 Globular Cluster Discrete Velocities

For the GC line-of-sight velocities $`v_z`$, we use data from Mould, Oke, & Nemec (1987), Huchra & Brodie (1987), Mould et al. (1990), CR, and Cohen (2000). Most of the positions are taken from Strom et al. (1981). We compare the velocities of common objects to determine the systematic offset between the data sets, and to estimate the measurement uncertainties ($`\mathrm{\Delta }v_z100`$ km s-1). We discard as foreground stars all objects with heliocentric velocities $`v_z<250`$ km s-1, and as background galaxies all those with $`v_z>2650`$ km s-1. We examine the colors and magnitudes of objects with $`v_z<500`$ km s-1 to distinguish GCs from stars. Objects of uncertain identification are discarded. We map the data to the intermediate axis $`m`$, using ellipticity estimates from McLaughlin, Harris, & Hanes (1994). Our final data set has 234 velocities from $`R=25^{\prime \prime }`$$`526^{\prime \prime }`$ (see Fig. 1).

FIG. 1.— Line-of-sight velocities for 234 globular clusters around M87, as a function of galactocentric radius $`R`$. The dotted line shows a systemic velocity of 1300 km s-1.

There is a paucity of low velocity measurements at large radii, which may be related to the rotation of the GCs along the galaxy’s major axis in the outer parts, as demonstrated by Kissler-Patig & Gebhardt (1998). For the GCs with $`R370^{\prime \prime }`$, there were only 4 measurements taken on the “approaching” side and 23 on the “receding” side. Thus, the large-radius velocity asymmetry is probably due to incomplete spatial coverage, which may make our mean velocity for the system systematically high (it is decreased from 1352 to 1293 km s-1 if the clusters with $`R370^{\prime \prime }`$ are omitted). This is also suggested by the global velocity distribution (Fig. 2), which appears more symmetric if the mean velocity is set to $`1300`$ km s-1. The best-fit Gaussian curve to the data has a central velocity of $`\widehat{v}_\mathrm{p}=1323\pm 28`$ km s-1 ($`1308\pm 29`$ km s-1 if the outer clusters are omitted). We adopt a systemic velocity of 1300 km s-1, which we subtract from the measurements to get the galactocentric velocities. This compares well with estimates of the stellar mean redshift from vdM (1277 km s-1), ST (1293 km s-1), and the RC3 ($`1282\pm 9`$ km s-1; de Vaucouleurs et al. 1991). The large-radius velocity asymmetry will not pose any further difficulties because we will model only the even part of the DF, so that all the discrete velocities are implicitly considered as present also at their reflected position about the mean velocity.

Fig. 2 also suggests that the GCs have a tangentially-biased orbit distribution, which typically produces a flat or double-peaked LOSVD. We fit a fourth-order Gauss-Hermite velocity moment (see §3) to the system, finding $`h_4=0.04\pm 0.02`$ ($`\widehat{\sigma }_\mathrm{p}=404\pm 16`$ km s-1), which is weakly suggestive of tangential anisotropy. A double-peaked LOSVD (especially evident at larger radii) has also been found in the GC system of the Fornax Cluster cD galaxy NGC 1399 (Kissler-Patig et al. 1999).

FIG. 2.— Distribution of globular cluster line-of-sight velocities, relative to a systemic velocity of 1300 km s-1. The histogram shows the number of data points in velocity bins. The solid curve shows a smoothed superposition of the measurements. The dotted curve shows the best-fit Gaussian curve to the data ($`\sigma =376`$ km s-1). Note the double peaks of low significance at $`\pm 200`$ km s-1 and at $`\pm 500`$ km s-1.

## 3 METHODS

The method of orbit modeling was pioneered by Schwarzschild (1979) to reproduce the observed density distribution of a galaxy from a weighted library of representative particle orbits. Given a fixed gravitational potential $`\mathrm{\Phi }(r)`$, the solutions are nonparametric — no functional form is assumed for the DF — so that one need not worry that the models cannot fully explore the possible solutions. The DF is guaranteed to be nonnegative and physical. This method was extended by Richstone & Tremaine (1984) to include projected velocity dispersion information, and by Rix et al. (1997) to include higher-order velocity moments. We used similar techniques to examine the uncertainty in a high-redshift galaxy’s mass as implied by a measurement of its central velocity dispersion — an important step in determining $`H_0`$ from certain gravitational lens systems (Romanowsky & Kochanek 1999, hereafter Paper I). Other recent orbit modeling efforts have concentrated on applying the method to three-integral axisymmetric systems (e.g., van der Marel et al. 1998; Cretton & van den Bosch 1999; Gebhardt et al. 2000; Cretton, Rix, & de Zeeuw 2000). Here we extend our spherical method to include the discrete velocities of GCs. The outline of the method is as follows.

We begin with an assumed radial profile for the mass density distribution $`\rho (r)`$. Our first model is a simple representation of a dark-matter dominated system, the singular isothermal sphere (SIS): $`\rho (r)=\sigma _0^2/2\pi Gr^2`$. The second model is the constant $`M/L`$ distribution described in §2.1, parameterized by $`\mathrm{{\rm Y}}_B`$: the ratio of the mass density $`\rho (r)`$ and the $`B`$-band luminosity density $`\nu _{}(r)`$, with units of $`M_{}/L_{,B}`$. We pick a random distribution of orbits that densely samples the phase space of energy and angular momentum $`(E,|L|)`$. The initial radii $`r_{0k}`$ of the orbits are logarithmically spaced in $`r_0`$, and the energy $`E_k`$ of each orbit is selected to correspond to that of a circular orbit at this radius, $`\mathrm{\Phi }(r_{0k})+v_\mathrm{c}^2(r_{0k})/2`$. For the SIS model, this procedure produces a sampling that is uniform in energy. The angular momentum $`L_k`$ of the orbit is selected randomly from the range $`[0,L_{\mathrm{max},k}]`$, where $`L_{\mathrm{max},k}=r_{0k}v_\mathrm{c}(r_{0k})`$. For both the stellar models and the GC models, we use 2500 particle orbits. For the stars, $`r_0`$ spans the range $`0`$.′′$`07`$$`4450^{\prime \prime }`$, resulting in a radial coverage of $`r=0^{\prime \prime }`$$`7330^{\prime \prime }`$ (SIS model) or $`r=0^{\prime \prime }`$$`8710^{\prime \prime }`$ (constant $`M/L`$ model). For the GCs, $`r_0=`$10.′′0–$`13360^{\prime \prime }`$, corresponding to $`r=0^{\prime \prime }`$$`22030^{\prime \prime }`$ (SIS model) or $`r=0^{\prime \prime }`$$`26420^{\prime \prime }`$ (constant $`M/L`$ model).

We next compute the orbit projection “kernels”, which correspond to the contribution of each orbit to each observable. The predictions of the model can then be expressed as follows:

$$y_i^\mathrm{m}=\underset{k}{}w_k^2K_{ik}_{t,\theta ,\varphi },$$
(1)
where the $`𝐰`$ are the orbit weights, and the kernels $`𝐊`$ have been averaged over time and over all spherical-polar viewing angles $`(\theta ,\varphi )`$. For example, the kernel for the angle-averaged surface density in an annulus between $`R_1`$ and $`R_2`$ of an orbit at an instantaneous radius $`r^{}`$ is given by the integral

$`K_{I(R_1,R_2)}_{\theta ,\varphi }(r^{})`$ $`=`$ $`{\displaystyle \frac{1}{\pi (R_2^2R_1^2)}}{\displaystyle _{R_1}^{R_2}}2\pi R𝑑R{\displaystyle _R^{\mathrm{}}}{\displaystyle \frac{\delta (rr^{})}{4\pi r^2}}{\displaystyle \frac{2rdr}{\sqrt{r^2R^2}}}{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\xi }{2\pi }}`$ (2)
$`=`$ $`{\displaystyle \frac{\sqrt{1R_1^2/r^2}\sqrt{1R_2^2/r^2}}{\pi (R_2^2R_1^2)}},`$ (3)
where the orbit’s instantaneous density is $`\delta (rr^{})/4\pi r^2`$, the variable $`\xi `$ represents the line of sight’s position in the tangential plane $`(\theta ,\varphi )`$, and the integration is carried out along the line of sight $`z=\pm \sqrt{r^2R^2}`$. The viewing-angle integral (along $`d\xi `$) is nontrivial for kernels that involve velocity measurements. Placed at the initial radius $`r_0`$, the orbit is run forward in time for one radial period $`T_r`$, and the final kernel is found by averaging over time:

$$K_{t,\theta ,\varphi }=\frac{1}{T_r}_0^{T_r}K_{\theta ,\varphi }(t)𝑑t,$$

where the integration is handled by a Bulirsch-Stoer type integrator.

We must calculate the model’s LOSVD, $`dL/dv_\mathrm{p}(v_\mathrm{p},R)`$, to fit the observed stellar velocity moments ($`\widehat{\sigma }_\mathrm{p},h_4,h_6`$) and the GC velocity measurements. The LOSVD kernels $`K_{dL/dv_\mathrm{p}(R),k}`$ for a measurement at radius $`R`$ are binned in velocity with $`v_\mathrm{p}=[0..v_{\mathrm{max}}]`$, where $`v_{\mathrm{max}}(R)`$ is the velocity of a particle at radius $`r=R`$ that has plunged inward on a radial orbit from the largest possible radius. For the constant $`M/L`$ model, this implies that $`v_{\mathrm{max}}v_{\mathrm{esc}}`$ (the escape velocity). Our models are symmetric in velocity, i.e., $`dL/dv_\mathrm{p}(v_\mathrm{p})=dL/dv_\mathrm{p}(v_\mathrm{p})`$. By construction, $`v_{\mathrm{max}}^2`$ scales linearly with the mass normalization of the galaxy model ($`\sigma _0^2`$ or $`\mathrm{{\rm Y}}_B`$), so that we can use the same binned LOSVD kernels when changing the mass scale.

The Gauss-Hermite velocity moments $`h_l`$, which measure deviations of a LOSVD from a Gaussian, are defined by

$$h_l\frac{\sqrt{2}\gamma _0}{\widehat{\gamma }_\mathrm{p}}_{\mathrm{}}^{\mathrm{}}\frac{dL}{dv_\mathrm{p}}(v_\mathrm{p})e^{\widehat{w}^2/2}H_l(\widehat{w})𝑑v_\mathrm{p},$$
(4)
where $`\widehat{w}=(v_\mathrm{p}\widehat{v}_\mathrm{p})/\widehat{\sigma }_\mathrm{p}`$, $`\gamma _0`$ is the line strength, $`(\widehat{\gamma }_\mathrm{p},\widehat{v}_\mathrm{p},\widehat{\sigma }_\mathrm{p})`$ are the coefficients for the best Gaussian fit to $`dL/dv_\mathrm{p}`$, and $`H_l(\widehat{w})`$ are the Hermite polynomials (van der Marel & Franx 1993). The moments ($`h_4`$, $`h_6`$) measure how flat ($`h_4<0`$, $`h_6>0`$) or how peaked ($`h_4>0,h_6<0`$) the LOSVD is compared to a Gaussian, where these signatures are typically produced by tangential and radial orbits, respectively. The two moments differ in how much weight they put in the high-velocity wings of the LOSVD. In Paper I, we fitted a Gaussian curve to the LOSVD to find ($`\widehat{\gamma }`$, $`\widehat{\sigma }_\mathrm{p}`$) at each step in our minimization routine, and then calculated $`h_4`$; this is equivalent to solving the nonlinear equations $`h_0=1`$ and $`h_2=0`$. Here, we expand the Gauss-Hermite series about the fixed data values $`\widehat{\sigma }_\mathrm{p}^\mathrm{d}`$. Rather than directly fitting $`\widehat{\sigma }_\mathrm{p}^\mathrm{d}`$, we equivalently fit the second-order moment $`h_2=0.00\pm \mathrm{\Delta }h_2`$, where $`\mathrm{\Delta }h_2\mathrm{\Delta }\widehat{\sigma }_\mathrm{p}/\sqrt{2}\widehat{\sigma }_\mathrm{p}^\mathrm{d}`$ . With this formulation, calculating the model’s velocity observables is now linear with respect to the orbit weights. We use 41 velocity bins from $`v_\mathrm{p}=0`$ to $`v_{\mathrm{max}}`$ to numerically integrate equation (4), which we have heuristically found to give very good accuracy. Note that because our model is completely spherically symmetric (simulating only the even part of the DF), it would be more appropriate to fit measurements of $`z_4`$ (expanded around $`\widehat{v}_\mathrm{p}=0`$) rather than $`h_4`$ (see van der Marel et al. 1994 §5.1); however, these measurements are not available for M87. Since the measured rotation is small ($`10`$ km s-1), the difference would be unimportant.

The kernels $`(𝐊_I,𝐊_{dL/dv_\mathrm{p}})`$ need be computed only once for each galaxy mass model. We then adjust the weights of the orbits so that the model’s projected observables $`𝐲^\mathrm{m}`$ (equation 1) best fit the data $`𝐲^\mathrm{d}`$. For most of the data ($`I,\widehat{\sigma }_\mathrm{p},h_4,h_6`$), we express the likelihood of the fit using the standard $`\chi ^2`$ statistic,

$$\chi ^2=\underset{i}{}\left(\frac{y_i^\mathrm{m}y_i^\mathrm{d}}{\sigma _i}\right)^2,$$
(5)
where the measurement uncertainty is $`\sigma _i`$, while the likelihood function for the discrete velocities and positions of the GCs is

$$_i(v_i,R_i)\frac{dL}{dv_\mathrm{p}}(v_\mathrm{p},R_i)e^{(v_iv_\mathrm{p})^2/2\sigma _i^2}𝑑v_\mathrm{p},$$
(6)
where $`v_i\pm \sigma _i`$ are the individual velocity measurements. In maximizing $`_i`$, we are forcing the LOSVD to peak at the measured velocities $`v_i`$ weighted by the measurement errors $`\sigma _i`$. This is schematically illustrated in Fig. 3. Given complete freedom, this method would produce a best-fit solution whose LOSVDs resembled $`\delta `$-functions at the measurements $`v_i`$. However, the model’s averaging in angle and time produces an intrinsic smoothing to the LOSVDs, and does not permit such unphysical solutions. For the current implementation, we use 15 velocity bins from $`v_\mathrm{p}=0`$ to $`v_{\mathrm{max}}`$ for the LOSVD at each radius in order to compute the integral in equation (6). Note that our method does not require binning the velocities in radius, nor computing velocity moments.

FIG. 3.— Schematic line-of-sight velocity distributions (solid curves), fitted to discrete velocity data (dotted curves). The double-peaked LOSVD on the left is indicative of tangential anisotropy, while the centrally-peaked LOSVD on the right indicates radial anisotropy.

The final function which we will minimize is

$$f(𝐰)\frac{1}{2}\chi ^2\underset{i}{}\mathrm{ln}_i+\lambda S,$$
(7)
where $`S`$ is a measure of entropy:

$$S=\underset{k}{}w_k^2\mathrm{ln}w_k^2.$$
(8)
We employ this smoothing function $`S`$ as a heuristic device for quickly reaching a rough solution; in our final solutions, the Lagrangian factor $`\lambda `$ is reduced to such a low value $`(10^5)`$ so as to make the entropy constraints inconsequential. With no such regularization imposed, this is an ill-conditioned inversion problem, so our DF solutions will be choppy in a way that real galaxies’ DFs presumably are not. However, our methods are the statistically correct way to handle the uncertainties, absent any other a priori smoothness conditions. We use a conjugate gradient method (see Press et al. 1992 §10.6) with first and second derivative information, to minimize $`f`$. We have tested our methods on an isotropic Hernquist (1990) model to demonstrate that we recover correctly the input mass, given constraints on $`I(R)`$ and $`\widehat{\sigma }_\mathrm{p}(R)`$ (see Paper I).

## 4 RESULTS

Here we present the results of our modeling. First we investigate the necessity for dark matter in the system by comparing two simple models (§4.1). Then we consider more generalized mass models (§4.2). Finally, we report constraints on the orbit anisotropies (§4.3).

### 4.1 Test for Dark Matter

We first fit the two simple mass models (SIS and constant $`M/L`$) to the stellar data \[$`I(R)`$, $`\widehat{\sigma }_\mathrm{p}(R)`$, $`h_4(R)`$, $`h_6(R)`$\] over a sequence of the mass parameters $`\sigma _0`$ and $`\mathrm{{\rm Y}}_B`$ (see Fig. 4). Table 1 lists the best-fit values for these parameters and the difference in log-likelihood between the two models. For the SIS model, we find $`\sigma _0=272_5^{+19}`$ km s-1; the fit of the best solution to the data is shown in Fig. 5. For the constant $`M/L`$ model, we find $`\mathrm{{\rm Y}}_B=8.1\pm 0.6`$. This model does not fit the data as well as the best SIS model does (1 $`\sigma `$ significance), providing weak evidence for a dark matter halo in this galaxy.

In order to see which constraints are important for indicating the presence of dark matter, we have tried models which lack either the $`h_6(R)`$ data or the large-radius $`\widehat{\sigma }_\mathrm{p}(R)`$ data, (Table 1). Dropping these constraints makes little difference in differentiating between mass models (though the mass normalization of the SIS model does change in the latter case). In particular, we note that the extra velocity dispersion data is unhelpful, owing to the flexibility of the orbital anisotropies to reproduce an arbitrary velocity dispersion profile. Only with additional LOSVD information at large radii can the stellar kinematical data by itself further constrain the mass distribution.

FIG. 4.— The log-likelihood of the best-fit singular isothermal model (left panel) and constant mass-to-light ratio model (right panel), as a function of the mass parameters $`\sigma _0`$ and $`\mathrm{{\rm Y}}_B`$. The filled squares show the fit to the stellar data, and the open squares the fit to the globular cluster data.

|  | Singular Isothermal |  | Constant $`M/L`$ | |
| --- | --- | --- | --- | --- |
| Constraints | $`\sigma _0`$ (km s-1) |  | $`\mathrm{{\rm Y}}_B`$ ($`M_{}/L_{,B}`$) | $`\mathrm{\Delta }\mathrm{ln}`$ |
| Stars: $`I`$, $`\widehat{\sigma }_{\mathrm{vdM}}`$, $`\widehat{\sigma }_{\mathrm{ST}}`$, $`h_4`$, $`h_6`$ | $`272_5^{+19}`$ |  | $`8.1\pm 0.6`$ | 0.8 (1.0 $`\sigma `$) |
| Stars: $`I`$, $`\widehat{\sigma }_{\mathrm{vdM}}`$, $`\widehat{\sigma }_{\mathrm{ST}}`$, $`h_4`$ | $`275_9^{+19}`$ |  | $`8.1_{0.6}^{+0.7}`$ | 0.6 (0.9 $`\sigma `$) |
| Stars: $`I`$, $`\widehat{\sigma }_{\mathrm{vdM}}`$, $`h_4`$, $`h_6`$ | $`254_{18}^{+22}`$ |  | $`8.3_{1.7}^{+1.0}`$ | 0.9 (1.0 $`\sigma `$) |
| Globular clusters | $`409_{24}^{+36}`$ |  | $`42.9_{3.1}^{+3.3}`$ | 12.9 (4.7 $`\sigma `$) |
| Stars + Globular clusters | $`301_{14}^{+5}`$ |  | $`9.3_{0.2}^{+0.1}`$ | 206.0 (20.1 $`\sigma `$) |

TABLE 1.— M87 model solutions, given different sets of observational constraints. The best-fit mass parameters ($`\sigma _0`$ and $`\mathrm{{\rm Y}}_B`$) are given for the SIS and constant $`M/L`$ models. The difference in log-likelihood between the two models is shown, along with the statistical significance at which the SIS model is preferred. The “Bayesian” probability for the SIS model is $`(1+e^{\mathrm{\Delta }\mathrm{ln}})^11`$, while that for the constant $`M/L`$ model is $`e^{\mathrm{\Delta }\mathrm{ln}}(1+e^{\mathrm{\Delta }\mathrm{ln}})^1e^{\mathrm{\Delta }\mathrm{ln}}`$ for $`\mathrm{\Delta }\mathrm{ln}>2`$.

FIG. 5.— Fits to the M87 stellar data (error bars) for the best singular isothermal model ($`\sigma _0=272`$ km s-1 : solid lines), the best constant mass-to-light ratio model ($`\mathrm{{\rm Y}}_B=8.1`$ : dotted lines), and the best NFW2 model (see §4.2; dashed lines). The data are the radial profiles of the surface brightness, the velocity dispersion, and the fourth- and sixth-order Gauss-Hermite velocity moments (from top to bottom). Note that $`R_{\mathrm{eff}}100^{\prime \prime }`$.

We fit the GC data to the same two mass models (see Fig. 4 and Table 1). For the SIS model, we find $`\sigma _0=409_{24}^{+36}`$ km s-1; the fit of the best SIS solution to the data is shown in Figs. 6–8. Note especially in Fig. 7 that some of the higher-order features of the observed LOSVDs are reproduced by the model. For the constant $`M/L`$ model, we find $`\mathrm{{\rm Y}}_B=43\pm 3`$. This model fits the data much worse than does the SIS model (5 $`\sigma `$ significance). Here, the constant $`M/L`$ model cannot reproduce the radially-rising GC velocity dispersion (see Fig. 6). Also, $`\mathrm{{\rm Y}}_B=43`$ is implausible for a standard stellar population, requiring an age $`17`$ Gyr and a metallicity $`[\mathrm{Fe}/\mathrm{H}]0.50`$ (Worthey 1994), so that even if mass traced light, the mass could not consist of a standard stellar population. Finally, note that $`\mathrm{{\rm Y}}_B=43\pm 3`$ is completely inconsistent with $`8.1\pm 0.6`$ as obtained via the stellar dynamics, unless there is a drastic radial gradient in the stellar population. We thus conclude that the constant $`M/L`$ model is definitively ruled out for M87, assuming standard gravitational dynamics. (Note that the 20 $`\sigma `$ significance given in Table 1 is formal only, and that in practice, a number of systematic uncertainties will reduce this significance — but it can still be said that the result holds to a high significance level.)

FIG. 6.— Fits to the M87 globular cluster data (error bars) for the best singular isothermal model ( $`\sigma _0=409`$ km s-1: solid lines), the best constant mass-to-light ratio model ($`\mathrm{{\rm Y}}_B=42.9`$: dotted lines), and the best NFW2 model (see §4.2; dashed lines). The data include the radial profiles of the surface density, the velocity dispersion, and the fourth-order Gauss-Hermite velocity moment (from top to bottom). Note that only the surface density is actually fit, while the other data are shown for comparative purposes only.

FIG. 7.— Line-of-sight velocity distributions for the M87 globular cluster system, in nine radial bins. The heavy solid lines show the LOSVDs of the best-fit singular isothermal solution ($`\sigma _0=409`$ km s-1), averaged in each bin. The light solid lines show the data (26 points in each bin). The dotted lines show simulated LOSVDs derived from the superposition of the data in each bin, symmetrized about $`v_\mathrm{p}=0`$. Compare Fig. 3.

FIG. 8.— Contour plot of line-of-sight velocity distribution ($`I^1dI/dv_\mathrm{p}`$) for the best-fit singular isothermal globular cluster system model ($`\sigma _0=409`$ km s-1), normalized by the surface brightness at each radius $`R`$. The last contour represents zero intensity. The data are shown as squares. Note that the model LOSVD is peaked at small radii ($`R<50^{\prime \prime }`$) and is flat and even double-peaked elsewhere.

### 4.2 Generalized Mass Models

With a constant $`M/L`$ model ruled out, the next question is what the actual radial distribution of mass is in the system. The simple SIS model is seen to be inaccurate: the $`\sigma _0`$ mass parameters derived separately from the stellar and GC constraints are inconsistent at the 99% confidence level (see Fig. 4 and Table 1). That the GC-derived mass is higher than the stellar-derived mass implies that the radial mass profile $`M(r)`$ rises more steeply than isothermal ($`r`$) in the outer parts. To make an initial estimation of $`M(r)`$, we use the spherical isotropic Jeans equations (see Sargent et al. 1978 equation 11 and compare CR Fig. 6), modeling the stellar and GC subsystems separately. The stellar constraints imply $`Mr^{0.9}`$ out to $`r150{}_{}{}^{\prime \prime }10`$ kpc, while the GC constraints imply $`Mr^{1.3}`$ outside this radius (see Fig. 9). These crude Jeans models thus suggest that it would be fruitful to further explore two-component mass models consisting of a constant $`M/L`$ galaxy plus a dark halo with $`M(r)`$ steeper than isothermal in the region of $`r100^{\prime \prime }`$$`500^{\prime \prime }`$ (7–36 kpc).

FIG. 9.— The radial mass profile of M87. The heavy solid lines show estimates from the Jeans equations using the stars and globular clusters separately. The long-dashed lines show confidence limits from the X-ray analysis of Nulsen & Böhringer (1995). The other lines show the best orbit models fit to the combined star and GC data, for the constant $`M/L`$ model (dotted), the SIS model (dashed), and the NFW2 model (light solid). The point with error bars shows the estimated mass of the central supermassive black hole (Macchetto et al. 1997).

We next use our orbit modeling methods to examine such a generalized mass model, where the form of the dark halo is based on the the predictions from simulations of cosmologically-motivated structure formation (e.g., Navarro, Frenk, & White 1996). This model has a total radial mass density profile of

$$\rho (r)=\mathrm{{\rm Y}}_B\nu _{}(r)+\frac{\rho _0r_s^3}{r(r+r_s)^2},$$

where $`r_s`$ is the break radius between the halo’s inner $`r^1`$ profile and outer $`r^3`$ profile. We construct mass models with a variety of the parameters $`(\mathrm{{\rm Y}}_B,\rho _0,r_s)`$, fitting the combined stellar and GC constraints, and report the results of our three best models in Table 2 (see also Figs. 5 and 6). Note that each distinct model begins with a set of kernels calculated for two fixed parameters $`(\mathrm{{\rm Y}}_B/\rho _0,r_s)`$, and may then have its overall mass normalization scaled linearly to best fit the data.

|  | Model Parameters | | | |  |
| --- | --- | --- | --- | --- | --- |
| Model name | $`\mathrm{{\rm Y}}_B`$ ($`M_{}/L_{,B}`$) | $`\sigma _0`$ (km s-1) | $`\rho _0r_s^3`$ ($`10^{12}M_{}`$) | $`r_s`$ (′′; kpc) | $`\mathrm{\Delta }\mathrm{ln}`$ |
| Constant $`M/L`$ | 9.3 |  |  |  | 208.3 (20.3 $`\sigma `$) |
| SIS |  | 301 |  |  | 3.7 (2.3 $`\sigma `$) |
| NFW1 | 5.6 |  | 0.9 | 500 ; 36 | 0.4 (0.8 $`\sigma `$) |
| NFW2 | 5.9 |  | 2.4 | 939.3 ; 68 |  |
| NFW3 | 6.5 |  | 4.8 | 1800 ; 131 | 0.6 (0.9 $`\sigma `$) |

TABLE 2.— M87 model solutions, given combined stellar and globular cluster constraints. The best-fit model parameters are given for each model (see text for description). The difference in log-likelihood between each model and the best overall model (NFW2) is shown, along with the statistical significance at which the model is dispreferred.

All these models (NFW1, NFW2, and NFW3) can produce better fits to the full data set than can the best SIS model (at the 98% significance level), and thus probably better approximate the overall mass distribution. The mass profiles of these models are consistent with $`M(r)`$ derived from the galaxy’s X-ray halo (Nulsen & Böhringer 1995), both in mass normalization and in density exponent (see Fig. 9). Both our NFW models and the X-ray model have density profiles that decline more gradually than isothermal in the galaxy’s outer parts (e.g., $`\rho r^{1.5}`$ at $`r200{}_{}{}^{\prime \prime }`$ 15 kpc), giving a circular velocity curve $`v_\mathrm{c}(r)`$ that continues rising at large radii (Fig. 10). Most elliptical galaxies have a constant or declining $`v_\mathrm{c}(r)`$ at these radii (Gerhard et al. 2001; the notable exception is the Fornax cD galaxy NGC 1399), so this is strong evidence that the dark matter at large radii is associated with the Virgo Cluster itself.

FIG. 10.— Circular velocity profile of M87, for different models. Lines styles are as in Fig. 9, with the addition of the NFW1 (dot-short dashed) and NFW3 (dot-long dashed) models. The vertical dashed line indicates the outer radius of the region constrained by the GC data.

A caveat about these new models is that the mass normalization from the stellar and GC data separately is again inconsistent (at the 1–2 $`\sigma `$ level). This may indicate that a still more accurate mass model $`M(r)`$ could be found to fit all the existing constraints (possibly one with a larger, less-concentrated cluster halo as in McLaughlin 1995), and/or that one or both of the dynamical subsystems (stars and GCs) would be better modeled without the assumption of spherical symmetry. Further exploration of more detailed models would be useful, but outside the scope of this project.

### 4.3 Orbit Anisotropies

We have not implemented any method of rigorously exploring the range of DF functional forms permitted by the data, so the following results are meant to be suggestive rather than conclusive. Fig. 11 illustrates the stellar orbital characteristics of the overall best-fit solutions, where the mean anisotropy is characterized by the parameter $`\beta 1v_\theta ^2/v_r^2`$. The NFW solutions have a somewhat radial DF in the region constrained by the kinematical data ($`r0`$.′′$`15`$$`170{}_{}{}^{\prime \prime }0.01`$–12 kpc). The downwards dive in $`\beta `$ near $`r=0`$.′′$`7`$ is due to the sharp turnover in $`\widehat{\sigma }_\mathrm{p}(R)`$ inside that radius. Real galaxies presumably have smooth $`\beta (r)`$ profiles, but with no smoothness imposed, our models are somewhat noisy because of the ill-conditioning of this inversion problem. Also, the region $`r<1^{\prime \prime }`$ is within the radius of influence of the central supermassive black hole (with mass $`M_{\mathrm{BH}}=3.2\pm 0.9\times 10^9M_{}`$; Macchetto et al. 1997), and thus has not been accurately treated by our models. In summary, our best models indicate mild radial orbit anisotropy ($`\beta =`$ 0.2–0.5) for the well-constrained stellar region of the galaxy ($`r1^{\prime \prime }`$$`30{}_{}{}^{\prime \prime }`$ 0.1–2 kpc) — a characteristic which has turned out to be quite typical for bright elliptical galaxies (see Gerhard et al. 2001 and references therein).

FIG. 11.— Stellar orbital anisotropy as a function of radius, for the overall best-fit solutions from five different mass models: constant $`M/L`$ (dotted line), SIS (solid line), NFW1 (dot-short dashed line), NFW2 (dashed line), and NFW3 (dot-long dashed line). The vertical dashed lines indicate the radial regions constrained by the stellar kinematical data.

Fig. 12 shows the orbital characteristics of the globular cluster system for the overall best-fit solutions. There is considerable variation between the best NFW solutions, so we cannot yet strongly constrain the GC anisotropy, except to note that it seems roughly isotropic at large radii ($`r100^{\prime \prime }`$$`500{}_{}{}^{\prime \prime }`$ 10–40 kpc). This seems to contradict the tangential anisotropy one would expect given the flat and double-peaked LOSVDs seen in the data and in the models that fit only the GC data (see §2.4 and Figs. 7 and 8); — illustrating further the need for a mass model that better fits both the stellar and GC constraints. Note that the solutions’ radial anisotropy at $`r<40{}_{}{}^{\prime \prime }<`$ 3 kpc should not be considered highly significant since there are only 8 GC velocities measured inside that projected radius.

FIG. 12.— Velocity anisotropy of globular cluster system as a function of radius, for different mass models (see Fig. 11 for line descriptions). The vertical dashed lines indicate the radial region constrained by the GC velocity data.

## 5 CONCLUSIONS

We have developed a very general spherical orbit modeling method that makes full use of all the information in a set of discrete velocity data, and have applied it to the giant elliptical galaxy M87 using all the available dynamical constraints (measurements of the spatial distribution and the kinematics of the stars and of the globular clusters). By comparing two simple mass models — a constant mass-to-light ratio galaxy and a singular isothermal sphere — we find that a constant $`M/L`$ is ruled out to a high level of significance.

We further implement some more generalized mass models in order to more accurately represent the galaxy as consisting of a constant $`M/L`$ stellar system plus a dark halo. We find improved solutions which have a dark halo density profile decreasing more slowly than $`r^2`$ at large radii (10–40 kpc). This very extended halo seems best interpreted as belonging to the Virgo Cluster itself, rather than to M87.

There is substantial room for improvement of our M87 models. First, measurements of higher-order velocity moments for the stars are available to only $`0.3R_{\mathrm{eff}}`$. Such measurements at larger radii would be quite helpful in better constraining the models, and can be most efficiently obtained by measuring discrete velocities of the galaxy’s planetary nebulae. Second, a larger range of mass models could be explored to find solutions that fit the stellar and GC data together more consistently. And third, future models could relax our assumption of spherical symmetry, providing a more accurate model of the system in its outer parts, where flattening is seen in both the stellar and GC systems.

The same techniques should also be applied to other, non-cD giant elliptical galaxies in clusters (such as NGC 4472), as well as to galaxies in the field. The different formational histories of these galaxies (via mergers, accretion, etc.) may result in significantly different mass distributions and dynamics. Furthermore, the colors, abundances, and spatial distributions of the GCs in M87 and other galaxies indicate that they comprise two or more distinct populations (e.g., Neilsen & Tsvetanov 1999; Gebhardt & Kissler-Patig 1999). It would be interesting to use our orbit modeling methods to look for any dynamical differences between these populations, which may shed light on their formation histories, for which several different scenarios have been proposed (e.g., Ashman & Zepf 1992; Forbes, Brodie, & Grillmair 1997; Harris, Harris, & McLaughlin 1998; Côté, Marzke, & West 1998). The available set of kinematical tracers observed in the outer parts of early-type galaxies is now increasing rapidly, and we are hopeful that a much clearer picture of galaxy halos is forthcoming.

We thank John Blakeslee, Judy Cohen, Dan Fabricant, Marijn Franx, Lars Hernquist, Konrad Kuijken, and the anonymous referee for helpful comments; Ken Sembach for providing his data in tabular form; John Huchra, Roeland van der Marel, and Dean McLaughlin for both; Arunav Kundu for providing his paper in advance of publication, and his data in tabular form. C.S.K. is supported by the Smithsonian Astrophysical Observatory.