text:
# Heavy Flavor Enhancement as a Signal of Color Deconfinement

\[

## Abstract

We argue that the color deconfinement in heavy ion collisions may lead to enhanced production of hadrons with open heavy flavor (charm or bottom). We estimate the upper bound of this enhancement.

\]

The production of open heavy flavor (HF) hadrons (charm and bottom) still remains a ’terra incognita’ of heavy ion physics: neither open charm nor open bottom yields have been measured so far. Open charm measurements are only planned in Pb+Pb collisions at the CERN SPS. The standard theoretical picture assumes that the average number of hadron pairs with open HF $`HF_{AB(b)}`$ created in a nucleus-nucleus ($`A+B`$) collision at given impact parameter $`b`$ is simply connected111Here we neglect (anti-)shadowing effects which are expected to be not very large at SPS and RHIC energies. with the probability to create a HF hadron pair in a nucleon-nucleon (N+N) collision $`HF_{NN}`$:

$$HF_{AB(b)}=N_{coll}^{AB}(b)HF_{NN}=N_{coll}^{AB}(b)\frac{\sigma _{NNHF+X}}{\sigma _{NN}^{inel}},$$
(1)
where $`N_{coll}^{AB}(b)`$ is the average number of primary nucleon collisions, which is determined by the geometry of the colliding nuclei, $`\sigma _{NNHF+X}`$ is the total cross section of the HF hadron pair production in $`N+N`$ collisions and $`\sigma _{NN}^{inel}`$ is the total inelastic cross section of $`N+N`$ interaction. (Note that in high energy collisions the HF production is dominated by the creation of hadrons with open HF. The HF quarkonia correspond to a tiny fraction of the total HF yield and can be safely neglected in our consideration.)

There are however some indirect indications that an essential deviation from the standard formula (1) may exist. Recent analysis of the dimuon spectrum measured in central Pb+Pb collisions at 158 A GeV by NA50 Collaboration  reveals a significant enhancement of the dilepton production in the intermediate mass region (1.5$`÷`$2.5 GeV) over the standard sources. The primary222Alternative explanations are also suggested .  interpretation attributes this observation to the enhanced production of open charm : about $`3`$ times above the direct extrapolation (1) from N+N data. Similar result has been recently obtained in the framework of the statistical coalescence model . This model connects the multiplicities of hadrons with open and hidden charm. It was found  that an enhancement of the open charm by the factor of about $`2÷4`$ over the direct extrapolation is needed to explain the data on the $`J/\psi `$ multiplicity. It was suggested in Ref. that this enhancement may appear due to the broadening of the phase space available for the open charm because of the presence of strongly interacting medium.

In the present letter we demonstrate that a deconfined medium (quark-gluon plasma (QGP) or its precursor) can make an essential influence on the hadronization of HF (anti-)quarks. This leads to an enhancement of the HF hadron production in A+B collisions in comparison to the direct extrapolation (1) from the N+N data. We restrict ourselves to a rough estimation of the upper bound of possible HF enhancement due to the color deconfined medium.

The process of production of a HF hadron pair can be subdivided into two stages: the hard production of a HF quark-antiquark pair ($`Q\overline{Q}`$) and its subsequent hadronization into observed particles. Therefore, there is an essential difference between HF hadron production and, e.g., hard dilepton production (the Drell-Yan process): created $`Q`$ and $`\overline{Q}`$ can and even have to interact with the surrounding quarks and gluons to be transformed into observed HF hadrons.

To get an intuitive picture of possible medium effects let us start from the open HF production in $`e^+e^{}`$ annihilation. The HF $`Q\overline{Q}`$ pair created at the first stage, hadronizes into observed particles. The hadronization has a nonperturbative nature. Its dynamics can be qualitatively understood in the framework of the string picture. When the distance between $`Q`$ and $`\overline{Q}`$ reaches the range of the confinement forces, a string connected these colored objects is formed. If the $`e^+e^{}`$ center-of-mass (c.m.) energy $`\sqrt{s}`$ (equal to the invariant mass of $`Q\overline{Q}`$ pair $`M_{Q\overline{Q}}`$) lies well above the corresponding HF meson threshold $`2m_M`$ (equal to $`2m_D`$ or $`2m_B`$ for $`c\overline{c}`$ and $`b\overline{b}`$ quarks, respectively), $`Q`$ and $`\overline{Q}`$ break the string into two (or more) peaces, so that the final state contains a HF hadron pair (and possibly a number of light hadrons). However, when the $`e^+e^{}`$ c.m. energy exceeds the heavy quark threshold ($`\sqrt{s}>2m_Q`$) but lies below the corresponding HF meson threshold $`2m_M`$ ($`\sqrt{s}<2m_M`$), the string cannot be broken and the open HF hadron pair can not be formed.

Let us imagine now the $`e^+e^{}`$ annihilation inside a deconfined medium. Due to the Debye screening, no string is formed between colored objects in this case. If the heavy $`Q`$ and $`\overline{Q}`$ are created, they can fly apart within the medium as if they were free particles. It does not matter whether their initial invariant mass $`M_{Q\overline{Q}}`$ exceeds the corresponding hadron threshold or not. The created $`Q\overline{Q}`$ pair will be able to form a HF hadron pair at the stage of QGP hadronization. This means that the $`e^+e^{}`$ annihilation inside the QGP would produce HF hadrons even if the collision energy is not sufficient for producing these hadrons in the vacuum.

In N+N or A+B collisions the HF $`Q\overline{Q}`$ pairs are produced due to hard parton interactions. The calculations in the leading order of the perturbative quantum chromodynamics (pQCD) show that a great fraction of $`Q\overline{Q}`$ pairs are created with invariant masses $`M_{Q\overline{Q}}`$ below the corresponding meson threshold $`2m_M`$ even at the largest RHIC energy. If this $`Q\overline{Q}`$ pair creation takes place in the deconfined medium, which is expected to be formed in high energy A+B collisions, the presence of such a medium makes possible a hadronization of these pairs. This should lead to an enhancement of the HF hadron production in A+B collisions in comparison to the standard result (1) obtained within the direct extrapolation of the N+N data.

There are of course essential differences between the open HF hadron production in the $`e^+e^{}`$ annihilation and in N+N or A+B collisions. Even in N+N collisions, when no deconfined medium is expected, the created $`Q\overline{Q}`$ pair can interact with the spectator partons and has therefore a chance to form HF hadron pair even if its primary invariant mass was insufficient for this process. Moreover, in contrast to the $`e^+e^{}`$ annihilation, the most of $`Q\overline{Q}`$ pairs are created in the color octet state and therefore they have to interact with the spectators to form a color neutral final state. Instead of breaking the string, the $`Q`$ and $`\overline{Q}`$ can form hadron states by means of coalescence with light spectator (anti-)quarks333 This mechanism is responsible, e.g., for the reaction $`ppp\mathrm{\Lambda }_c\overline{D}`$..

As no theoretical descriptions of this complicated process exist, we restrict ourselves to a rough estimation of the upper bound of possible HF hadron enhancement due to the color deconfined medium. We assume that

* In the case of N+N collisions, no subthreshold $`Q\overline{Q}`$ pairs contribute to the HF hadron production444 In other words, we assume that the interaction with the spectators does not change the energy of the $`Q\overline{Q}`$ pair and no coalescence with the spectator (anti-)quarks takes place..
* In the case of A+B collisions, provided that the deconfined medium is formed, all $`Q\overline{Q}`$ pairs hadronize into particles with the open HF555A small fraction of them form also the hidden heavy-flavor mesons, but it can be safely neglected..

The first assumption looks reasonable at low collision energies, whereas to justify the second one high energies are evidently preferable. This means that assuming validity of the both statements we overestimate the expected HF enhancement effect and, therefore, the above assumptions give its upper bound.

We make now the numerical estimates which follow from the above assumptions. The total cross section of heavy $`Q\overline{Q}`$ pair production by colliding nucleons is given by the formula (see e.g. Ref.)

$`\sigma _{NNQ\overline{Q}+X}(s)=`$ (2)
$`{\displaystyle \underset{(1,2)}{}}{\displaystyle _0^1}𝑑x_1{\displaystyle _0^1}𝑑x_2f_1(x_1,\mu _F)f_2(x_2,\mu _F)\widehat{\sigma }_{12Q\overline{Q}+\widehat{X}}(\widehat{s}),`$ (3)
where $`s`$ is the squared c.m. energy of the colliding nucleons, $`x_1`$ ($`x_2`$) is the fraction of the momentum of the first (second) nucleon carried by the parton $`1`$ ($`2`$), $`f_1`$ and $`f_2`$ are the fractional-momentum distribution functions or structure functions, $`\mu _F`$ is the factorization scale, $`\widehat{\sigma }_{12Q\overline{Q}}(\widehat{s})`$ is the cross section of heavy quark-antiquark pair production by interacting partons at squared center-of-mass energy $`\widehat{s}`$. For ultrarelativistic nucleons, $`\widehat{s}`$ is given by the formula $`\widehat{s}=x_1x_2s`$. The sum in the right hand side of Eq.(2) runs over all the pairs of parton types, that give nonzero contribution to the production cross section.

We restrict ourselves to the leading order of pQCD. In this case, two basic processes of heavy flavor creation have to be taken into account: the gluon fusion $`ggQ\overline{Q}`$ and the light quark-antiquark annihilation $`q\overline{q}Q\overline{Q}`$. So the sum in Eq.(2) includes $`(1,2)=(g,g),(\overline{q},q),(q,\overline{q})`$, where q in its turn runs over the light flavors $`q=u,d,s`$. The corresponding parton cross sections are given by the formulas :

$`\widehat{\sigma }_{ggQ\overline{Q}}(\widehat{s})={\displaystyle \frac{\pi \alpha ^2(\mu _R)}{3\widehat{s}}}`$ (4)
$`\times \left[\left(7+{\displaystyle \frac{31m_Q^2}{\widehat{s}}}\right){\displaystyle \frac{1}{4}}\chi +\left(1+{\displaystyle \frac{4m_Q^2}{\widehat{s}}}+{\displaystyle \frac{m_Q^4}{\widehat{s}^2}}\right)\mathrm{log}{\displaystyle \frac{1+\chi }{1\chi }}\right]`$ (5)
and

$$\widehat{\sigma }_{q\overline{q}Q\overline{Q}}(\widehat{s})=\frac{8\pi \alpha ^2(\mu _R)}{27\widehat{s}}\left(1+\frac{2m_Q^2}{\widehat{s}}\right)\chi ,$$
(6)
where $`\chi =\sqrt{14m_Q^2/\widehat{s}}`$, $`\mu _R`$ is the renormalization scale and $`m_Q`$ is the mass of the heavy quark. The masses of light quarks are neglected.

Eq.(2) can be rewritten in the form

$$\sigma _{NNQ\overline{Q}+X}(s)=_{(2m_Q)^2}^s𝑑\widehat{s}\frac{d\sigma _{NNQ\overline{Q}+X}}{d\widehat{s}},$$
(7)
where the differential cross section with respect to the squared invariant mass $`\widehat{s}=M_{Q\overline{Q}}^2`$ of the $`Q\overline{Q}`$ pair is given by the formula

$`{\displaystyle \frac{d\sigma _{NNQ\overline{Q}+X}}{d\widehat{s}}}=`$ (8)
$`{\displaystyle \frac{1}{s}}{\displaystyle \underset{(1,2)}{}}\widehat{\sigma }_{12Q\overline{Q}+\widehat{X}}(\widehat{s}){\displaystyle _{\widehat{s}/s1}^{1\widehat{s}/s}}𝑑x_L{\displaystyle \frac{f_1(x_1,\mu _F)f_2(x_2,\mu _F)}{x_1+x_2}},`$ (9)
where

$`x_1`$ $`=`$ $`\sqrt{\left({\displaystyle \frac{x_L}{2}}\right)^2+{\displaystyle \frac{\widehat{s}}{s}}}+{\displaystyle \frac{x_L}{2}}`$ (10)
$`x_2`$ $`=`$ $`\sqrt{\left({\displaystyle \frac{x_L}{2}}\right)^2+{\displaystyle \frac{\widehat{s}}{s}}}{\displaystyle \frac{x_L}{2}}.`$ (11)
The probability distributions of $`Q\overline{Q}`$ pairs with respect to $`\widehat{s}`$

$$w_{Q\overline{Q}}(\widehat{s};s)=\frac{d\sigma _{NNQ\overline{Q}+X}/d\widehat{s}}{\sigma _{NNQ\overline{Q}+X}(s)}$$
(12)
are shown in Fig. 1 and Fig. 2 for charm and bottom, respectively. The computation were done using the CERN library of parton distribution functions PDFLIB . The default set of structure functions MRS (G)  was chosen. The HF quark masses are fixed as $`m_c=1.25`$ GeV for charm and $`m_b=4.2`$ GeV for bottom, the c.m. energy of the colliding parton pair was used as the renormalization and factorization scales: $`\mu _F=\mu _R=\sqrt{\widehat{s}}`$.

We estimate now the upper bound of the HF enhancement in A+B collisions. We assume that in N+N collisions the HF $`Q\overline{Q}`$ pairs cannot hadronize, unless its c.m. energy exceeds the corresponding HF hadron threshold. Therefore, to calculate the total HF hadron production cross section we cut the integral in Eq.(7) at its lower bound by the corresponding meson threshold666 One has to cut the integral also at its upper bound, to respect the baryonic number conservation. Our calculations, however, show that this cut is important only at very low collision energies: $`\sqrt{s}2m_N+2m_M`$.:

$$\sigma _{NNHF+X}=_{(2m_M)^2}^{(\sqrt{s}2m_N)^2}𝑑\widehat{s}\frac{d\sigma _{NNQ\overline{Q}+X}}{d\widehat{s}},$$
(13)
where $`m_M`$ is the mass of the lightest meson containing corresponding HF quark ($`D`$-meson for the charm and $`B`$-meson for the bottom), $`m_N`$ is the nucleon mass.

In contrast, when two nucleons interact in the deconfined medium (as in high energy A+B collision), our assumption states that all $`Q\overline{Q}`$ pairs survive and form HF hadrons at the stage of the QGP hadronization. Therefore, the cross section $`\sigma _{NNHF+X}`$ in the formula (1) should be replaced by the cross section $`\sigma _{NNQ\overline{Q}+X}`$. Hence for the upper bound of the enhancement factor we use the formula

$$E^{\text{max}}(s)=\frac{\sigma _{NNQ\overline{Q}+X}(s)}{\sigma _{NNHF+X}(s)}.$$
(14)
The behavior of $`E^{\text{max}}(s)`$ for charm and bottom is shown in Fig.3. It is seen that the largest effect is expected at low energies. Therefore an experimental study of the effect should be done at the minimum energy, where the deconfinement medium is expected to be formed and the inclusive cross-section of HF production is large enough to make its measurement feasible.

The upper bound of open charm enhancement at SPS energy is by the factor of about $`5÷6`$. This means that the enhanced production of open charm hadrons by the factor $`2÷4`$ found in Ref. and Refs. can be explained by the influence of the deconfined medium.

We conclude that the deconfined medium, which is expected to be formed in nucleus-nucleus collisions can influence the process of hadronization of heavy quarks, this leads to the enhanced production of hadrons with open heavy flavors (charm and bottom). The rough estimation of the upper bound of the effect at SPS energies is found to be large enough to explain the indirect experimental data  and the phenomenological evaluations . We consider the enhancement of the heavy flavor yield as a possible signal of the color deconfinement. The direct measurement of open charm and open bottom in nucleus-nucleus collisions could be important for the confirmation of the quark-gluon plasma formation.

###### Acknowledgements.

The authors are thankful to K.A. Bugaev, M. Gaździcki, L. McLerran, I.N. Mishustin, G.C. Nayak, K. Redlich, E.V. Shuryak and H. Stöcker for fruitful discussions and their interest to the work. We appreciate criticism and useful comments of L. Frankfurt and L. Gerland. The financial support of DAAD, Germany is acknowledged. The research described in this publication was made possible in part by Award # UP1-2119 of the U.S. Civilian Research and Development Foundation for the Independent States of the Former Soviet Union (CRDF).


format: text

text:
# Discrete Solitons and Breathers with Dilute Bose-Einstein Condensates

## Abstract

We study the dynamical phase diagram of a dilute Bose-Einstein condensate (BEC) trapped in a periodic potential. The dynamics is governed by a discrete non-linear Schrödinger equation: intrinsically localized excitations, including discrete solitons and breathers, can be created even if the BEC’s interatomic potential is repulsive. Furthermore, we analyze the Anderson-Kasevich experiment \[Science 282, 1686 (1998)\], pointing out that mean field effects lead to a coherent destruction of the interwell Bloch oscillations.

Localization phenomena are ubiquitous in physics and in biology . Intrinsically localized excitations, as solitons (shape preserving) and breathers (characterized by internal oscillations), are important channels for energy transport in non-linear media, such as optical fibres and waveguides , polaronic materials  and biological molecules . Intense theoretical research is now focusing on the existence of solitons and breathers in a lattice (often named discrete solitons and discrete breathers (DSB) ), such as those appearing in quantum systems governed by a discrete non-linear Schrödinger equation (DNLSE) . Current approaches include the search for exact solutions in some limits ; effective (point) particle and variational approaches ; perturbation around the linearized case and, of course, numerical solutions . Although intensely studied, DSB have been experimentally observed only quite recently in superconducting ladders of Josephson junctions , in antiferromagnet systems , in optical waveguides  and in low dimensional materials .

The discrete solitons/breathers are characterized by a dynamical, self-maintained energy localization, due to both the discreteness and the non-linearity of the underlying equations of motion. The discreteness provides a band structure of the excitation spectrum, while the non-linearity allows for the tuning of the DSB energy outside the band. The finite energy gap guarantees the (meta-)stability of the DSB’s. These, obviously, have a different nature to the ”Anderson localizations”, created by impurities or imperfections of the lattice ; the incorporation of disorder into non-linear excitations, and the tracing out of its dynamical effects are also important theoretical problems.

Bright solitons can occur in spatially homogeneous, dilute Bose-Einstein condensates (BEC) with an attractive interatomic interaction (s-wave scattering length $`a<0`$) . Dark solitons, propagating density dips, have been predicted and experimentally observed in BEC’s with a repulsive interaction ($`a>0`$) . The dynamics of a BEC trapped in a spatially periodic potential , on the other hand, can be mapped, in the tight binding approximation, to a DNLSE. Bright discrete solitons and breathers can be created in DNLSE even with a repulsive BEC interatomic interaction.

In Ref.  a one-dimensional, vertical optical array was created by two counterpropagating laser beams. A weakly interacting Bose-Einstein condensate was trapped in $`30`$ wells, situated at the antinodes of the standing optical wave. Each well contained approximately $`1000`$ condensate atoms, with the peak densities matching a gaussian profile. The lowest Bloch band dynamics of this system maps on a DNLSE, obtained by discretizing the Gross-Pitaevskii equation governing the condensate dynamics in the periodic potential. Since the array is oriented vertically, the atoms undergo coherent Bloch oscillations, driven by the interwell gravitational potential. At the edge of the Brillouin zone, a fraction of atoms can Zener tunnel in the higher energy band which, in this specific case, is in the continuum. A coherent leakage from the trap was observed at each Bloch period. However, an increase of the on-site particle densities washed out the coherent signal. This has been understood as a drift of the relative phases between wells due to mean field effects .

In this Letter we study the DNLSE using a variational approach: we show that a (complex) gaussian variational ansatz (describing a soliton / breather profile, combined with a Lagrangian optimization) yields a coupled dynamics of the profile parameters. A variety of discrete solitons / breathers are investigated analytically and compared with numerical calculations. This comparison is surprisingly successful in describing even details of the quite complex dynamical and collisional behaviour. Stability phase diagrams for such states are obtained by inspection of the profile dynamics equations.

The full Bose condensate dynamics satisfies the Gross-Pitaevskii equation (GPE) :

$$i\mathrm{}\frac{\mathrm{\Phi }}{t}=\frac{\mathrm{}^2}{2m}^2\mathrm{\Phi }+[V_{ext}+g_0\mathrm{\Phi }^2]\mathrm{\Phi }$$
(1)
where $`V_{ext}`$ is the external potential and $`g_0=\frac{4\pi \mathrm{}^2a}{m}`$, with $`m`$ the atomic mass. For a tilted trap as in  the external potential is given by the sum of the gravitational potential and the laser field

$$V_{ext}(\stackrel{}{r})=mgz+U_L(x,y)\mathrm{sin}^2[2\pi z/\lambda ]$$
(2)
where $`\lambda `$ is the wavelength of the lasers (the spacing in the lattice is $`\lambda /2`$) and $`U_L(x,y)`$ is determined by the transverse intensity profile of the (nearly gaussian) laser beams. In  $`\lambda =850`$ nm and the $`1/e^2`$ radius of the transverse profile is $`80\mu `$m, an order of magnitude larger than the transverse radius of the condensate. The well depths scale linearly with the intensity of the beam. At the center of the beam the trap depths are $`1.4E_R`$ where $`E_R=\frac{\mathrm{}^2k^2}{2m}`$ is the recoil energy ($`k=\frac{2\pi }{\lambda }`$).

In the tight-binding approximation the condensate order parameter can be written as:

$$\mathrm{\Phi }(\stackrel{}{r},t)=\sqrt{N_T}\underset{n}{}\psi _n(t)\varphi (\stackrel{}{r}\stackrel{}{r}_n),$$
(3)
where $`N_T`$ is the total number of atoms and $`\varphi (\stackrel{}{r}\stackrel{}{r}_n)`$ is the condensate wave function localized in the trap $`n`$ with $`𝑑\stackrel{}{r}\varphi _n\varphi _{n+1}0`$, and $`𝑑\stackrel{}{r}\varphi _n^2=1`$. $`\psi _n=\sqrt{\rho _n(t)}e^{i\theta _n(t)}`$ is the $`n`$-th amplitude ($`\rho _n=N_n/N_T`$, where $`N_n`$ and $`\theta _n`$ are the number of particles and phases in the trap $`n`$). Replacing the ansatz (3) in (1) we find that the GPE reduces to a DNLSE:

$$i\frac{\psi _n}{t}=\frac{1}{2}(\psi _{n1}+\psi _{n+1})+(ϵ_n+\mathrm{\Lambda }\psi _n^2)\psi _n$$
(4)
where $`ϵ_n=\frac{1}{2K}𝑑\stackrel{}{r}\left[\frac{\mathrm{}^2}{2m}(\stackrel{}{}\varphi _n)^2+V_{ext}\varphi _n^2\right]`$, $`\mathrm{\Lambda }=\frac{g_0N_T}{2K}𝑑\stackrel{}{r}\varphi _n^4`$ and $`K𝑑\stackrel{}{r}\left[\frac{\mathrm{}^2}{2m}\stackrel{}{}\varphi _n\stackrel{}{}\varphi _{n+1}+\varphi _nV_{ext}\varphi _{n+1}\right]`$; the time has been rescaled as $`t\frac{\mathrm{}}{2K}t`$. Eq. (4) is the equation of motion $`\dot{\psi _n}=\frac{}{(i\psi _n^{})}`$, where $``$ is the Hamiltonian function

$$=\underset{n}{}[\frac{1}{2}(\psi _n\psi _{n+1}^{}+\psi _n^{}\psi _{n+1})+ϵ_n\psi _n^2+\frac{\mathrm{\Lambda }}{2}\psi _n^4]$$
(5)
with $`i\psi _n^{},\psi _n`$ canonically conjugate variables. Both the Hamiltonian $``$ and the norm $`_n\psi _n^2=1`$ are conserved.

To analyze the Anderson-Kasevich experiment, we study the dynamical evolution of a gaussian profile wave packet, which we parametrize as $`\psi _V^n(t)=\sqrt{k}\mathrm{exp}\left\{\frac{(n\xi )^2}{\gamma ^2}+ip(n\xi )+i\frac{\delta }{2}(n\xi )^2\right\}`$ where $`\xi (t)`$ and $`\gamma (t)`$ are, respectively, the center and the width (in the lattice units) of the density $`\rho _n=\psi _n^2`$, $`p(t)`$ and $`\delta (t)`$ their associated momenta and $`k(\gamma ,\xi )`$ a normalization factor. The wave packet dynamical evolution can be obtained by a variational principle from the Lagrangian $`=_ni\dot{\psi }_n\psi _n^{}`$, with the equations of motion for the variational parameters $`q_i(t)=\xi ,\gamma ,p,\delta `$ given by $`\frac{d}{dt}\frac{}{\dot{q}_i}=\frac{}{q_i}`$. After some algebra we obtain  $`=p\dot{\xi }\frac{\gamma ^2\dot{\delta }}{8}\frac{\mathrm{\Lambda }}{2\sqrt{\pi }\gamma }+\mathrm{cos}pe^\eta V(\gamma ,\xi )`$ where $`\eta =\frac{1}{2\gamma ^2}+\frac{\gamma ^2\delta ^2}{8}`$ and $`V=k\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑nϵ_ne^{\frac{2(n\xi )^2}{\gamma ^2}}`$. The variational equations of motion become:

$$\{\begin{array}{cc}\dot{p}=\frac{V}{\xi }\hfill & \\ \dot{\xi }=\mathrm{sin}pe^\eta \hfill & \end{array}$$
(6)
$$\{\begin{array}{cc}\dot{\delta }=\mathrm{cos}p\left(\frac{4}{\gamma ^4}\delta ^2\right)e^\eta +\frac{2\mathrm{\Lambda }}{\sqrt{\pi }\gamma ^3}\frac{4}{\gamma }\frac{V}{\gamma }\hfill & \\ \dot{\gamma }=\gamma \delta \mathrm{cos}pe^\eta \hfill & \end{array}$$
(7)
with the pairs $`\xi ,p`$ and $`\frac{\gamma ^2}{8},\delta `$ which are canonically conjugate dynamical variables with respect to the effective Hamiltonian

$$H=\frac{\mathrm{\Lambda }}{2\sqrt{\pi }\gamma }\mathrm{cos}pe^\eta +V(\gamma ,\xi ).$$
(8)
The wave packet group velocity is given by $`v_g\frac{H}{p}=\dot{\xi }=\mathrm{tan}p/m^{}`$ with an inverse effective mass $`(m^{})^1\frac{^2H}{p^2}=\mathrm{cos}pe^\eta `$. The quasi-momentum dependence of the effective mass allows a rich variety of dynamical regimes. Solitonic solutions with a positive non-linear parameter $`\mathrm{\Lambda }>0`$, for instance, are allowed by a negative effective mass. A regime with a diverging effective mass $`m^{}\mathrm{}`$ leads to a self-trapping of the wave packet.

In the following we will study the dynamical regimes of Eq. (4) in two particular cases. We first consider a tilted (washboard) periodic potential describing the vertical optical trap created in the Anderson-Kasevich experiment. In particular, we will show that non-linear effects tend to destroy the Bloch oscillations, consistently with the experimental results. The external potential for the tilted trap is given by Eq. (2) with $`ϵ_n=\omega n`$ where $`\omega =\frac{mg\lambda /2}{2K}`$. We find $`V=\omega \xi `$ and $`\dot{p}=\omega `$ .

Then we consider the case of an horizontal array, in which gravitation only provides a constant energy shift, with $`\dot{p}=0`$. We will classify four different regimes which include discrete breathers and solitons.

$`Tiltedtrap`$. It is well known that single atoms in a tilted washboard potential oscillate among sites at the Bloch frequency. This regime is described by Eqs. (6, 7) with $`\mathrm{\Lambda }=0`$ (corresponding to a negligible mean field condensate interaction). This is, precisely, the regime investigated in  in which a coherent output was observed. Indeed a variational estimate gives $`\mathrm{\Lambda }0.5`$ and $`\omega 2`$ (the scaled time is in units of $`\mathrm{}/2K=0.35ms`$). In this limit the equations of motion can be solved exactly giving $`p(t)=\omega t+p_0`$, $`\xi (t)=A[\mathrm{cos}p_0\mathrm{cos}(\omega tp_0)]`$, $`\gamma ^2(t)=4A^2\mathrm{log}A\omega \{\mathrm{cos}[2(\omega tp_0)]\mathrm{cos}2p_04\mathrm{sin}p_0\mathrm{sin}(\omega tp_0)4\mathrm{sin}^2p_0\}+\gamma _0^2`$ and $`\delta (t)=\frac{8A\mathrm{log}A\omega [\mathrm{sin}(\omega tp_0)+\mathrm{sin}p_0]}{\gamma ^2(t)}`$ where $`A=H_0/\omega \mathrm{cos}p_0`$, $`\delta (0)=\delta _0=0`$ and $`H_0`$ is the (conserved) initial energy. In the inset of Fig. 1 we show that the Bloch oscillations described by the variational ansatz (solid line) are in excellent agreement with the full numerical solution of the DNLSE (dashed line). The numerical average position is defined as $`_nn\psi _n^2`$.

The effect of non-linearity on the Bloch oscillations is dramatic. This has been studied experimentally in  by increasing the density in each well, and observing a degradation in the interference pattern. With $`\mathrm{\Lambda }0`$ we have:

$$\ddot{\xi }+\frac{\mathrm{\Lambda }\delta }{2\sqrt{\pi }\gamma }\dot{\xi }+\omega ^2\xi =\omega H_0\frac{\mathrm{\Lambda }\omega }{2\sqrt{\pi }\gamma }.$$
(9)
The Eq. (9) displays an effective damping term proportional to the velocity $`\dot{\xi }`$. We stress that the dynamics is fully Hamiltonian, and real dissipative processes are absent. For $`t\mathrm{}`$, $`\gamma `$ tends to a constant value $`\gamma _{fin}`$ and $`\delta \frac{2\mathrm{\Lambda }}{\sqrt{\pi }\gamma _{fin}^3}t`$, so the term $`\mathrm{\Lambda }\delta `$ has the correct positive sign (for large $`t`$). The apparent damping is the consequence of a diverging effective mass of the wave packet $`m^{}e^{\frac{\mathrm{\Lambda }^2}{2\pi \gamma _{fin}^4}t^2}`$, which stops the Bloch oscillations. This effect manifests, in the numerical analysis of Eq. (4), as a distortion of the on-site phases. In Fig.1 (solid line) we show the variational Bloch dynamics with a non-linear parameter $`\mathrm{\Lambda }=10`$ and initial values $`\xi _0=0`$, $`p_0=0`$ and $`\delta _0=0`$. The oscillation roughly decreases as $`\xi (t)A(1e^{\frac{\mathrm{\Lambda }^2t^2}{2\pi \gamma _{fin}^4}}\mathrm{cos}\omega t)`$. The dashed line shows the full numerical solution of the DNLSE, in good agreement with the analytical result (solid line). The discrepancy at $`t>10`$ is due to the breaking of the gaussian wave packet in the numerical simulation.

$`Untiltedtrap`$. The momentum $`p(t)=p_0`$ is conserved. We note that the equations of motion are invariant with respect to the replacement $`\mathrm{\Lambda }\mathrm{\Lambda }`$, $`p_0p_0+\pi `$ and $`tt`$. In Fig. 2 we report the dynamical phase diagram as a function of $`\mathrm{\Lambda }`$ versus $`\mathrm{cos}p`$, for $`\mathrm{\Lambda }>0`$ and with initial values $`\xi _0=0`$ and $`\delta _0=0`$. This phase diagram has been calculated analytically from Eqs. (6-8) and checked numerically. In the region $`\mathrm{cos}p>0`$ there are two distinct regimes. When $`H_0>0`$, $`\gamma (t)<\gamma _{max}`$: this is the self-trapped regime in which the boson wave packet remains localized around few sites. The self-localization is a genuine non-linear effect, characterized by a diverging effective mass. In particular, the self-trapped wave packet cannot translate along the array: this is a major difference with respect to soliton-like solutions. The limit values for $`t\mathrm{}`$ are $`\gamma \frac{\mathrm{\Lambda }}{2\sqrt{\pi }H_0}`$, $`\delta \mathrm{}`$ and $`\dot{\xi }0`$. We note that a non-linear self-trapping occurs also in a two-site problem .

A diffusive regime occurs when $`\mathrm{cos}p_0<H_00`$. In this case $`\gamma (t\mathrm{})\mathrm{}`$ and $`\dot{\xi }H_0/\mathrm{tan}p_0=const`$. The transition between the two regimes occurs at $`\mathrm{\Lambda }_c=2\sqrt{\pi }\gamma _0\mathrm{cos}p_0e^{1/2\gamma _0^2}`$. With $`\mathrm{\Lambda }>\mathrm{\Lambda }_c`$, the ratio between the initial value of the width $`\gamma _0`$ and the limit width $`\gamma _{max}(t\mathrm{})`$ is given by

$$\frac{\gamma _0}{\gamma _{max}}=\frac{\mathrm{\Lambda }\mathrm{\Lambda }_c}{\mathrm{\Lambda }}.$$
(10)
In Fig. 3 we plot $`\gamma _0/\gamma _{max}`$ vs. $`\mathrm{\Lambda }/\mathrm{\Lambda }_c`$ while in the inset we report the variational and numerical values of the width of the density vs. $`\mathrm{\Lambda }/\mathrm{\Lambda }_c`$ after a time $`\tau =10`$ (scaled units). The discrepancy at large $`\mathrm{\Lambda }/\mathrm{\Lambda }_c`$ is due to a slight deviation of the numerical density profile from a gaussian shape. We checked the stability of self-trapping, also considering different initial forms of the wave packet.

In the region $`\mathrm{cos}p<0`$ the self-trapping condition is given by $`H_0>\mathrm{cos}p_0`$: in this case $`\mathrm{\Lambda }_c=2\sqrt{\pi }\gamma _0\mathrm{cos}p_0(1e^{1/2\gamma _0^2})`$. A soliton solution can be determined by imposing $`\dot{\gamma }=\dot{\delta }=0`$. We find $`\mathrm{\Lambda }_{sol}=2\sqrt{\pi }\frac{\mathrm{cos}p_0}{\gamma _0}e^{1/2\gamma _0^2}`$. For $`\mathrm{\Lambda }=\mathrm{\Lambda }_{sol}`$ the center of the wave packet moves with a constant velocity $`\dot{\xi }`$ and its width remains constant in time. For $`\mathrm{\Lambda }_c<\mathrm{\Lambda }<\mathrm{\Lambda }_{sol}`$, $`\xi \mathrm{}`$ while $`\gamma (t)`$ oscillates, corresponding to a breather solution. The numerical solution of Eq. (4) confirms the variational predictions. In Fig. 4(a) we report the numerical density profile for $`\mathrm{\Lambda }`$ in the breather region with $`p_0=3\pi /4`$; in Fig. 4(b) the numerical density profile is drawn for the same $`\mathrm{\Lambda }`$ and $`p_0=\pi /4`$: the change of the sign of $`\mathrm{cos}p`$ gives a breather solution in the first case, and the spreading of the density in the second one. In Figs. 4(c) and 4(d) we show that for $`\mathrm{\Lambda }=\mathrm{\Lambda }_{sol}`$ the soliton after hitting a wall rebounds and regains its shape.

In conclusion, we have studied the dynamics of a dilute Bose condensate trapped in a periodic potential. We have analyzed the Anderson-Kasevich experiment, pointing out that mean field effects lead to a coherent destruction of the interwell Bloch oscillations. We have shown that intrinsically localized excitations, like discrete breathers and solitons, can exist also with a repulsive interatomic potential by studying analytically their dynamical phase diagram.

Discussions with S.R. Shenoy and P. Sodano are acknowledged. This work has been supported by the Cofinanziamento MURST.


text:
# A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy

## 1 Introduction

The center of our Galaxy has been extensively studied at wavelengths between the near infrared and radio regions in the electromagnetic spectrum, (Sandqvist & Genzel, 1992). The dynamical center of our Galaxy is considered to be coincident with a compact synchrotron radio source, Sgr A. Sgr A is associated with a black hole of mass of the order of 10$`{}_{}{}^{6}\mathrm{M}_{}^{}`$, (Bower & Backer, 1998). Kinematics in the neighborhood of this black hole would provide a direct estimate of the mass interior to the orbit. With this in mind, observational studies of gas and stellar motions in the very central regions of our Galaxy have been undertaken over the last few years.

This paper develops a new non-parametric method for mapping the distribution function and the underlying potential of a spherical isotropic stellar system, from discrete kinematic data, and applies it to published radial velocity and proper motion data in the Galactic center region to re-estimate the mass of the central black hole.

## 2 Observational studies

The notion that our own Galaxy could harbor a massive black hole at its center was alluded to by theories that suggested the presence of central massive black holes in galactic nuclei. Such a hypothesis regarding the Milky Way gained stronger footing with the availability of favorable observational evidence, starting with the report on Ne II $`12.8\mu `$m emitting high velocity ionized gas at the Galactic center by Wollman et al. (1976). An upper limit of $`4\times 10^6\mathrm{M}_{}`$ was imposed on the mass enclosed within a radius of $`0.8`$pc from the Galactic center. Similar observations of the $`12.8\mu `$m Ne II emission from the Galactic center were reported in subsequent years in Lacy et al. (1979) and Lacy et al. (1980). Later it was realized that the emission was from large-scale flows of ionized gas, (Lo & Claussen, 1983; Serabyn & Lacy, 1985; Genzel et al., 1985; Gusten et al., 1987). The streamer velocity and position could be fitted by orbital motions. Such fits indicated the presence of a mass of $`3\times 10^6\mathrm{M}_{}`$ within the inner $`0.5`$pc.

However stellar motions were still needed to get a more reliable estimate of the mass in the central regions in the Galaxy since gas dynamics can be dictated by factors other than gravity. Observations of the early type stars that form a smaller cluster found in the central $`10^{\prime \prime }`$ have also been made. Such an observational programme led Krabbe et al. (1995) to predict the existence of a central mass of $`24\times 10^6\mathrm{M}_{}`$ within the central $`0.15`$pc. Genzel et al. (1996) reported the results of a radial motion study of $`222`$ stars, within the central 1pc of the dynamical center of our Galaxy. This study predicted a central core of radius $`0.06`$pc and mass 2.2–3.2$`\times 10^6\mathrm{M}_{}`$. The work done in Eckart & Genzel (1997), was the first study of proper motions of stars in the Galactic center. Combining the radial motion data of Genzel et al. (1996) and the proper motion data from Eckart & Genzel (1997), Genzel and Eckart reached the conclusion that the central core of our Galaxy has a mass of $`2.45\pm 0.4\times 10^6\mathrm{M}_{}`$ and lies within $`0.015`$pc of Sgr A. Ghez et al. (1998) estimated the mass of the central dark core of the Milky Way to be $`2.6\pm 0.2\times 10^6\mathrm{M}_{}`$ and its volume to be $`10^6`$pc3. This conclusion was reached on the basis of the proper motion data of 90 stars in their sample. The latest report of stellar dynamics in the Galactic center is in Genzel et al. (2000). In this analysis, anisotropy of stellar motions is considered for the first time. The authors conclude that while isotropy is a good assumption, there are distinct groups of stars that display significant anisotropy in their motions. They estimate a central mass of $`2.63.3\times 10^6\mathrm{M}_{}`$ depending on the modeling of their data.

### 2.1 Radial and proper motion study by Genzel et al. and Eckart & Genzel

Genzel et al. (1996) observed radial velocities of 223 stars within $`0.1`$pc. Eckart & Genzel (1997) reported proper motions of stars in the central gravitational field of the Galaxy. The observational programme spanned over a total of four years, at a resolution of $`0.15^{\prime \prime }`$. It was found that the velocity dispersions along the three mutually orthogonal spatial axes were very similar, implying a very low degree of central anisotropy.

A strong Keplerian fall-off of velocity dispersion with projected radius \[$`\sigma (r_p)r_p^{0.5}`$\] was established, thus strongly implying the existence of a central black hole. This provided a strong physical motivation to model the center as a point mass and the velocity field as isotropic. The mass estimate of the central core of the Galaxy could then be found by fitting the parameters of the model to the radial and proper motion data. A simple virial estimate was supplemented by a Bahcall & Tremaine (1981) mass estimate for the Galactic core. For this latter scheme the point mass was considered to be embedded in an isotropic cluster of dispersions of the order of $`50`$km/s at infinity. Virial and Bahcall-Tremaine mass estimators gave similar values ($`2.12\pm 0.7\times 10^6\mathrm{M}_{}`$) for the mass within a projected radius of $`0.24`$pc. Using the Jeans equations with the radial and proper motion data from Genzel et al. (1996) and Eckart & Genzel (1997) implied the existence of a central mass of $`2.4\times 10^6\mathrm{M}_{}`$ with a total (systematic and fitting) uncertainty of about $`\pm 0.64\times 10^6\mathrm{M}_{}`$. This mass is considered to be concentrated within a radius of 0.015$`\pm `$0.02pc east of Sgr A.

### 2.2 Proper motion study by Ghez et al.

In this work Ghez and co-workers observed proper motions of the central cluster of the Galaxy. They followed the time evolution of the positions of $`90`$ stars over a period of two years at a resolution of $`0.05^{\prime \prime }`$. The radial velocities were measured for stars at the largest radii only, ($`0.1`$pc). These few values of the line-of-sight velocity dispersion, when compared to the projected radial and tangential velocity dispersions indicated a good degree of isotropy in the central regions. Modeling the center as a point mass, they estimated a central mass of $`(2.6\pm 0.2)\times 10^6\mathrm{M}_{}`$ interior to a radius of $`0.015`$pc.

### 2.3 Anisotropy incorporated analysis of stellar dynamics by Genzel et al.

In Genzel et al. (2000), the data set used includes line of sight velocity measurements reported earlier in Genzel et al. (1996), improved values of proper motions reported in Eckart & Genzel (1997), and proper motions reported in Ghez et al. (1998). This “homogenized best” data set includes proper motions of $`104`$ stars, with their errors as well as 227 radial velocities and errors in these line-of-sight (LOS) velocities. For 32 stars both sky and line-of-sight velocities are reported. The assumption of isotropy is concluded to be broadly good, though motions of individual groups of stars are spotted to show significant anisotropy, (examples of both radial and tangential anisotropy have been spotted). The whole sample of stars is divided into $`7`$ annuli and the projected radial and tangential velocity dispersions corresponding to each annulus are calculated. The global expectations of these projected velocity dispersions can be related to the anisotropy parameter, $`\beta `$. Thus for each radial annulus, the anisotropy parameter is obtained. The significance of the observed anisotropy is checked against data generated by Monte Carlo clusters defined by simple anisotropic distribution functions. Using the value of $`\beta `$ for each annulus, Leonard & Merritt (1989), Bahcall-Tremaine and virial mass estimates are made for the very central region, the respective results being $`2.9\pm 0.4\times 10^6\mathrm{M}_{}`$, $`3.1\pm 0.3\times 10^6\mathrm{M}_{}`$ and $`2.52.6\times 10^6\mathrm{M}_{}`$.

## 3 Motivation for a non-parametric scheme

The observational works discussed above obtain the mass in the central region of the Galaxy using either the Jeans equation or projected mass estimators like Bahcall-Tremaine, Leonard-Merritt, or the virial theorem. The shortcomings of these methods are as follows.

* The implementation of Jeans equation requires smooth approximations to the projected dispersions and the surface density, that can be subsequently related to the total dispersion and the spatial density via the Abel integral equation. This suffers from the problem of shot noise in radial bins, which tends to make the Abel inversion unstable (Genzel et al., 1996). In addition, the Jeans equations do not constrain the distribution function to be non-negative.
* The mass estimators are derived by taking moments of the Jeans equations, and relate $`r^nM(r)`$ (where $`M(r)`$ is the enclosed mass, $`n`$ is an integer, and the average is over the stellar density) to observed velocity dispersions. Being moments, they contain only limited information about the mass distribution. Moreover, for $`n=0`$ (Bahcall-Tremaine and Leonard-Merritt estimators) stars far from the center get more weight, which is undesirable when estimating a black hole mass. Also, kinematic data sets in practice sample only part of the potential, so some modeling or extrapolation is necessary.

A non-parametric method, which can reconstruct both density and distribution function without assuming functional forms for either, would avoid all these problems.

## 4 Practical algorithms used to solve the spherical inverse problem

Here we present a scheme for studying the spherical inverse problem under the assumption of isotropy. This algorithm is inspired by the work in Merritt (1993), but differs in some important ways: it uses both radial velocity and proper motion data rather than just the former, and uses the Metropolis algorithm as distinguished from maximum penalized likelihood. These differences are discussed in detail in the next section. The scheme reported in Merritt (1993) is a generic algorithm that does not involve the assumption of isotropy; a simple anisotropic distribution function that depends on both energy and the line-of-sight component of the angular momentum is considered. However, in this report, we solve the simpler case of the isotropic distribution function. (We buttress our analysis with a check of consistency of the data with isotropy). The algorithm in Merritt (1993), which uses discrete data, is discussed below.

### 4.1 A generic non-parametric method

Let a stellar system be represented by an equilibrium distribution function (DF), $`f`$, which is a function of energy $`E`$ and angular momentum, $`L`$. One way of projecting $`f`$ into the observable space is by first projecting it on the plane of the sky and then along the line of sight velocity axis in the phase space. The resulting projected DF will then be dependent on the plane of sky coordinates and also on the LOS velocity coordinate. This projected DF, say $`\nu _p`$, is related to $`f`$ via the integral equation

$$\nu _p(x,y,v_p)=𝑑z𝑑v_x𝑑v_xf(E,L),$$
(1)
where the LOS is along the $`z`$-direction while the plane of the sky is scanned by the $`x`$-$`y`$ plane. If the stellar system being studied exhibits spherical geometry, then the last equation becomes

$$\nu _p(r_p,v_p)=𝑑z𝑑v_x𝑑v_yf[v_r^2+v_t^2+2\mathrm{\Phi }(r),L].$$
(2)
Here $`r`$ is the spherical radius and $`r_p`$ the cylindrical radius on the plane of the sky, so that $`r_p^2=x^2+y^2=r^2z^2`$; $`v_r`$ and $`v_t`$ are components parallel and tangential to the radius $`𝐫`$; $`\mathrm{\Phi }(r)`$ is the gravitational potential of this system. To simplify matters we could choose the $`y`$-axis to lie completely in the plane containing the radius vector $`𝐫`$ and the LOS. This implies that $`v_r^2=(v_y\mathrm{sin}\theta )^2+(v_z\mathrm{cos}\theta )^2`$, $`\theta `$ being the polar angle, i.e., $`\mathrm{cos}\theta =z/r`$.

A non-parametric solution determines the potential and the distribution function by searching for a pair of functions $`f`$ and $`\mathrm{\Phi }`$ that provide the best-fit to the complete velocity data. The first step is therefore to identify a way by which a trial DF could be projected into observable space, according to eqn. 2. An easy way to do this would be to approximate $`f(E,L)`$ as a 2-D histogram.

Thus $`f`$ is approximated to be a constant over any integral cell defined around a pair of $`E`$ and $`L`$ values. If, of course, the assumption of isotropy is invoked, the integral space becomes one dimensional so that $`f`$ is a function of $`E`$ only. $`f`$ is then approximated to be a constant over any energy bin in this integral space. The contribution to the projected distribution function from each such energy bin is given by

$$\nu _p^{\mathrm{cell}}(r_p,v_p)=𝑑zA(r,r_p,v_p),$$
(3)
where $`A(r,r_p,v_p)=𝑑v_x𝑑v_y`$ is the area that the energy integral bin occupies in the $`(v_x,v_y)`$ space. The total projected distribution function is a simple sum over all the energy bins, i.e.,

$$\nu _p(r_p,v_p)=\underset{i}{}f_i\nu _{pi}^{\mathrm{cell}}(r_p,v_p).$$
(4)
To get the total distribution function $`f`$ from the projected DF, eqn. 4 has to be inverted, (eqn. 4 can be treated as a set of linear equations). If $`\nu _p(r_p,v_p)`$ is a continuous function of $`r_p`$ and $`v_p`$, the inversion is straightforward but if the data set consists of a large number of discrete data points, then the projected DF can be estimated by binning the data with respect to $`r_p`$ and $`v_p`$.

The probability that the observed kinematical data was drawn from $`f(E)`$ could be measured by the likelihood function ($`L`$) that can be defined as the product of all the projected distribution functions, each obtained for a pair of $`(r_p,v_p)`$

$$\mathrm{log}(L)=\underset{\mathrm{i}=1}{\overset{\mathrm{N}}{}}\nu _{pi},$$
(5)
where $`N`$ is the total number of pairs of phase space points in the observed data set and $`\nu _{pi}`$ is the projected distribution function corresponding to the $`i^{\mathrm{th}}`$ pair of apparent position and velocity in the data set. The optimization can be carried out in the presence of linear constraints on the total number of stars, presented as a penalty function. Thus in Merritt (1993), the maximum penalized likelihood was employed.

### 4.2 Unique features of our algorithm

Our method is essentially an isotropic version of the above, but there are a few important differences between the two algorithms, as discussed below.

1. We derive the potential from a discretized density distribution. The mass (including any dark matter) is assumed to be in spherical shells with density decreasing outwards but otherwise free-form.
2. As discussed in the last section, in Merritt (1993), a maximum penalized likelihood was suggested to find out the distribution function that the observed data was most likely to have been drawn from. When the likelihood is highly non-linear, it cannot be optimized by linear constraints or even by quadratic programming, (Merritt (1996)). In this case, the maxima in the likelihood can be found by searching for it using brute force (Merritt & Saha (1993)). However this method is rather unsatisfactory in terms of probability of success and required CPU time. We use the Metropolis algorithm, and attempt not only to recover the model with the maximum likelihood, but to generate an ensemble of models distributed according to the likelihood. This ensemble will of course be dominated by models close to the maximum likelihood model. The advantage of having such an ensemble is that we can estimate uncertainties simply by measuring the spread of values across the ensemble. Moreover, we can use discrete data without ever having to bin them.
3. Since we have both radial and proper motion data, we realized that we could project the distribution function into observable space in three ways. Firstly we could project the DF on the plane of the sky and then project it along the line-of-sight velocity axis in the phase space. This would give the projected DF, $`\nu _p`$ which will have dependence on the plane of sky coordinates and the line of sight velocity coordinate. We could also project the DF first on the plane of the sky and then on the plane which is perpendicular to the line of sight velocity axis in the velocity space. The resulting projected distribution function is then dependent on the plane of sky coordinates and also on the velocity coordinates perpendicular to the line of sight velocity coordinate. We name such a projected distribution function $`\nu _{}`$. There is a small sample of stars for which we know both radial and transverse velocities. This allows us to project the distribution function along each of the three velocity axes of the phase space, subsequent to the projection on the plane of the sky. The resulting projected distribution function will be dependent on the apparent position on the sky and the radial and tangential velocities. We call such a projected distribution function $`\nu _{3\mathrm{D}}`$.
4. We incorporate the effects of the errors in the velocity measurements by convolving the calculated projected distribution functions with error distributions that are assumed to be Gaussian.
5. It is rightly argued that finding the data to be inconsistent with anisotropy is a surer check of the assumption of an isotropic distribution function than checking for the consistency of the data with isotropy. But the former is a much harder test to execute. Hence we carry out our analysis by starting with the assumption of isotropy and at the end of analysis, perform a goodness-of-fit test to check for the validity of the assumption for the observed data sets.

### 4.3 The Metropolis algorithm

The algorithm used in our code to identify the distribution of likelihoods is the Metropolis algorithm, (Binney et al., 1992). The likelihood function $`L`$ is a function of the potential and of the distribution function histograms, i.e., $`L`$ is a function in $`n`$-dimensional space, where $`n`$ is the number of energy bins in the DF histogram times the number of mass shells. Let us represent the likelihood as $`L=L(𝐱)`$ where $`𝐱(\rho ,f)`$.

The code samples the function $`L(𝐱)`$ through a series of iterations $`L(𝐱_n)`$. Two consecutive iterations are related by a transition probability $`p(𝐱_n𝐱_{n+1})`$. The transition probabilities are chosen as

$$p(𝐱𝐱^{})=\mathrm{max}(\frac{L(𝐱^{})}{L(𝐱)},1),$$
(6)
although in principle any choice satisfying detailed balance

$$L(𝐱^{})p(𝐱^{}𝐱)=L(𝐱)p(𝐱𝐱^{})$$
(7)
would serve.

The algorithm is schematically represented in Figure A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy. At each iteration, the code tries a new $`𝐱`$, i.e., new $`\rho `$ and $`f`$. These histograms are always non-negative and monotonic, and with unit normalization for $`f`$, but otherwise arbitrary. The trial $`𝐱`$ is accepted or rejected with the transition probability; in case of rejection the code keeps the old $`𝐱`$ for the next iteration.

In its early phase Metropolis searches out the region around the maxima in the likelihood and then the more “equilibrium” phase sets in, when Metropolis wanders about in the region of the multi-dimensional energy-density space, where the likelihood is close to maximal. The extent of this wandering provides the uncertainties in the mass estimates.

As the iterations proceed, the likelihood is monitored. Initially, there is trend towards increasing likelihood. When this trend disappears, Metropolis is in the equilibrium phase. In practice we waited till iteration $`k`$ such that $`L(𝐱_k)<L(𝐱_{k799})`$, and then took $`\{𝐱_{k799},\mathrm{},𝐱_k\}`$ as our likelihood-weighted sample of $`\rho `$ and $`f`$. We recorded $`\{L(𝐱_{k799}),\mathrm{},L(𝐱_k)\}`$. The mass in the innermost radial bin, ($`M_{\mathrm{in}}`$) was calculated at iteration numbers $`k799`$, …, $`k`$, from the likelihood-weighted $`\rho `$ obtained at the end of each iteration. Thus a sample of these 800 values of $`L`$ and $`M_{\mathrm{in}}`$ were recorded.

For the test problem (Section 6 below) we compared the likelihood-weighted $`f`$ and $`\rho `$ thus generated, with $`f`$ and $`\rho `$ corresponding to a known model.

For real data, we added an extra precaution. We generated samples of 800 values of $`L`$ and $`M_{\mathrm{in}}`$ from 10 different initial configurations. To each sample we applied a goodness-of-fit test (see below); the samples that passed that test were concatenated into a master sample. The item of interest from this master sample is of course the mass in the innermost bin, histograms of which are displayed later in this paper. The median of such a histogram is quoted as the mass estimate for the black hole in our Galaxy while the 68% confidence limits are given by the $`16\%`$ and $`84\%`$ levels.

The Metropolis algorithm must not be confused with the related technique of simulated annealing. Simulated annealing could be used to locate the global maximum of $`L(𝐱)`$, but that is insufficient for our purposes because we are interested not just in the most-probable black hole mass but in its uncertainty.

### 4.4 Goodness of fit

To check our assumptions (especially isotropy of the velocity distribution) and the fitting procedure, it is necessary to have a goodness-of-fit test.

The first step in a goodness-of-fit test is to choose a statistic, or measure of goodness-of-fit. There are many possibilities, but the obvious one is the likelihood $`L`$.

The second step is to compute the statistic for the real data, and compare its value with the expected distribution of that statistic from the model. If 99% of synthetic data that is drawn from the model, fit better than the actual data then the model is rejected at 99% confidence, and so on.

To implement such a test, we took $`f`$ from the last iteration of Metropolis and generated 100 synthetic kinematical data sets of the same size as the real data set, and evaluated the statistic (i.e., $`L`$) for each of these. We then computed the fraction of the synthetic data sets that fit less well than the real data. We call this fraction $`P`$,111This is sometimes called the $`p`$-value of a statistic. and quote it as a percentage. Low $`P`$ implies that the fit must be rejected.

If we had been doing least-squares, the above prescription would reduce to a $`\chi ^2`$ test, since $`\chi ^2=\frac{1}{2}\mathrm{ln}L`$ in least-squares. The expected distribution of $`\chi ^2`$ is model-independent (as is the expected distribution of the KS statistic), so there is no need to generate synthetic data when using it.

For acceptable fits, $`P`$ is not uniformly distributed between 0 and 100%, but tends to be on the high side. This is because the fitting procedure has already in some sense optimized $`P`$. In least-squares, this effect is corrected for by subtracting the number of fitted parameters from $`\chi ^2`$; we do not know if this correction can be generalized.

The value of $`P`$ has been pictorially depicted in Figures A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy and A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy; for a given initial model, $`P`$ is measured and plotted in a stacked bar diagram, with the three different columns corresponding to the three different observed kinematical data-sets used in the work.

## 5 Implementation details of our algorithm

### 5.1 Calculation of the potential

The mass distribution is a series of concentric spherical shells each with different densities. The potential

$$\mathrm{\Phi }=\frac{1}{r}_0^r\rho (r^{})r_{}^{}{}_{}{}^{2}𝑑r^{}+_r^{\mathrm{}}𝑑r^{}\rho (r^{})r^{}.$$
(8)
can then be expressed as follows. If $`r`$ falls within the $`m`$-th shell (of $`n`$ shells in all) then

$$\mathrm{\Phi }=\frac{1}{r}\underset{\mathrm{i}=1}{\overset{\mathrm{m}1}{}}\rho _i(i^2i+\frac{1}{3})\delta ^3+\underset{\mathrm{i}=\mathrm{m}+1}{\overset{n}{}}\rho _i(i\frac{1}{2})\delta ^2+\rho _m\left[\frac{r^2}{6}+\frac{m^2\delta ^2}{2}\frac{(m1)^3\delta ^3}{3r}\right],$$
(9)
and if $`r`$ is outside all the shells then

$$\mathrm{\Phi }=\frac{1}{r}\underset{\mathrm{i}=1}{\overset{n}{}}\rho _i(i^2i+\frac{1}{3})\delta ^3.$$
(10)
Here $`\delta `$ is the thickness of each shell.

We constrain $`\rho `$ to be non-increasing as $`r`$ increases, but apart from that the density shells can take any positive values.

It needs to be emphasized here that we seek the mass of the black hole in the Galactic center as the mass that is included inside a pre-fixed radius, $`r_{\mathrm{in}}`$. This is the projected radius at which the innermost velocity data measurement is reported, in the used data sets. In the radial velocity measurement by Genzel et al. (1996), $`r_{\mathrm{in}}0.046`$pc while in the proper motion work by Ghez et al. (1998), it is about 0.0044pc. At $`r<r_{\mathrm{in}}`$, the code will spread the mass uniformly. Hence the best that we can do with the available data is to ask for the mass inner to the radius at which the innermost measurement has been recorded. As pointed out in Ghez et al. (1998), the mass estimate could be totally attributed to a central black hole only if resolution improves to the extent that $`r_{\mathrm{in}}`$ corresponds to the Schwarzschild radius of the black hole. This is merely wishful at this stage as the minimum radius at which a measurement is recorded in Ghez et al. (1998) is about 40,000 times less than the Schwarzschild radius of a $`2.6\times 10^6\mathrm{M}_{}`$ black hole.

### 5.2 The projected distribution functions

In the coordinate system that we adopt, the plane of the sky is scanned by the $`(x,y)`$ plane while the line-of-sight is along the $`z`$ direction.

The projected distribution function $`\nu _p`$ is given as

$$\nu _p(r_p,v_z)=𝑑zf(E)𝑑v_x𝑑v_y.$$
(11)
As discussed before, $`\nu _{}`$ is in general dependent on the apparent position and the velocity coordinates that are perpendicular to the line of sight velocity axis in velocity space, i.e., on $`v_x\&v_y`$. If we assume axisymmetry in velocity space, on the plane perpendicular to the $`𝐯_𝐩`$ axis, then we can represent $`\nu _{}`$ as the equation below, (where $`v_\mu =\sqrt{v_x^2+v_y^2}`$),

$$\nu _{}(r_p,\sqrt{v_x^2+v_y^2})=𝑑zf(E)𝑑v_z,$$
(12)
and $`\nu _{3\mathrm{D}}`$ is given as

$$\nu _{3\mathrm{D}}(r_p,v_\mu ,v_z)=𝑑zf(E).$$
(13)
The bounds on individual phase space coordinates that correspond to each such energy bin are identified. These bounds are then used as limits for the projection integrals that relate the projected DF for the $`i^{\mathrm{th}}`$ energy bin to the constant value (say $`f_i`$) of the total DF over this energy bin. It must be noted that the bounds on the phase space variables may depend on how the projection is executed.

Let the $`i^{\mathrm{th}}`$ energy bin run from energy $`=E_{1i}`$ to $`E_{2i}`$. Since $`r^2=r_p^2+z^2`$, the minimum value that $`r`$ could take over the current energy bin, (for a given value of the observed apparent position) is $`r_p^2`$, obtained by setting $`z`$ to zero. Let the maximum value that $`r`$ could have in this energy bin be $`r_{2i}`$. Now $`r`$ is maximum for a given energy when the kinetic energy is a minimum. For the projection of the distribution function along the vector $`𝐯_𝐩`$, $`v_p`$ is a measured quantity. Hence the minimum kinetic energy is obtained by setting $`v_\mu `$ to zero. This implies that when $`f`$ is projected along the LOS velocity, the upper limit on $`r`$ for a specified energy bin is: $`r_{2i}`$, where $`r_{2i}`$ is the root of

$$E_{2i}=\frac{1}{2}v_p^2+\mathrm{\Phi }(r_{2i}).$$
(14)
Similarly when $`f`$ is projected perpendicular to the LOS velocity, $`v_\mu `$ is a measured quantity. Therefore, minimizing kinetic energy implies setting $`v_p`$ to zero. Thus $`r_{2i}`$ in this case is defined as

$$E_{2i}=\frac{1}{2}v_\mu ^2+\mathrm{\Phi }(r_{2i}).$$
(15)
The roots are found in our algorithm with one of the standard root-finding routines.

For a given $`r`$, the speed $`v=\sqrt{v_p^2+v_\mu ^2}`$ is extremal when the energy is either a maximum or a minimum. The limits on speed $`v`$ ($`v_1`$ being the lower and $`v_2`$ the upper limit) are such that

$$\frac{1}{2}v(r)_{ji}^2=E_{ji}\mathrm{\Phi }(r),\text{where }j=1,2.$$
(16)
Using the transformations

$$r_p^2=r^2z^2,v_\mu ^2=v^2v_p^2,$$
(17)
we can write the integrals for the projected DFs over the current energy bin as

$`\nu _i(r_p,v_\mu )=f(E_i){\displaystyle _{r_p}^{r_{2i}(v_\mu )}}{\displaystyle \frac{r[\sqrt{v_{2i}^2(r)v_\mu ^2}\sqrt{v_{1i}^2(r)v_\mu ^2}]}{\sqrt{r^2r_p^2}}}𝑑r`$ (18)
$`\nu _{ip}(r_p,v_p)=f(E_i){\displaystyle _{r_p}^{r_{2i}(v_p)}}{\displaystyle \frac{r[v_{2i}^2(r)v_{1i}^2(r)]}{\sqrt{r^2r_p^2}}}𝑑r`$ (19)
$`\nu _{i3}(r_p,v_p,v_\mu )=f(E_i){\displaystyle _{r_p}^{r_{2i}(v_p,v_\mu )}}{\displaystyle \frac{r}{\sqrt{r^2r_p^2}}}𝑑r`$ (20)
To get the projected DFs, we need to sum over all the energy bins. (We constrain $`f`$ to be non-increasing as $`E`$ increases, i.e., $`f`$ is greatest for the most tightly bound orbits.) In this whole scheme $`r_p\&v_p`$ and $`r_p\&v_\mu `$ are read from the supplied data set. The sky is seen at the different data points, each data point being independent. Thus the probability that the positions and velocities in the data set have been drawn from the currently chosen distribution function will be measured by the product of the values of the projected DF obtained with the apparent position and velocity values from the whole data set. This product is then defined to be our likelihood function, $`L`$. Thus

$$\mathrm{log}L=\underset{\mathrm{data}\mathrm{pts}}{}\mathrm{log}\nu _k\text{where }k\text{ “}\text{” or “}p\text{” or “3D”}\text{}$$
(21)
The aim of this scheme is to start with an arbitrary pair of density and distribution function histograms and sample according to the likelihood function.

## 6 Tests of the Method

We tested our code on a Plummer model, with simulated data sets of similar size to the real data sets.

From the Plummer distribution function at the Plummer potential (Binney & Tremaine, 1987), we generated a data set containing 100 proper motion values, another containing 200 radial velocity values and a third with 30 values of both radial and transverse velocities. These numbers correspond roughly correspond to the populations in the observed data sets that we use.

The Plummer model is the polytrope with a polytropic index of 5. The distribution function for a general polytrope with index $`n`$ is given as:-

$$f=F(E)^{(n\frac{3}{2})}.$$
(22)
In our work we decided to set the total mass in the Plummer model to unity.

The aim is to start with an arbitrarily chosen pair of DF and density histograms, use them in the code that was described above, along with the data drawn from the Plummer distribution function. The forms of the density and distribution function profiles corresponding to the convergence of the Metropolis algorithm maximum value of the likelihood function were then compared with the Plummer density and distribution function forms.

The comparison between the final distribution function and density profiles thrown up by our code and the actual Plummer density and distribution function profiles are presented below. The initial distribution function and density profile for which Figure A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy was obtained were chosen as

$$f=F(E)^6,\rho =(1+r)^{0.9}.$$
(23)
Figure A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy shows that the comparison is quite good between the actual Plummer model and the profiles suggested by our algorithm. Encouraging comparisons were obtained even when the initial density profiles were of the Dehnen model type and the initial distribution function had a power-law dependence on $`1+(E)^2`$.

## 7 Results

### 7.1 Initial models used with observed data

The initial density profiles that we chose to input into Metropolis, while working with the observed data sets, correspond to the Dehnen models,

$$\rho =\frac{(1+\alpha )}{4\pi }r^{\alpha 2}(1+r)^{\alpha 2}$$
(24)
in the notation of Saha (1993). The density profiles in eqn. 24 can represent different mass distributions depending on the value of the parameter $`\alpha `$. For example, $`\alpha =1`$ implies a point mass while $`\alpha =0`$ represents a Jaffe model—see Binney & Merrifield (1998) for details. In our work $`\alpha `$ was varied in the range $`[0.1,1.0]`$.

As mentioned before, the distribution function in this work has been modeled as a function of energy only. The initial form of the distribution function was chosen to be a power law.

$$f(E)^\beta $$
(25)
with $`\beta `$ either 1.0 or 2.0.

### 7.2 Distribution functions and density profiles

Figure A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy shows the density profiles and distribution functions generated by Metropolis to fit the three data sets. This particular run started with a relatively steep DF and density. We find that $`f`$ does not vary much over the Metropolis sample—there is a strong tendency for $`f`$ to stay close to the initial state. For this reason we decided to collate the results from 10 initial states.

### 7.3 Working with proper motion data

Our toy data obtained from the Plummer model gave us confidence to use observed proper motion data presented in Ghez et al. (1998) to obtain an estimate of the mass of the central black hole in our Galaxy.

The proper motion values cited in Ghez et al. (1998) are accompanied by observational errors. While a series of runs were carried out, without taking into account the reported errors in the transverse velocities, another series of runs were done in which the errors were convolved with the projected distribution functions obtained with the corresponding pair of transverse velocity and apparent position. The error distribution was assumed to be a Gaussian. The initial models for both types of runs were the same, (Sec.7.1).

The histogram of the likelihood weighted estimate of the mass within $`r_{\mathrm{in}}=0.0044`$pc is shown in Figure A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy. The histograms were constructed for the sample of mass estimates in the way described in Sec.4.3, with those models that passed our goodness-of-fit test (Sec.4.4). For both the error-convolved and error-ignored runs, the initial models with the less steep distribution functions, passed the test satisfactorily while the initial DFs $`(E)^2`$ gave very poor values of the parameter $`P`$ defined in Sec.4.4. The dependence of $`P`$ on the initial models in the runs done with the proper motion data are represented in Figure A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy and Figure A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy.

From the error-convolved distribution of mass estimates shown in Figure A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy, we find a median of $`2.0\times 10^6\mathrm{M}_{}`$. At the lower 16% of the distribution, the mass is $`1.3\times 10^6\mathrm{M}_{}`$ while at the upper 84% the mass is $`2.7\times 10^6\mathrm{M}_{}`$.

### 7.4 Working with the radial velocity data

We have used $`204`$ of the radial velocity data and the errors in the same, reported in Genzel et al. (2000) to obtain another mass estimate of the black hole at the center of the Galaxy. For this data set, $`r_{\mathrm{in}}`$ is 0.046pc.

The initial distribution function and density profiles used in conjunction with the radial velocity data were as mentioned in eqn. 25 and eqn. 24 respectively. For the runs in which observational errors were incorporated, the distribution functions were convolved with the errors, the error distribution being assumed to be Gaussian. The histograms representing the values of likelihood weighted mass, are shown in Fig A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy. The histograms were obtained using the same scheme described in Sec.4.3. The histograms were constructed with only those initial models, which passed the goodness-of-fit test (Sec.4.4) satisfactorily. All the considered initial models pass the goodness-of-fit test quite well when observational errors are considered, but only half of the initial models do so when errors are ignored. Hence the sample of mass estimates for the error-convolved runs is much bigger than the same when errors are ignored. This also explains the longer tail of the error-convolved histogram.

We find that in the case of the error-convolved runs done with the radial velocity data set, $`P`$ (defined in Sec.4.4) is smaller for the runs done with the steeper initial distribution function ($`(E)^2`$) than with the flatter ones ($`(E)`$), but is never zero. The variation in $`P`$ with initial model is shown in Figure A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy and Figure A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy.

From the error-convolved histogram in Figure A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy, we obtain a median of $`2.2\times 10^6\mathrm{M}_{}`$ while the mass estimates at the 16% and the 84% levels are $`1.2\times 10^6\mathrm{M}_{}`$ and $`3.8\times 10^6\mathrm{M}_{}`$ respectively.

### 7.5 Working with the 3-D motion data

In Genzel et al. (2000), $`32`$ stars were reported with their proper motion and radial velocities as well as the errors in the same. For this data set, $`r_{\mathrm{in}}`$ was 0.046pc. Histograms were constructed with the same initial models as mentioned in Sec.7.1 for runs done both by convolving the projected distribution functions with errors as well as by ignoring the observational errors. These histograms can be seen in Figure A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy. All the considered initial models were found to pass the goodness-of-fit test, (Sec.4.4) for both types of runs, (with and without errors).

From the error-convolved histogram in Figure A Non-Parametric Estimate of the Mass of the Central Black Hole in the Galaxy, we obtain a median of $`1.8\times 10^6\mathrm{M}_{}`$ while the mass estimates at the 16% and the 84% levels are $`1.5\times 10^6\mathrm{M}_{}`$ and $`2.2\times 10^6\mathrm{M}_{}`$ respectively.

### 7.6 Results explained

Some of the results presented above have been qualitatively discussed below.

The distribution of the mass estimates has a non-zero width for two reasons.

* the mass estimate sample is recorded while Metropolis is wandering around the maximal region of the likelihood function during the “equilibrium phase” though it makes occasional journeys to lower likelihoods too.
* the mass estimate sample is composed of results obtained from different initial models. Though the initial model is not expected to dictate the convergence to the true distribution function (from which the observed data is drawn), the limited number of data points does limit the independence of the results from the initial model.

The results obtained with any initial model could fail the goodness-of-fit test if there is any systematic error (or bias) in the code (or the data) or if indeed the observed data is not consistent with isotropy or if Metropolis was unable to seek out the optimal fit to the observed data. If the test is passed, it implies that the used kinematical data agrees with isotropy. As described above, we carried out this test with a variety of initial distribution function and density profiles. We find that while the test is passed by some of the models, (usually the models characterized by the relatively flatter initial DF) when the radial or transverse velocity data is used, all the models easily pass the goodness-of-fit test when the 3-D motion information is used. If indeed there was anything systematically wrong, the analysis done with all the initial models would give similar results during the check for isotropy. This leads us to conclude that the reason for the failure of any initial model to pass the goodness-of-fit test is not a systematic error.

The error-convolved histograms usually imply a smaller mass for the black hole, for any given data set. The effect of errors in the velocity measurements is equivalent to introduction of noise, which implies greater energy for a star at a given phase space location than when the noise is absent. This in turn could (if uncorrected) imply a higher enclosed mass. When the error distribution is convolved with the projected distribution function, we compensate for this spurious enhancement in mass. This effect of noise in the data is conjectured to be increasingly overwhelmed with increasing number of used observational data points. Hence we would expect in our results that the broken histogram, (representing the distribution of mass estimates in absence of the error convolution) is furthest from the solid histogram, (representing the distribution of mass estimates obtained from the error-convolved runs) when the mere $`32`$ data point strong 3-D velocity data is used.

A steeper distribution function implies lower value of $`f`$ at most energies than does a flatter distribution function. Hence the flatter DF will imply higher likelihoods. Thus the likelihoods collected during the final “equilibrium” phase of Metropolis, with the steeper initial models, fall short of the simulated likelihoods. Since the limitation in the number of data points is the chief factor responsible for putting a brake on the mutual exclusivity of the results and the initial models, all the considered models pass the goodness-of-fit test with the bigger radial velocity sample during the error-convolved runs than the smaller proper motion sample in Ghez et al. (1998). However, when the number of data points goes up, the “true” fit to the observed data, as sought by Metropolis, deteriorates since the maxima of $`L`$ gets progressively narrower with increasing number of kinematical data points, (all of the considered models pass the goodness-of-fit test when the much smaller 3-D motion data set is used). Thus there is a trade-off between the goodness of fit achieved by Metropolis and the independence of the results from the initial models.

It is also observed that a greater number of the initial models appear to confirm consistency of the data with isotropy (via the goodness-of-fit test), when the projected distribution functions are error convolved. As explained above, the introduction of noise in the data in the form of errors in the velocity measurement causes an enhancement in the stellar energies. Thus for a given initial model, noise pushes down the value of the likelihood. When this noise is compensated for, (via convolution of the error distribution and the projected DF), the likelihood function is expected to rise to its correct value again, which may fall within the range of the simulated likelihood sample.

## 8 Summary

We propose a new method for analyzing discrete measurements of stellar radial velocities and/or proper motions and apply it to Galactic center stellar dynamics, where there are currently $`200`$ radial velocities, $`100`$ proper motions, and $`30`$ three-dimensional velocities now available. Both the distribution function and the underlying potential are recovered non-parametrically. This problem was previously intractable because of the technical difficulty of searching through a huge space of possible potentials and distribution functions; we overcome this difficulty by applying the Metropolis algorithm.

In this paper we consider the simplest case of a spherical isotropic system. We have designed a goodness-of-fit test to check whether the data exclude this case, based on which we conclude that this assumption is adequate under current data.

From proper motions reported by Ghez et al. (1998), we estimate $`2.0\pm 0.7\times 10^6\mathrm{M}_{}`$ within a radius of $`0.0044`$pc when observational errors are incorporated into the analysis. From radial velocities reported by Genzel et al. (2000), we estimate $`2.2_{1.0}^{+1.6}\times 10^6\mathrm{M}_{}`$ within $`0.046`$pc, and from three-dimensional velocities reported again in Genzel et al. (2000), the mass obtained is $`1.8_{0.3}^{+0.4}\times 10^6\mathrm{M}_{}`$ within $`0.046`$pc.