Non-parametric reconstruction of distribution functions from observed galactic discs

Key words: methods: data analysis – methods: numerical – galaxies: general – galaxies: kinematics and dynamics. 1 I N T RO D U C T I O N In years to come, accurate kinematical measurement of nearby disc galaxies will be achievable with high-resolution spectroscopy. Measurement of the observed line profiles will yield relevant data with which to probe the underlying gravitational nature of the interaction holding the galaxy together. Indeed the assumption that the system is stationary relies on the existence of invariants, which put severe constraints on the possible velocity distributions. This is formally expressed by the existence of an underlying distribution function which specifies the dynamics completely. The determination of realistic distribution functions which account for observed line profiles is therefore required in order to understand of the structure and dynamics of spiral galaxies. Inversion methods have been implemented for spheroids (globular clusters or elliptical galaxies) by Merrifield (1991), Dejonghe (1993), Merritt (1996, 1997), Merritt & Tremblay (1993, 1994), Emsellem, Monnet & Bacon (1994), Dehnen (1995), Kuijken (1995) and Qian (1995). Indeed, for spheroids, the surface density alone yields access to the even component of a two-integral distribution function which may account for the internal dynamics (while the odd component can be recovered from the mean azimuthal flow). However, the corresponding recovered distribution might not be consistent with higher Jeans moments, since the equilibria may involve three (possibly approximate) integrals. The inversion problem corresponding to a flattened spheroid which is assumed to have two or three (Stakel-based) integrals has been addressed recently by Dejonghe et al. (1996) and is illustrated by NGC 4697. Non-parametric approaches have in particular been used with success by Merritt & Gebhardt (1994) and Gebhardt et al. (1996) to solve the dynamical inverse problem for the density in spherical geometry. If the spheroid is seen exactly edge on, Merritt (1996) has devised a method which allows one to recover simultaneously the underlying potential. Mon. Not. R. Astron. Soc. 301, 419–434 (1998) q 1998 RAS *E-mail: [email protected] Here the inversion problem for thin and round discs is addressed for cases where symmetry ensures integrability. In this context, the inversion problem is truly two-dimensional and requires special attention for the treatment of quasi-radial orbits in the inner part of the galaxy. By Jeans’ theorem the steady-state mass-weighted distribution function describing a flat galaxy must be of the form f 1⁄4 f ð«; hÞ, where the specific energy, «, and the specific angular momentum, h, are given by « 1⁄4 2ðv 2 R þ v 2 fÞ 1 w ; h 1⁄4 R vf : ð1Þ Here vR and vf are the star radial and angular velocities respectively of stars confined to a plane and wðRÞ is the gravitational potential of the disc. The azimuthal velocity distribution, FfðR; vfÞ, follows from this distribution according to FfðR; vfÞ 1⁄4 … f ð«; hÞ dvR ; ð2Þ where the integral is over the region 1=2 × ðvR þ vfÞ < w corresponding to bound orbits. Pichon & Lynden-Bell (1996) demonstrated that, in the case of a thin round galactic disc, the distribution can be analytically inverted to yield a unique f ð«; hÞ provided the potential wðRÞ is known. The velocity distribution FfðR; vfÞ can be estimated – within a multiplicative constant – from line-of-sight velocity distribution (LOSVD) data obtained by long-slit spectroscopy when the slit is aligned with the major axis of the galactic disc projected on to the sky. Similarly, the rotation curve observed in H i gives in principle access to the underlying potential. More generally, simultaneous measurements of velocity distributions are derived with slits presenting arbitrary orientations with respect to the major axis, as discussed in Appendix C. The inversion of equation (2) is known to be ill-conditioned: a small departure in the measured data (e.g. caused by noise) may produce very different solutions since these are dominated by artefacts corresponding to the amplification of noise. Some kind of balance must therefore be found between the constraints imposed on the solution, in order to deal with these artefacts on the one hand and the degree of fluctuations consistent with the assumed information content of the signal on the other hand (i.e. the worse the data quality, the lower the informative content of the solution and the greater the constraint on the restored distribution so as to avoid an over-interpretation of the data). Finding such a balance is called the ‘regularization’ of the inversion problem (e.g. Wahba & Wendelberger 1979) and methods implementing adaptive level of regularization are described as ‘non-parametric’. Under the assumption that these discs are axisymmetric and thin, the proposed non-parametric methods described in this paper yield in principle a unique distribution: the smoothest solution consistent with all the available observables, the knowledge of the level of noise in each measurement and some objective physical constraints that a satisfactory distribution should fulfil. Section 2 presents all relevant theoretical aspects of regularization and non-parametric inversion for galactic discs distributions. Section 3 present the various algorithms and the corresponding numerical techniques, which we implemented in steps to carry efficiently this two-dimensional minimization. It corresponds in essence to an extension of the work of Skilling & Bryan (1984) for maximum entropy to other penalizing functions that are more relevant in this context. All techniques are implemented in Section 4 on simulated data arising when the slit of the spectrograph is aligned with the long axis of the projected disc. A discussion follows. 2 N O N PA R A M E T R I C I N V E R S I O N F O R F L AT A N D RO U N D D I S C S The non-parametric inversion problem involves finding the best solution to equation (2) for the distribution function when only discretized and noisy measurements of FfðR; vfÞ are available. A distinction between parametric and non-parametric descriptions may seem artificial: it is only a function of how many parameters are needed to describe the model with respect to the number of independent measurements. In a parametric model there is a small number of parameters compared with the number of data samples. This makes the inversion for the parametric model somewhat regularized, i.e. well-conditioned. Once the model has been chosen, however, there is no way to control the level of regularization and the inversion will always produce a solution, whether the parametric model and its implicit level of regularization is correct or not. In a non-parametric model, as a result of the discretization, there is also a finite number of parameters but it is comparable to and usually larger than the number of data samples. In this case, the amount of information extracted from the data is controlled explicitely by the regularization. Here the latter non-parametric method is therefore preferred, because no particular unknown physical model for disc distributions is to be favoured. 2.1 The discretized kinematic integral equation Since « is an even function of vR and since the relation between vR and « is one-to-one on the interval vR [ 1⁄20;∞Þ and for given R and vf, equation (2) can be rewritten explicitly as FfðR; vfÞ 1⁄4  2 p …0 1YðR;vfÞ f ð«;R vfÞ  « þ YðR; vfÞ p d« ; ð3Þ where the effective potential is given by YðR; vfÞ 1⁄4 wðRÞ 1 1 2 vf : ð4Þ For a given angular momentum h the minimum specific energy is «minðhÞ 1⁄4 min R[1⁄20;∞Þ n h 2R2 1 wðRÞ o