Spectroastrometry of rotating gas disks for the detection of supermassive black holes in galactic nuclei
I. Method and simulations
A. Gnerucci^{1}  A. Marconi^{1}  A. Capetti^{2}  D. J. Axon^{3,4}  A. Robinson^{3}
1  Dipartimento di Fisica e Astronomia, Università degli Studi di
Firenze, Firenze, Italy
2  INAF  Osservatorio Astronomico di Torino, Strada Osservatorio 20,
10025 Pino Torinese, Italy
3  Physics Department, Rochester Institute of Technology, 85 Lomb
Memorial Drive, Rochester, NY 14623, USA
4  School of Mathematical & Physical Sciences, University of
Sussex, Falmer, Brighton, BN2 9BH, UK
Received 19 May 2009 / Accepted 11 December 2009
Abstract
This is the first in a series of papers in which we study the
application of spectroastrometry in the context of gas kinematical
studies aimed at measuring the mass of supermassive black holes. The
spectroastrometrical method consists in measuring the photocenter of
light emission in different wavelength or velocity channels. In
particular we explore the potential of spectroastrometry of gas
emission lines in galaxy nuclei to constrain the kinematics of rotating
gas disks and to measure the mass of putative supermassive black holes.
By means of detailed simulations and test cases, we show that the
fundamental advantage of spectroastrometry is that it can provide
information on the gravitational potential of a galaxy on scales
significantly smaller (1/10)
than the limit imposed by the spatial resolution of the observations.
We then describe a simple method to infer detailed kinematical
informations from spectroastrometry in longslit spectra and to measure
the mass of nuclear mass concentrations. Such method can be applied
straightforwardly to integral field spectra, which do not have the
complexities due to a partial spatial covering of the source in the
case of longslit spectra.
Key words: line: profiles  techniques: high angular resolution  techniques: spectroscopic  galaxies: active  galaxies: nuclei  galaxies: kinematics and dynamics
1 Introduction
One of the fundamental open questions of modern astrophysics is understanding the physical processes that transformed the nearly homogeneous primordial medium into the presentday universe, characterized by a wealth of complex structures such as galaxies and clusters of galaxies. Understanding how galaxies formed and how they become the complex systems we observe today is therefore a major theoretical and observational effort. In the last few years, strong evidence has emerged for the existence of tight links between supermassive black holes (BH), nuclear activity and galaxy evolution. These links reveal what is now called as coevolution of black holes and their host galaxies. Strong evidence is provided by the discovery of 'relic' BHs in the center of most nearby galaxies, and that BH masses () are tightly proportional to structural parameters of the host spheroid like mass, luminosity and stellar velocity dispersion (e.g. Kormendy & Richstone 1995; Gebhardt et al. 2000; Ferrarese & Merritt 2000; Marconi & Hunt 2003; Ferrarese & Ford 2005; Graham & Driver 2008, and references therein). Moreover, while it is widely accepted that Active Galactic Nuclei (AGN) are powered by accretion of matter on a supermassive BH, it has been possible to show that BH growth is mostly due to accretion of matter during AGN activity, and therefore that most galaxies went through a phase of strong nuclear activity (Soltan 1982; Yu et al. 2004). It is believed that the physical mechanism responsible for this coevolution of BHs and their galaxies is probably the feedback by the AGN, i.e. the accreting BH, on the host galaxy (Silk & Rees 1998; Fabian 1999; Granato et al. 2004; Di Matteo et al. 2005; Menci et al. 2006; Bower et al. 2006).
In order to proceed further it is important to secure the most evident sign of coevolution, the correlations between BH mass and galaxy properties which can be achieved by increasing the number, accuracy and mass range of existing measurements. Supermassive BHs are detected and their masses measured by studying the kinematics of gas or stars in galaxy nuclei and, currently, about 50 BH mass measurements most of which in the range and only very few measurements below and above those limits (see, e.g., the most recent compilation by Graham 2008).
This paper, the first in a series, deals with BH mass measurements from gas kinematics and presents a new method, based on spectroastrometry, which can provide a simple but accurate way to estimate BH mass and which partly overcomes the limitations due to spatial resolution which plague the ``classical'' gas (or stellar) kinematical methods.
In Sect. 2 we introduce the standard method for gas kinematical studies, that is based on the gas rotation curves, and briefly discuss its characteristics and limitations. In Sect. 3 we introduce the new gas kinematical method based on spectroastrometry. We explain the basis of this approach that consist on measuring ``spectroastrometric curves'' (Sect. 3.1) and show simulations based on a model of the gas dynamics from which we want to learn how the spectroastrometric curve changes in function of some parameters of the model (Sect. 3.2). In Sect. 4 we explain the practical application of the method. We start by presenting a method for using simultaneously several spectroastrometric curves of the same source (Sect. 4.1) and we consider its application on noisy data (Sect. 4.2). In Sect. 5 we present a trivial fitting method for recovering the values of the various model parameters using long slit spectroscopy with noisy data. Finally, in Sect. 6 we describe the practical application of the method with Integral Field Units (IFU's) and in Sect. 8 we draw our conclusions. In Appendix A we discuss in details of the spectroastrometric measurements with special regard to the determination of the light centroids and in Appendix B we describe in detail the method to recover 2D spectroastrometric maps from multiple slit spectra.
2 The standard gas kinematical method
Figure 1: Example of the measurements on which the construction of rotation and spectroastrometric curves is based. Upper left panel: the longslit spectrum of a given continuumsubtracted emission line (isophotes) with superimposed two particular sections at fixed slit position and velocity (dotted lines). Bottom left panel: the profile relative to the section at fixed slit position with superimposed the position of its centroid (dashed line). Upper right panel: the profile relative to the Sect. at fixed velocity with superimposed the position of its centroid (dashed line). Bottom right panel: the positionvelocity diagram (PVD) of the upper left panel with superimposed rotation (blue dots with error bars) and spectroastrometric curves (red points with error bars). In all panels velocities are measured relative to the systemic velocity of the galaxy which is set to the zero point of the velocity scale. 

Open with DEXTER 
The standard method of gas kinematics is based on spatially resolved spectra of the nuclear region of a galaxy. This method consists in recovering the rotation curve of a given gas emission line from a longslit spectrum. The modeling assumes that the gas is rotating in a thin disc configuration (neglecting hydrodynamical effects) under the effect of the gravitational potential of the stellar mass and of a pointlike dark mass that is the BH. The value of and other unknown parameters of the model are obtained by fitting the rotation curves (see Marconi et al. (2006) and references therein for details). Because of the many unknown parameters of the model, to better constrain the fit, usually many spectra of the galactic nucleus are obtained each with a different orientation of the spectroscopic slit and a simultaneous fit of all of them is performed.
Clearly, the ability to detect the presence of a BH and measuring its mass strongly depends on the signaltonoise ratio of the data and only poor constraints on can be obtained from low S/N data. However, even with high S/N data, the fundamental limit of the standard gas kinematical method (``rotation curves method'' hereafter) resides in the ability to spatially resolve the region where the gravitational potential of the BH dominates with respect to the contribution of the stars. Similar considerations apply equally to stellar dynamical measurements.
As a rule of thumb, this region corresponds to the socalled
sphere of influence of the BH.
The radius of the sphere of influence
(
)
can be estimated as (Binney J.
1987 ):
where G is the gravitational constant, the BH mass and is the velocity dispersion of the stars in the galaxy. For a galaxy with BH mass and stellar velocity dispersion km s^{1} we obtain pc which, for a very nearby galaxy at distance Mpc, provides an apparent size of . If a typical very good spatial resolution available from ground based observations is of the order of (this Full Width Half Maximum (hereafter FWHM) value should then be compared with twice ), we can notice that the sphere of influence is only marginally resolved even for very nearby galaxies with moderately large BHs. For a galaxy distance of 30 Mpc the sphere of influence becomes marginally resolved even with the Hubble Space Telescope (HST) which provides the best spatial resolution currently available from space. Therefore, the rotation curves method can detect only BHs with moderately high masses and located in nearby galaxies; this is a strong limitation from the point of view of a ``demographic'' study of BHs in galactic nuclei, because the adopted ``investigation tool'' cannot reveal the entire BH population.
3 The spectroastrometric gas kinematical method
The use of spectroastrometry was originally introduced by Beckers (1982), Christy et al. (1983) and Aime et al. (1988) to detect unresolved binaries. However these earlier studies required specialist instrumentation, and it was not until the work of Bailey (1998) that this method was exploited using standard common user instrumentation: a longslit CCD spectrograph.
Subsequently spectroastrometry has been used by several authors to study pre main sequence binaries (Porter et al. 2004,2005; Baines et al. 2004) and the presence of inflow or outflow or the disk structure on the gas surrounding pre main sequence stars (Whelan et al. 2005; Takami et al. 2003). More recently, Brannigan et al. (2006) discussed the presence and the detection of artifact in the output of the method and Pontoppidan et al. (2008) used spectroastrometry on m CO lines to study the kinematical properties of protoplanetary disks.
The fundamental advantage of the spectroastrometric method is that, in principle, it can provide position measurements on scales smaller than the spatial resolution of the observations.
We will now explain the general principle of the spectroastrometric method with a simple example: consider two pointlike sources located at a distance smaller than the spatial resolution of the telescope; these sources will be seen as spatially unresolved with their relative distance not measurable from a conventional image. However, if in the two sources are present spectral features, such as absorption or emission lines at different wavelengths, the light profiles extracted from a longslit spectrum at these wavelengths will show the two sources separately. From the difference in the centroid of the light profiles at these two wavelengths one can estimate the separation between the two sources even if this is much smaller that the spatial resolution. This ``overcoming'' of the spatial resolution limit is made possible by the ``spectral'' separation of the two sources.
While this method has been applied to unresolved pointlike sources as binary stars and protostellar systems, it has never been applied to the problem of measuring BH masses from the nuclear gas emission in galaxies. Motivated by the need to overcome the shortcomings of the standard gas kinematic method described above in Sect. 2, in this paper we study the application of the spectroastrometric method to gas kinematical studies of BH's in galactic nuclei. We will start presenting the spectroastrometric method based on the same longslit spectra used for the rotation curves method but instead of the rotation curves we recover the ``spectroastrometric curve''. However, as shown in Sect. 6, the application of spectroastrometry to integral field data is even more simple and straightforward than to longslit spectra.
Here we focus on the theory of the method and the development of a practical framework for its application, exploring its capabilities and limitations using simulated data in order to understand how the spectroastrometric curve is affected by the object itself or by the instrumental setup. In subsequent papers we will apply the method to the kinematical data of real galaxies to yield new improved BH mass measurements.
In this paper we concentrate on the application of the spectroastrometric method to continuum subtracted spectra in order to focus exclusively on the gas kinematics. Indeed, among other things, the underlying continuum only dilutes or modifies the spectroastrometric signal expected for the spatial distribution and kinematics of emission line gas.
3.1 The spectroastrometric curve
Here we explain how the spectroastrometric curve is obtained and its main differences with ``classical'' rotation curves.
A longslit spectrum of a continuum subtracted emission line provides a pixel array whose axes map the dispersion and slit directions, the socalled positionvelocity diagram (hereafter PVD). The slit axis maps the observed position along the slit. The dispersion axis maps the wavelength of emission from which the line of sight velocity can be derived. The upper left panel in Fig. 1 shows the isophotes of an emission line from a PVD. Ideally, for infinite S/N and perfectly circularly rotating gas, the rotation curve denotes the mean gas velocity as a function of the position along the slit. In practice, that curve is obtained by fitting the observed line profiles along the slit with Gaussian functions which provide an estimate of the average velocity after discarding components which are clearly not circularly rotating. The lower left panel of Fig. 1 displays the line profile extracted at the slit position marked by the horizontal dotted line in the PVD. The spectroastrometric curve provides the mean position of the emitting gas as a function of velocity. In ideal data that curve is derived by taking the light profile of the line at given velocities along the dispersion direction and measuring the corresponding mean emission centroids. In practice, the finite S/N of real data requires more complex measurements, which are described below. The upper right panel of Fig. 1 displays the light profile extracted at the velocity marked by the vertical dotted line in the PVD. The two curves can be compared in lower right panel superimposed on the isophotes in the PVD diagram. Clearly, the two methods analyze the same spectrum from complementary points of view.
3.2 Simulations
In the following we show the results of the tests based on simulations of rotating gas disks, outlining the effect of varying the free parameters of the simulation and of the spectroscopic observations.
The model we use in our simulations is based on a thin gas disk rotating in a plane under the effect of the gravitational potential produced by a mass distribution of stars in a galaxy nucleus and by a pointlike dark mass, the central BH . Therefore, we assume that the gas is circularly rotating in the disk plane with the rotational velocity uniquely determined by the combination of and the stellar mass distribution. This model depends on several parameters: the dynamical parameters (the BH mass, the shape of the mass density function of the stars and the systemic velocity of the galaxy) and the geometrical parameters that establish the position and orientation of the gas disk (the distance of the galaxy or the angular distance scale, the inclination of the disk plane with respect to the line of sight and the orientation of the line of nodes of the disk). The model takes into account the effect of the shape of the intrinsic light distribution of the emission line on the sky plane, that is modeled analytically with a combination of Gaussian or exponential functions. The model takes also into account the effect of the instrumental Point Spread Function that is modeled with a Gaussian function with a given FWHM, and the effect of the others instrumental setup parameters: the spectral resolution, the slit width and position angle, the detector's pixel size. With this model we can simulate a longslit spectrum of a gas emission line in a particular galaxy nucleus (see Marconi et al. 2006, and references therein for a detailed description).
Figure 2: Spectroastrometric curve (black solid line) obtained from a simulated spectrum for a rotating gas disk with various values of the central BH and stellar mass. In gray we plot the isophotes of the simulated line spectrum. Upper left panel: case of BH without mass contribution from stars. Upper right panel: stars only, with . Lower panels: ( left) and ( right) with included the same stellar mass distribution as in model in the upper right panel. Note that the differences in the kinematics observed in the left panels (both BH dominated, but with and without stellar contributions) are marginal, and only visible i.e. beyond from the nucleus. 

Open with DEXTER 
Figure 3: Left panel: spectroastrometric curves from spectra differing only in spatial resolution. Solid line: spatial resolution of . Dotted line: spatial resolution of . Dashed line: spatial resolution of . The vertical longdashed lines denote the ``high velocities''. Right panel: corresponding rotation curves. 

Open with DEXTER 
For simplicity, in order to illustrate the technique, we make use of a basic reference model which is meant to closely match the physical parameters of a BH of mass in an elliptical galaxy in the local universe observed with a typical longslit spectrograph, like ISAAC at the VLT (Moorwood et al. 1999). The key model parameters are:
 the distance of the galaxy is set to 3.5 Mpc that corresponds to an angular distance scale of 17 pc/ (e.g. the same of the galaxy CentaurusA);
 the disk inclination is set to ;
 the disk line of nodes position angle (with respect to the North direction) is set to ;
 the amplitude of rotational velocity in the disk due to the stellar mass component at is set to 200 km s^{1};
 the FWHM of the spatial PSF is set to as for typical high quality ground based observations;
 the spectral resolution is set to 10 km s^{1};
 the detector's angular pixel size is set to , resulting in a spatial oversampling of 3.5;
 the spectrograph slit is usually set parallel to the disk line of nodes unless otherwise specified.
Figure 4: Left panel: spectroastrometric curves for spectra differing in the value BH mass. Solid line: . Dotted line: . Right panel: rotation curves for the same models. 

Open with DEXTER 
3.2.1 BH mass
We first consider the effect of changing the value of the central BH mass on the spectroastrometric curve.
We simulated the spectrum of a particular rotating gas disk model with different values of the mass of the central BH, and then we derived the spectroastrometric curves from these simulated spectra (see Fig. 2).
The spectroastrometric curve in Fig. 2 is centered at the systemic velocity of the galaxy. In these noiseless simulations, used for illustration purposes, the spectroastrometric curves are drawn only when the mean flux of the light profile along the slit is larger than 10^{3} of the maximum flux of the spectrum. This is meant to give a representation of the low signaltonoise regions that should be used in real data for the measurement of the centroid of the light profile, without considering unrealistically low flux levels. We defer a discussion of the practical cut off fluxes imposed by the signaltonoise in real data to Sect. 4.2.
In the case (Fig. 2 upper right panel) the gas kinematics is due only to the gravitational potential of the stellar mass distribution and the spectroastrometric curve is monotonically decreasing at increasing velocity. In all other cases (Fig. 2 other panels), the points at high velocities ( km s^{1} and km s^{1}) tend to approach the value for increasing velocities (we consider ``increasing'' referring to ; for this simulation), giving the curve the characteristic ``S''shape. The reason for this behavior is that when the gas kinematics is dominated by the BH's gravitational potential, the gas at high velocity is located close to the BH and the light centroid of the high velocity gas has to be near to the BH position.
Operatively the ``high velocity'' points are those where the line emission along the slit is spatially unresolved. In Appendix A we describe how to find these points as a by product of centroid determination. We can also recall that the presence of a turnover in the spectroastrometric curve is the signature of the presence of a BH. Therefore, the ``high velocity'' range is also characterized by a drop of the spectroastrometric signal with respect to the maximum shift reached at the turnover point.
The velocity in a circular orbit of radius r around a pointlike mass (the BH) is well approximated by a keplerian law with in the innermost region of the gas disk where we can neglect the contribution of the gravitational potential of the stellar mass distribution. Combined with the above considerations, this also suggests that the spectroastrometric curve asymptotic value at high velocities provides an estimate of the BH position along the slit.
For increasing values the spectroastrometric curve extends to higher velocities (higher values of ) since this has the effect of increasing the amount of emission in the high velocity bins. However the limited extension in the velocity axis of the spectroastrometric curve with real data is due to the presence of noise. At velocities where the flux of the line is too low with respect to noise we will not be able to calculate a reliable value of the centroid of the light profile along the slit.
In conclusion the spectroastrometric curve reveals the presence of a pointlike mass contribution to the gravitational potential when it shows an ``Sshaped'' structure, with a turnover of the high velocities components that get closer to the . It can then be concluded that the information about the BH resides predominantly in the ``high velocity'' part of the curve.
3.2.2 Spatial resolution
Here we consider the effects of the spatial resolution on the spectroastrometric curve. The spatial resolution is the width (FWHM) of the point spread function which, in the model, is approximated by a Gaussian function.
The results of this test are shown in Fig. 3 where we compare spectroastrometric and normal rotation curves. In the left panel we present the spectroastrometric curves for models differing only in spatial resolution while in the right panel we display the corresponding rotation curves. The spectroastrometric curves differ in the ``low velocity'' range ( km s^{1}), as they show a steeper gradient and a large amplitude.
However, in the ``high velocities'' range differences are negligible, 0.1 pix at most. In Sect. 3.2.1 we showed how the information on the BH mass is encoded in the high velocity range and, in particular, the presence of the BH is revealed by the fact that the spectroastrometric curve approaches at high velocities. The spectroastrometric curve at the ``high velocities'' is almost unchanged by worsening the spectral resolution, leaving the BH signature unaltered.
For the ``standard'' rotation curves the information about the BH is also encoded in the presence of points at high velocities at low distances from the center. However, a lower spatial resolution the BH signature is effectively canceled in the rotation curves when the sphere of influence of the BH is not resolved. In this simulation the apparent radius of the sphere of influence is , which is resolved in the case of spatial resolution, partially resolved for case, and unresolved for .
Clearly, even in the case of the spectroastrometric curves a better spatial resolution is desirable. In fact, a poorer spatial resolution results in a broadening of the light profile along the slit and in a decreased accuracy with which the centroid position (and consequently the BH mass) can be measured. Note, however, that the spatial resolution is not a free parameter, since it can be measured directly from the data. Its effects can be modeled and taken into proper account.
We have just shown that the the rotation curve changes drastically when varying the spatial resolution while the spectroastrometric curve does not. Now we conclude this section with a last set of simulations to showing in a qualitative but more accurate way how this method can really allow us to overcome the spatial resolution limit.
The results of these simulations are shown in Fig. 4. In the left panel we display the spectroastrometric curves for two models with the same spatial resolution ( ) but different BH mass; in the right panel we display the corresponding rotation curves.
By considering only the rotation curves, the case is effectively indistinguishable from that with and indeed the sphere of influence is unresolved in this case (the apparent radius of the sphere of influence is ).
Instead, for this particular set of simulations, a mass appears to be still distinguishable from the case; in the high velocity range the differences between the two models reach a value of , i.e. 1.3 pixels. Such difference can be measured in real data, assuming that we can achieve a reasonable accuracy of 1 pixel in the measure of the photocenter at these velocities. In conclusion, according to these simulations we are able to detect a BH whose apparent size of the sphere of influence is less than 1/10 of the spatial resolution.
3.2.3 Slit position angle
Here we consider how the spectroastrometric curve is affected by changing the position angle of the slit. We simulated the spectrum of a particular rotating gas disk model with different values of the position angle of the slit and then we get the spectroastrometric curves from this simulated spectra.
In the left panel of Fig. 5 we can see the comparison of the spectroastrometric curves relative to different values of the slit position angle referred to the disk line of nodes position angle ( where and are the position angles of the slit and disk line of nodes respectively).
We can see that for (slit aligned with the disk line of nodes) the spectroastrometric curve has the maximum amplitude. For increasing values of the amplitude decreases reaching a null value for . For values larger than the curve invert itself. This is clearly a geometrical projection effect just like the one affecting the amplitude of normal rotation curves as a function of . In the right panel of Fig. 5 we show the curve amplitude as a function of which is well approximated by a cosine function (dashed line).
Figure 5: Left panel: simulated spectroastrometric curves from spectra differing only in slit position angle (referred to the disk line of nodes ). Right panel: measured amplitude of the spectroastrometric curve (distance between the maximum and the minimum) as a function of the slit position angle referred to the disk line of nodes; the dashed line denotes a cosine function opportunely rescaled to match the value of the amplitude of . 

Open with DEXTER 
Figure 6: Left panel: isovelocity contour map of the line of sight velocity field on the plane of sky for a rotating gas disk with the of line of nodes along the x axis, with superposed a wide slit. Right panel: comparison of spectroastrometric curves for simulated spectra differing only in slit width (we remind that we assumed a spatial resolution of 0 ). Solid line: wide slit. Dotted line: wide slit. Dashed line: wide slit. 

Open with DEXTER 
3.2.4 Slit width
Another parameter which influences the spectroastrometric curve is the width of the slit with which the spectra are obtained. Each point of the spectroastrometric curve represents the centroid of the light profile along the slit at a given velocity but one can only sample the fraction of light emitted by the gas at that specific velocity that is intercepted by the slit.In Fig. 6, left panel, we show the map of line of sight velocity field of a rotating gas disk. The gas in a given velocity bin lies in the locus delimited by two subsequent isovelocity contours. The impact of superimposing the slit on the isovelocity contour map is to artificially truncate the spatial regions contributing to line emission at a given velocity that are not confined within the slit extension. This results in distorted and displaced photocenters in velocity space.
In Fig. 6, right panel, we show the spectroastrometric curves obtained from models differing only in slit width. As long as the slit width is smaller or equal to the spatial resolution of the observations (which is, in this case, FWHM) then differences among spectroastrometric curves in the high velocity range are negligible (e.g., slits with , width). On the contrary the spectroastrometric curve is significantly affected for larger slit widths, even in the high velocity range, because of the inclusion of more extended emission. This comparison indicates that one should select for the observations a slit width smaller or at most equal to the spatial resolution of the observations. As already noted discussing the effects of spatial resolution, the small residual differences in the spectroscopic curves corresponding to various slit widths can be effectively modeled out.
3.2.5 Spectral resolution
Figure 7: Comparison of spectroastrometric curves for spectra differing only in spectral resolution. Solid line: spectral resolution of 10 km s^{1}. Dotted line: spectral resolution of 50 km s^{1}. Dashed line: spectral resolution of 100 km s^{1}. Dotdashed line: spectral resolution of 150 km s^{1}. 

Open with DEXTER 
We now focus on the effects of the finite spectral resolution of the
observations on the spectroastrometric curves. In Fig. 7 we show
models differing only in spectral resolution. By decreasing the
spectral
resolution, the amplitude of the spectroastrometric curves is decreased
and
they show, at a given spatial offset, higher velocities. Effectively,
due to
the spectral convolution, the curves are stretched along the velocity
scale.
To understand this effect we can approximate the true velocity profile
and the
instrumental line profile with Gaussian functions. The resulting line
profile
is the convolution of these two functions and is therefore a Gaussian
function
with standard deviation
where is the resulting standard deviation, is the standard deviation of the velocity profile and is the standard deviation of the instrumental response function.
We here consider as a source of the line broadening the unresolved rotation of the gas that originates from the fact that one is observing with finite spatial resolution and is then not able to spatially resolve the high velocity regions close to the BH. In our simulations, the line width due to unresolved rotation reaches km s^{1} at the galaxy's center.
As expected, the spectroastrometric curves are almost unchanged at high spectral resolution ( km s^{1}) but they are significantly altered for the lower spectral resolutions considered ( km s^{1}) when approaches the value of . As a consequence, the BH mass estimates derived from data of insufficient spectral resolution are systematically overestimated, since at a given position, one measures an artificially increased velocity.
This result underscores the importance of using data of as high as possible spectral resolution. The optimal value must be however derived tradingoff with the level of signaltonoise necessary to build well defined spectroastreometric curves.
Finally, we note that the value of can be estimated by modeling the classical rotation curves. It is then possible to establish, a posteriori, whether the BH measurement is affected by such an effect and eventually to validate its value.
Figure 8: Comparison of spectroastrometric curves for spectra differing only in intrinsic velocity dispersion. Solid line: intrinsic velocity dispersion of 10 km s^{1}. Dotted line: intrinsic velocity dispersion of 50 km s^{1}. Dashed line: s intrinsic velocity dispersion of 100 km s^{1}. Dotdashed line: intrinsic velocity dispersion of 150 km s^{1}. 

Open with DEXTER 
3.2.6 Intrinsic velocity dispersion
Another source of line broadening might be due to non circular or chaotic motions in the gas which can be modeled by adding an ``intrinsic'' velocity dispersion to the gas motions (e.g. see the discussion in Marconi et al. 2006). In Fig. 8 we show the spectroastrometric curves for models differing only in intrinsic velocity dispersion. The comparison with Fig. 7 clearly shows that, when adopting the same values of , no significant difference is found regardless of whether the source of line broadening is poor spectral resolution or an intrinsic velocity dispersion in the source.
In general, when spurious line broadening is present, the high velocity points in the spectroastrometric curve are not entirely due to the gravitational potential of the BH, and the curve is artificially stretched in the velocity direction leading to a possible overestimate of the BH mass with the method described in Sect. 5.
From the simulations presented in Figs. 7 and 8 we can verify that
the
effect of line broadening on spectroastrometric curves is indeed
approximable as a xaxis ``stretching''. Following
Eq. (2)
we verified that can then recover the ``destretched'' velocities as:
thus correcting the spectroastrometric curve. A detailed analysis of the effects of spurious line broadening is beyond the scope of this paper and will be presented elsewhere but, briefly, we can in principle estimate the spurious line broadening and that due to unresolved rotation by modeling the classical rotation curves (e.g. Marconi et al. 2006 and references therein) and then correct the spectroastrometric curves by destretching the velocity axis as indicated in Eq. (3). Clearly, the most accurate approach will be that of a combined fit of the classical rotation and spectroastrometric curve, much beyond the simple analysis presented here in Sect. 5.
3.2.7 Flux distribution
In this section we present simulations based on simple intrinsic flux distribution models, with the aim of showing how subresolution variations of the flux distribution influence the spectroastrometric curve.
The basic model adopted for the flux distribution in these simulations is an exponential function where r is the radial distance from the symmetry center located at the position and r_{0} is a characteristic radius which we assume equal to . At the spatial resolution of our simulations ( ), this emission is spatially unresolved.
Figure 9: Comparison of spectroastrometric curves for spectra differing only in intrinsic flux distribution. Solid line: base model flux distribution (explained in the text). Dotted line: model with . Dashed line: model with . Dotdashed line: base model with a central hole of radius . 

Open with DEXTER 
Figure 10: Comparison of spectroastrometric curves for spectra differing only in intrinsic flux distribution. Solid line: base model flux distribution (explained in the text). Dotted line: base model centered on position. Dashed line: base model centered on position. Dotdashed line: base model centered on position. 

Open with DEXTER 
In the first set of simulations, shown in Fig. 9, we vary the characteristic radius r_{0} in size and introduce a central hole in the flux distribution with radius . In the second set of simulations, shown in Fig. 10, we vary the position of the center of the flux distribution.
In the ``high velocity'' range the variation in the centroid position with respect to the base model is lower than and for the first and second set of simulations, respectively; these values should be compared with the pixel size, , and the spatial resolution, . As we will show in Sect. 5, these variations will produce only small changes (up to 0.2 dex) on the BH mass which can be inferred from spectroastrometric data. On the contrary, much larger variations are present in the ``low velocity'' range.
This behavior, observed also in previous simulations, is explained by the fact that that the gas with lower velocities is located in an extended (spatially resolved) region, and the light profiles along the slit are influenced by the shape of the flux distribution. Conversely the gas at higher velocities is confined in a small (spatially unresolved) region close to the BH and the light profiles along the slit are not influenced by the shape of the flux distribution.
Only in the case of a central hole in the flux distribution does the largest variation in the spectroastrometric curve occur in the high velocity range. This behavior is explained by the fact that that the gas with higher velocities (located in a small region close to the BH) is not illuminated because of the central hole of the flux distribution and the centroid position is influenced by the flux distribution at larger distance from the BH.
3.2.8 Summary of simulation results
We can now summarize the results from the above simulations and the considerations that can be derived.
 The presence of a black hole is revealed by a turnover in the spectroastrometric curve, with the high velocity components approaching a null spatial offset from the center of the galaxy. Conversely, in the case of no black hole, the curve shows a monotic behavior;
 the information about the BH is predominantly encoded in the ``high velocity'' range of the spectroastrometric curve which is generated by spatially unresolved emission;
 the ``high velocity'' range of the curve is not influenced strongly by the spatial resolution of the observations, leaving the BH signature unaltered. According to our simulations we are able to detect a BH whose apparent size of the sphere of influence is as small as 1/10 of the spatial resolution. Instead, for the ``standard'' rotation curves method the information about the BH is effectively canceled when the sphere of influence of the BH is not spatially resolved;
 the amplitude of the spectroastrometric curve decreases by increasing the angle between the slit and the line of nodes of the gas disk, according to a cosine law;
 the slit width affects the spectroastrometric curve, due to the truncation of the isovelocity contours. This effect can be reduced to a negligible level selecting for the observations a slit width smaller than (or equal to) the spatial resolution;
 the effect of the finite spectral resolution is to strech (in velocity space) the spectroastrometric curve. A robust BH mass estimates requires a velocity resolution smaller than the line width due to (spatially) unresolved rotation;
 variations of the intrinsic line flux distribution at subresolution scales, do not greatly affect the spectroastrometric curve in the high velocity range.
4 Practical application of the method
4.1 The use of multiple spectra of the source: the spectroastrometric map
In order to improve the constraints on model parameters with the ``standard method'', it is common practice to obtain several longslit spectra at different slit orientations. Also spectroastrometry benefits from this approach as it is possible to recover for the same source several spectroastrometric curves, one for each position angle (PA) of the slit. With respect to the standard rotation curves method, where spectra taken along paraller slits have been often used, here the fundamental requirement is that the various slits must be centered on the expected BH position.
We simulated the case of three longslit spectra, oriented parallel (PA1), perpendicular (PA2) and forming an angle of with the line of nodes (PA3). In Fig. 11 we can see the spectroastrometric curves obtained from the three spectra.
Figure 11: Spectroastrometric curves for three simulated spectra of the same model with different PA of the slit. Solid line: slit with PA (referred to the disk line of nodes position angle) . Dotted line: slit with . Dashed line: slit with . 

Open with DEXTER 
Figure 12: Spectroastrometric 2D map derived from the three spectroastrometric curves of Fig. 11. Left panel: derived photocenter positions on the sky plane, the black open diamonds are the points actually used for the minimization. Right panel: the 3D plot of the map, where the z axis is velocity. 

Open with DEXTER 
Each spectroastrometric curve provides the photocenter position along one slit, i.e. the position of the photocenter projected along the axis identified by the slit. Combing the spectroastrometric curves we can thus obtain the map of photocenter positions on the plane of the sky for each velocity bin. In principle, the spectroastrometric curves from two nonparallel slits should suffice but we can use the redundant information from the three slits to recover the 2D sky map as described in more detail in Appendix B.
In Fig. 12 (left panel) we present the 2D spectroastrometric map on the plane of the sky. All the points of the map lies on a straight line which identifies the line of nodes. Indeed, whenever the intrinsic distribution of line emission is circularly symmetric and centered on the BH position (as assumed in our simulations), the light centroids at a given velocity are located on the disk line of nodes. In the right panel of Fig. 12 we show the threedimensional representation of the spectroastrometric curve combining the 2D spatial map of the left panel with the velocity axis. In the 3D representation all the points lie on a plane parallel to the velocity axis and aligned with the disk line of nodes. On that plane, the spectroastrometric curve present the usual S shape symmetric with respect to the zero point of velocities.
We remark that in this paper we only apply the spectroastrometric method to continuum subtracted spectra. When dealing with real data an accurate continuum subtraction will be an important requirement to obtain an accurate spectroastrometric curve, especially for the ``high velocity'' points where the linecontinuum contrast is lower. Any residual continuum emission can alter the spatial profile of the emission line, and the ``turn over'' signature of a point mass in the spectroastrometric curve. In the latter case, the spectroastrometric displacement of the emission line will be shifted towards the continuum photocentre in the line wings no matter the actual spatial distribution of the high velocity gas. The accuracy of continuum subtraction, as well as the signal to noise ratio at high velocities, is therefore a fundamental limit to the ability of detecting the central BH. A more detailed analysis of the effects of continuum emission will be presented in a forthcoming paper where we will apply the spectroastrometric method to real data.
4.2 The effect of noise on spectroastrometry
Figure 13: Left row: spectroastrometric curves for the model of Fig. 11 with added noise for a S/N of 20, 50, and 100 ( from top to bottom). Solid line: slit at PA = . Dotted line: slit at PA . Dashed line: slit at PA = . For simplicity we plot the error bars of the points only for the PA curve. Right row: spectroastrometric 2D map derived from the three spectroastrometric curves of the left row. The red filled circle marks the derived position of the BH, while the red solid line represents the position angle of the line of nodes ( ) obtained by fitting a straight line to the XY position of the high velocity points. 

Open with DEXTER 
We now consider the effect of noise on spectroastrometric curves and maps, as in real data. We simulate noise by adding normally distributed random values to the simulated positionvelocity spectral maps. These random numbers are characterized by zero mean and standard deviation which is chosen in the following way. We first extract a ``nuclear'' spectrum by coadding line emission over an aperture equal to the PSF FWHM, centered on the nucleus position. The S/N ratio of the spectrum is then defined on the peak of the line profile, and the actual value per pixel is selected to provide a given S/N.
In Fig. 13 we show the spectroastrometric curves and maps for the same model presented in Fig. 11 to which we added an artificial noise in order to obtain S/N of 20, 50 and 100, respectively.
A key point to emphasize is that the curves (and maps) with higher S/N tend to have more measured points. This is because at the ``high velocities'' we have a low flux with respect to the line core and the presence of noise limits our ability to obtain a reliable light centroid position. In detail we stop estimating the light centroids when the mean flux of the light profile for that particular velocity bin is lower than the noise level. Thus, a spectrum with higher S/N enables us to obtain more accurate measurements on individual datapoints, but also extends the coverage at higher velocities, crucial to detect a BH and to measure its mass.
The 2D spectroastrometric maps can be used to estimate the geometrical parameters of the gas disk. In the following, we will always make use of simulations with noise.
As noted in Sect. 3.2.1, if the gas kinematics is dominated by rotation around a pointlike mass (the BH), the position of the light centroid at the highest velocities approaches the position of the BH. This consideration, which is also valid for the 2D spectroastrometric map, allows us to estimate the position on the plane of the sky of the BH as the average position of the ``high velocity'' points. In practice, we calculate this position by first taking the averages of the coordinates of the points in the blue and red ``high velocity'' ranges and then taking the average coordinates of the ``blue'' and ``red'' positions. The inferred BH position is marked by a red filled circle in the right row diagrams of Fig. 13.
Also, the high velocity points in Fig. 13 (and also Fig. 12) lie on a straight line which identifies the direction of the disk line of nodes (see Fig. 6 and Sect. 3.2.4). The red solid line shown in the right row of Fig. 13 represents the position angle of the line of nodes ( ) obtained by fitting a straight line to the high velocity points.
In the Table 1 we report the BH position and disk line of nodes position angle values obtained from the 2D map of the model of Fig. 13. We can observe that the accuracy of these measurements scales approximately with the square root of the S/N of the spectrum. Moreover in the maps with lower S/N the scatter of the points around the line of nodes increases (the points are almost perfectly aligned on the line of nodes direction in the noise free model of Fig. 10). This fact decreases the accuracy of the BH position and line of nodes PA estimates for lower S/N.
Table 1: Lines of Nodes and BH position estimates from noisy data.
5 Estimate of the BH mass from the spectroastrometric map
In this section we present a simple and straightforward method to recover the BH mass from the spectroastrometric measurements. The method is based on the following assumptions: (i) the gas is circularly rotating in a thin disk and (ii) each point of the spectroastrometric map represents a test particle on the disk line of nodes and is characterized by a line of sight (``channel'') velocity which is given by the center of the velocity bin where the light centroid was estimated. Under these assumptions it is trivial to relate the spectroastrometric map to the BH mass.
The circular velocity of a gas particle with distance r
from the BH is given
by:
where is the stellar mass enclosed in a sphere of radius r(assuming spherical symmetry) and can be written as:
with L(r) representing the radial luminosity density distribution in the galactic nucleus and M/L is the mass to light ratio of the stars.
The line of sight velocity (hereafter
for ``channel velocity'') of a
test particle located on the disk line
of nodes at distance r from the
BH is then given by
where i is the inclination of the disk plane and is the systemic velocity of the galaxy.
The centroid positions on the spectroastrometric map are
denoted by
and, in general, they are not perfectly aligned along the line of
nodes as described in Sects. 4.1
and B. We
select the ``high
velocity'' points and we then identify the line of nodes by a linear
fit of the
spectroastrometric
map. To reduce the correlation between
the values of intercept and slope of the fitted line we consider
where is the mean of the positions. and bdenote respectively the position angle and the intercept of the line of nodes to be determined by the fit. Then and b are recovered by minimizing the following where we take into account both errors in x and y:
where and are the uncertainties on the position of the points of the spectroastrometric map and, as observed in Sect. 4.1, the sum is extended only over the ``high velocities'' range.
The second step requires to associate each velocity bin to its
position
within the rotating disk. As described in Sect. 4.1, the 2D map points
lie with good
approximation on the line of nodes, particularly in the ``high
velocities'' range.
We then project the
position
of the 2Dmap points on the
line of nodes and calculate their coordinate ``S''
with respect to this
axis of reference, with ``S'' defined as:
Each has an uncertainty estimated by taking into account the uncertainties on the 2D map points position and as well as the uncertainties on the line of nodes parameters resulting from the fit. The distance r of a test particle from the BH used in Eq. (4) differ from only by a constant S_{0}, the unknown coordinate of the BH along the line of nodes, i.e.
Combining Eqs. (4)(6), and (10) we finally obtain the expected velocity in channel ch:
The unknown model parameters are:
and can be determined by minimizing the following
where is the uncertainty of the numerator computed as a function of the unknown parameters values (par) and the uncertainties on the positions along the line of nodes. We remark that the channel velocity has no associated unceratinty since it is not a measured quantity but is the central value of a velocity bin where the spectroastrometric curve is measured. As discussed in the previous sections, we restrict the fit (i.e. the sum over the velocity channels) to the ``high velocities'' range.
Figure 14: Fit results for models of Fig. 13 with S/N=20 ( top panel), S/N=50 ( middle) and S/N=100 ( bottom). Points with error bars are the result of ``observations''. The red line is the best fit model and the dotted line is the expected contribution from the stellar mass. 

Open with DEXTER 
Figure 15: Observed and model spectroastrometric curves at S/N=100 (e.g. vs. ). The horizontal dotted line is the BH ``S'' coordinate S_{0}. The vertical dotted line is the systemic velocity . 

Open with DEXTER 
The factor in Eq. (12) is much smaller for the points at lower velocities (i.e. smaller ) which are closer to the peak of the line profile and have therefore much larger S/N than the points at ``high velocity''. However, from the discussion in the previous sections, we know that the spectroastrometry is less reliable for the determination of BH properties like position and mass. Therefore, to avoid being biased by potentially faulty points we add in quadrature a constant error to . The value of is found by imposing that the reduced is equal to 1. Therefore, if the points at lower velocity with higher S/N are problematic they will not bias the final fit results, because their high weight will be greatly reduced by the addition od of a much larger .
After determining the best fit values of the free parameters,
the position of the BH on the plane of the sky is simply given
by
We have performed the fit of the spectroastrometric data for the models of Fig. 13, which are computed with the same parameters but different spectrum S/N. The results of the fitting procedure are presented in Fig. 14 where we plot the high velocity points of the spectroastrometric map in the vs. r ( ) diagram. The solid red lines represent the curves expected from the rotating disk model. In Fig. 15 we present the direct comparison between the observed spectroastrometric curve along the line of nodes (e.g. vs. ) with the best fit model (for S/N=100).
In all cases, for each ``test particle'', we are able to measure photocenter positions on the plane of the sky with an accuracy of better than 0 025, corresponding to 1/20 of the spatial resolution ( ).
The comparison between best fit and simulation parameters is shown in Table 2. In particular, we can recover the BH mass value with an accuracy of 0.1 dex. The M/L value is not well constrained by the models because with the adopted parameters the contribution to the gravitational potential of the stellar mass is negligible.
In the test cases considered here ( , stellar velocity dispersion of km s^{1} and Mpc), the radius of the BH sphere of influence is pc corresponding to , barely resolved with the adopted resolution. With the spectroastrometric technique, we have instead probed down to radii of the order of (1 pc), a factor 10 smaller than . This clearly demonstrates that with spectroastrometry we can recover information on spatial scales which are much smaller than the limit imposed by the spatial resolution, thus opening the possibility of probing the gravitational potential inside smaller BH spheres of influence.
Table 2: Fit results from the baseline model with different S/N.
Figure 16: Effect of varying the intrinsic flux distribution at subresolution scales: results for various intrinsic flux distributions. Top left panel: exponential flux distribution centered on with . Top right panel: exponential flux distribution centered on with . Bottom left panel: exponential flux distribution centered on with and a central hole of radius . Bottom right panel: exponential flux distribution centered on with . 

Open with DEXTER 
We now present a set of simulations where we vary the intrinsic flux distribution of the gas at subresolution scales and observe how this affects the recovery of the input model parameter values from the fit. As in Sect. 3.2.7 we adopt a basic model for the flux distribution that is an exponential function ( ) where the symmetry is center located at the position and the characteristic radius is set to . In Fig. 16 we show the cases of: the basic flux distribution, flux distribution with characteristic radius , flux distribution with central hole of radius , flux distribution with characteristic radius but centered at the position , all with a simulated noise for a S/N of 100. In Table 3 we report the best fit values of the free parameters. These simulations are chosen to represent extreme cases of sub resolution variation of the intrinsic flux distribution.
We can observe that with the fit we can recover the correct BH mass value with an accuracy of better than 0.16 dex, showing that even ``extreme'' variations of the flux distribution at subresolution scales, do not greatly affect the BH mass estimate.
Finally we consider a set of simulations where we decrease the value. In Fig. 17 we show the cases of BH masses and respectively, both with a simulated noise for a S/N of 100. In Table 4 we report the best fit values of the free parameters.
The case of lower black hole mass is particularly significative. We choose a value because, as observed in Sect. 3.2.8, the standard rotation curve method cannot distinguish between a case of BH and no BH. By looking at Fig. 17 we note that the BH signature is still present in a realistic spectroastrometric curve.
Table 3: Effect of varying the intrinsic flux distribution at subresolution scales: fit results.
Figure 17: Application of the method to low BH masses: results from simulations with S/N=100, ( left) and ( right). Top: spectroastrometric curves. Middle: 2D spectroastrometric maps. Bottom: fit results. 

Open with DEXTER 
Table 4: Application of the method to low BH masses: fit results.
Nonetheless, the fit is not well constrained by the few datapoints. Furthermore, we obtain a value of of (Table 4), substantially larger than the model mass of .
This is due to the fact that in our simulation the size of the velocity bins in the PVD diagram is 50 km s^{1}, matched to the assumed spectral resolution (). This value is insufficient to adequately resolve the emission line profile. As observed in Sect. 3.2.5, the spectroastrometric curves are stretched along the velocity axis with artificially high velocity values which result in an overestimated BH mass.
This problem can be overcome by increasing the spectral resolution of the simulated data to a higher value, , still well within reach of existing spectrographs. We then repeated the simulation setting the size of the velocity bin of the longslit spectra to 10 km s^{1}. In Fig. 18 we show the spectroastrometric curves and the result of the fit for this model. The sampling of the spectroastrometric curves is largely improved and the fit is consequentely better constrained. The derived BH mass value, is now in better agreement with the model value. The results of these simulations confirm that to assess the accuracy of BH masses from the spectroastrometric method one should first check that the line widths in the central emission region are well resolved spectrally.
In the last fits reported in Table 4 we hold fixed the masstoligth ratio M/L to the model value. The number of ''high velocity'' points in the spectroastrometric curves is small (only 2 for the model with 50 km s^{1} velocity bin) and it is not possible to constraint the radial dependence of the mass distribution. In general, when dealing with real data with very few ''high velocity'' points, one can only measure the total mass enclosed in the smaller radius which can be estimated from spectroastrometry. If it is possible to analyze the ``classical'' rotation curve, one could determine the masstoligth ratio M/L from that and use it in the spectroastrometric analysis. The fixed masstoligth ratios used in the above analysis should be interpreted as resulting from a classical analysis of the extended rotation curves.
6 Integral field spectroscopy
The extension of the spectroastrometric technique from multislit observations to Integral Field Units (IFU) is straightforward, and carries with it many advantages. With IFUs the complex issue of the limited spatial coverage of the slits (leading to the problem of velocity truncation) is partially removed and a full 2D map can be directly derived. Another substantial advantage is the reduced observing time requested to obtain a given S/N level, since there is no need to obtain observations of the same galaxy at different silt position angle.The analysis of IFU data now reduces to fitting a 2D Gaussian to each channel map in turn yielding the X, Y photocenters as a function of velocity, i.e. the spectroastrometric map. Armed with these derived photocenters one proceeds directly to the application of Sect. 5.
In practice it is important to determine the accuracy of each photocenter. This can easily be accomplished by using monte carlo realizations drawn from the original data in each channel giving a distribution of i values for Xand Y from which one can determine a median and interquartile spread.
In order to provide a quick analysis of advantages and disadvantages in the use of IFUs ve longslit spectra, we first recall the mandatory requirements needed to obtain an accurate and useful spectroastrometric curve. The spectra must be characterized by good spectral resolution, signaltonoise ratio and spatial sampling. The emission line profile must be well spectrally resolved otherwise, as discussed in Sect. 3.2.5, one would overestimate the BH mass. High signaltonoise and spatial sampling are required for a robust and accurate estimate of the spatial centroid of the emission line, as discussed in detail in Appendix A.
Compared to longslit spectrographs, IFUs have the advantage of a twodimensional spatial covering of the source which provides a more direct and accurate determination of the spectroastrometric map on the plane of the sky. However, for a given number of detector pixels, this is usually done at the expense of spectral resolution thus limiting the application of the spectroastrometric analysis to larger BHs. As a compromise between fieldofview and spectral resolution, an IFU can have a lower spatial sampling compared to longslit spectra which limits the accuracy of the photocenter determination. However, with IFUs it is possible to integrate longer, since one does not need to obtain spectra at different position angles.
As an example, let us compare the seeinglimited performances of ISAAC and SINFONI (both at the ESO VLT) for spectroastrometric use. ISAAC has a spatial sampling along the slit of , where for SINFONI in the seeing limited mode spatial pixels (spaxels) have angular dimensions of . There exist modes with finer spatial sampling ( ) but at the expense of a smaller field of view which might be insufficient in seeing limited observations. Overall, the spatial sampling in ISAAC and SINFONI are not extremely different, at least for moderately good seeing ( ). Regarding the spectral resolution, let consider the case of the K (J) band, at m (m). ISAAC, in the Short Wavelength Medium Resolution configuration offer a spectral resolution of 8900 (10 500) while SINFONI with the K (J) grating offers a spectral resolution of 4000 (2000). Clearly, the advantage of the 2D coverage provided by SINFONI is obtained at the expense of spectral resolution, and therefore SINFONI is not a good choice to detect small BHs (or large BHs in more distat galaxies).
In conclusions, IFUs have the advantage of the twodimensional spatial covering which allows longer integrations on source and more accurate spectroastrometric maps. On the other hand, longslit spectrographs can provide a better combination of high spectral resolution and spatial sampling for the detection of smaller BHs.
Regardless of the use of IFUs or longslit spectrographs, the superresolution provided by spectroastrometry allows BH mass determinations in galaxies at distances substantially higher than those that can be studied with the ``standard'' rotation curves method.
Figure 18: Application of the method to low BH masses: results from the simulations with BH mass of and velocity bins of 10 km s^{1}. Top: spectroastrometric curves (for simplicity we plot the error bars of the points only for the PA curve). Middle: spectroastrometric map. Lower panel: the result of the fit. 

Open with DEXTER 
7 Combining spectroastrometric and rotation curves
In this paper we have shown that by means of spectroastrometric curves we can recover information at scales smaller than the spatial resolution of the observations and we have provided a simple method to estimate BH masses. However, spectroastrometry is not a replacement of the standard rotation curve method. As shown in Sect. 3 (e.g. Fig. 1), the spectroastrometric curve is complementary to the rotation curve and it is clear that any kinematical model must account for both curves at the same time. For example, any contribution from extended mass distributions (e.g. stars) can be well constrained with rotation curves but not with spectroastrometric curves which sample very small scales. A detailed analysis of the constraints posed simultaneously by spectroastrometric and rotation curves will be discussed in forthcoming papers.
The reader might wonder why a combined modeling of spectroastrometric and rotation curves should be better than than modeling the full PositionVelocity diagram. The reason is very simple and is related to the fact that the full PVD depends on the unknown intrinsic flux distribution of the emission line. Marconi et al. (2006) showed that line profiles do depend on the assumed flux distribution while, on the contrary, such dependence is much weaker on the first moment of the line profile (i.e. the mean velocity). In this paper we have shown that the spectroastrometric curve depends on the assumed flux distribution in the low velocity range, while such dependence almost disappears in the high velocity range. Therefore, by selecting the rotation curve and the spectroastrometric curve we can greatly diminish the effects of the unknown intrinsic flux distribution which plague the full PVD.
8 Summary and conclusions
In this paper we have discussed the application of the spectroastrometric method in the context of kinematical studies aimed at measuring the masses of supermassive black holes in galactic nuclei and we have presented its advantages compared to classical ``rotation curves'' method.
We have conducted an extensive set of simulations which have shown that the presence of a supermassive black hole is revealed by a turnover in the spectroastrometric curve, with the high velocity component approaching a null spatial offset from the location of the galaxy nucleus. All the relevant information about the BH is encoded in the ``high velocity'' range of the spectroastrometric curve, which is almost independent of the spatial resolution of the observations. According to our simulations, the use of spectroastrometry can allow the detection of BHs whose apparent size of the sphere of influence is as small as 1/10 of the spatial resolution.
We have then provided a simple method to estimate BH masses from spectroastrometric curves. This method consist in the determination of the spectroastrometric map, that is the positions of emission line photocenters in the available velocity bins. This is trivially obtained from IFUs but can also be obtained by combining longslit spectra centered on a galaxy nucleus and obtained at different position angles. From the ``high velocity'' points in the spectroastrometric map one can then obtain a rotation curve to trivially estimate the BH mass. We have provided practical applications of this method to simulated data but considering noise at different levels. From this analysis, we confirm that with seeing limited observations ( ) we are able to detect a BH with mass , where D is the galaxy distance. This is a factor 10 better than can be done with the classical method based on rotation curves.
Finally, we have discussed advantages and disadvantages of IFU vs. longslit spectrographs concluding that IFUs have the advantage of the twodimensional spatial covering which allows longer integrations on source and more accurate spectroastrometric maps. On the other hand, longslit spectrographs can provide a better combination of high spectral resolution and spatial sampling for the detection of smaller BHs. Regardless of the adopted type of spectrograph, the superresolution provided by spectroastrometry allows BH mass determinations in galaxies at distances substantially higher than those that can be studied with the ``standard'' rotation curves method.
Appendix A: The measure of the light profile centroid
The first fundamental step in the application of the spectroastrometric
method is to find the best possible way to estimate the centroid
position of the light profile for a given wavelength (or velocity).
In principle one could simply estimate the centroid position as the
weighted mean of the position with the corresponding emission line
flux. If y_{i}
denotes the position and I_{i}
the flux of the ith pixel along the slit (for a
given velocity), the centroid position
is
simply:
The problem in using Eq. (A.1) resides in the presence of noise and in the fact that the value of depends on the choice of the position range used to compute the weighted mean. In Fig. A.1 we exemplify the effect of choosing different ranges to calculate the centroid of real data. This effect will not be negligible especially because one is looking for precisions which are much lower than the spatial resolution (e.g., 2.5 pixels for the data shown in the figure).
Figure A.1: Example of a light profile along the slit for a given velocity bin for real data. Three different position ranges are indicated as horizontal lines and the corresponding centroid positions are indicated with vertical lines with the same style. 

Open with DEXTER 
Another possibility to estimate the centroid position is provided by a
parametric fitting of the light profile with an appropriate function.
For instance, we have tried a GaussHermite expansion:
where and H_{3} and H_{4} are the GaussHermite polynomials of the third and fourth order.
The advantage of the fitting method is that with the parametric fit the results are much less sensitive to noise and choice of the position range for the analysis but it is difficult to find an appropriate parametric function which can reproduce the profile shape without a large increase of the free parameters which, in turn, decreases the accuracy of the centroid position. The light profile along the slit has often a complex shape, sometimes with multiple peaks and there is not any physical reason to suggest a particular parametric function.
Moreover, the shape of the light profile for a given velocity bin is determined not only by velocity field of the rotating gas disk but also by many other factors, like the intrinsic light distribution of the emission line, the slit width and position and many others. Therefore one has to investigate what is the signature which is most related to the velocity field.
The piece of information we seek is the average position of the gas rotating at a given observed velocity which can be converted into the gravitational potential in the nuclear region, with the usual assumption of a thin, circularly rotating disk.
As explained above, the use of the entire light profile along the slit for a given velocity bin to calculate the photocenter position is extremely complex, moreover such position might be influenced by many features which do are not connected with the rotation of the disk. After testing the method with real and simulated data, we have reached the conclusion that the quantity we need to measure is not the centroid of the whole light profile but only of the central spatially unresolved component which is related to gas located very close to the galaxy nucleus (e.g. gas rotating around the BH). Thus, the problem of measuring the light centroid along the slit can be transformed into measuring the position of the peak of the central unresolved emission.
Figure A.2: Example of centroid determination on a real spectrum. The dashed line denotes the fitted Gausshermite functions and the relative centroid position (vertical dashed line). The solid lines denotes the Gaussian function which has been fitted in a window around the principal peak as described in the text. The vertical solid lines represent the edges of the window and the Gaussian center respectively. 

Open with DEXTER 
The profile peak at the location of the nucleus is due to the spatially unresolved emission from the innermost parts of the rotating gas disk. Because it is spatially unresolved, the width of central peak is expected to be of the order of the spatial resolution, and its shape can be well approximated with a Gaussian function. This component is easy to identify and fit. Moreover this unresolved component is emitted from the innermost parts of the gas disk, where the gravitational potential of the BH is strongest with a larger influence on the gas kinematics. Conversely, the spatially resolved features are very complex to fit and are influenced by many factors not directly connected with the gravitational potential.
In practice, we first obtain an initial guess of the principal peak position and we select the light profile within a box centered on this initial guess and whose width is of the order of the spatial resolution. We then accurately locate the peak position by fitting a simple Gaussian function to this selected portion of the light profile. In details, we usually set the width of the box as where is the spatial resolution (standard deviation of the Gaussian approximating the point spread function) and k a constant between 1 and 3. Using values of k too close to 1 results in too few fit points but k should also not be too large if we want to consider only the unresolved component. From our tests, the best choice has come out to be .
This method used to obtain the light centroid provides us with an important control parameter, that is the width of the fitted Gaussian: if this is much larger than the spatial resolution, this is an indication that line emission is dominated by a spatially extended component and the obtained value of the centroid is less reliable.
In Fig. A.2 we show an example of centroid determination obtained by fitting a Gausshermite function to the whole profile and a simple Gaussian in a wide window centered around the principal peak. The centroid values obtained with the two methods might or might not be consistent but it is clear that the centroid value with the first method is shifted to the right by the contribution of the profile wings, which, moreover, are not even well reproduced by the Gausshermite function. In the bottom panel of Fig. A.2 none of the two methods can well reproduce the profile peak, but the width of the fitted Gaussian turns out to be 2.5 times larger than the spatial resolution indicating that the principal peak is partially resolved, and the value of the centroid is less reliable.
Figure A.3: FWHM of the Gaussian fitted to the principal peak of the light profile for spectroastrometric curves of Fig. 11. Solid line: PA1 curve. Dotted line: PA2 curve. Dotdashed lines: PA3 curve. The horizontal dashed line denotes 1.1 times the PSF FWHM. The gray dashed lines are the estimated limits of the ``high velocities''. 

Open with DEXTER 
Appendix B: Building the spectroastromertic map from multiple slits observations
In this section we describe how to build a 2D spectroastrometric map, i.e. as to estimate the photocenter position on the plane of the sky for each velocity bin, from multiple slits observations.
We first consider a reference frame in the plane of the sky
which is centered
on the center of PA1 slit and which has the X axis
along the North direction.
For a given velocity bin v_{i}
the position of the light centroid on the sky
plane is [x(v_{i}),
y(v_{i})];
the expected photocenter position along a given
slit (i.e. a point of the spectroastrometric curve) is simply
the position
[x(v_{i}),
y(v_{i})]
projected along the slit direction that is:
where and are the coordinates of the slit centre and is the position angle of the slit referred to the X axis (i.e. the North direction as in the usual definition of PA).
The spectroastrometric curve provides the measured centroid
position along a
given slit .
Thus, in order to determine the free
parameters
(
is
known), one can minimize the following :
It is worth noticing that the relevant quantities are not the absolute positions of the slit centers, but the relative ones. Indeed we have chosen the center of the reference frame coincident with the center of the PA1 slit. The number of free parameters is then 6 (2 slit centers and the photocenter position on the sky, i.e. 6 coordinates) compared to 3 data points (the photocenter positions from the 3 slits). However, the problem is not undefined since many velocity bins are available and the slit center positions must be the same for all velocity bins. Therefore, for each set of slit center positions ( ) we minimize separately all . The best values are then those which the minimize
if N points are available from each spectroastrometric curve, the number of degrees of freedom is then , where 4 is the number of unknown slit center coordinates and 2N is the number of unknown photocenter positions. Since N is larger than 4 the problem is well posed. The final spectroastrometric map on the plane of the sky is that given by the best fitting set of slit centers.
The sum over v_{i} in (B.3) is extended only over the ``high velocity'' range. Indeed, in Sect. 3.2 we concluded that the ``high velocity'' range of the spectroastrometric curve is more robust, and less affected by slit losses which artificially change the photocenter position in the low velocity range. As explained in Appendix A, where we discuss in detail, centroid determination, we find the spectroastrometric curve by fitting a Gaussian to the principal peak of line emission along the slit (that centered on the nucleus position). Then the ``high velocity range'' is selected by imposing that the FWHM of the fitted Gaussian is lower than 1.1 times the spatial resolution (FWHM of the PSF). This ensures that we consider oinly velocities where the line emission along the slit is spatially unresolved.
References
 Aime, C., Borgnino, J., Lund, G., & Ricort, G. 1988, in European Southern Observatory Astrophysics Symposia, ed. F. Merkle, 29, 249 [Google Scholar]
 Bailey, J. 1998, MNRAS, 301, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Baines, D., Oudmaijer, R. D., Mora, A., et al. 2004, MNRAS, 353, 697 [NASA ADS] [CrossRef] [Google Scholar]
 Beckers, J. M. 1982, Optica Acta, 29, 361 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J. T. S. 1987, Galactic Dynamics (Princeton University Press) [Google Scholar]
 Brannigan, E., Takami, M., Chrysostomou, A., & Bailey, J. 2006, MNRAS, 367, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Christy, J. W., Wellnitz, D. D., & Currie, D. G. 1983, Lowell Observatory Bulletin, 9, 28 [NASA ADS] [Google Scholar]
 Marconi, A., Pastorini, G., Pacini, F., et al. 2006, A&A, 448, 921 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moorwood, A., Cuby, J.G., Ballester, P., et al. 1999, The Messenger, 95, 1 [NASA ADS] [Google Scholar]
 Pontoppidan, K. M., Blake, G. A., van Dishoeck, E. F., et al. 2008, ApJ, 684, 1323 [NASA ADS] [CrossRef] [Google Scholar]
 Porter, J. M., Oudmaijer, R. D., & Baines, D. 2004, A&A, 428, 327 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Porter, J. M., Oudmaijer, R. D., & Baines, D. 2005, in The Nature and Evolution of Disks Around Hot Stars, ed. R. Ignace, & K. G. Gayley, ASP Conf. Ser., 337, 299 [Google Scholar]
 Takami, M., Bailey, J., & Chrysostomou, A. 2003, A&A, 401, 655 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Whelan, E. T., Ray, T. P., Bacciotti, F., et al. 2005, Nature, 435, 652 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
All Tables
Table 1: Lines of Nodes and BH position estimates from noisy data.
Table 2: Fit results from the baseline model with different S/N.
Table 3: Effect of varying the intrinsic flux distribution at subresolution scales: fit results.
Table 4: Application of the method to low BH masses: fit results.
All Figures
Figure 1: Example of the measurements on which the construction of rotation and spectroastrometric curves is based. Upper left panel: the longslit spectrum of a given continuumsubtracted emission line (isophotes) with superimposed two particular sections at fixed slit position and velocity (dotted lines). Bottom left panel: the profile relative to the section at fixed slit position with superimposed the position of its centroid (dashed line). Upper right panel: the profile relative to the Sect. at fixed velocity with superimposed the position of its centroid (dashed line). Bottom right panel: the positionvelocity diagram (PVD) of the upper left panel with superimposed rotation (blue dots with error bars) and spectroastrometric curves (red points with error bars). In all panels velocities are measured relative to the systemic velocity of the galaxy which is set to the zero point of the velocity scale. 

Open with DEXTER  
In the text 
Figure 2: Spectroastrometric curve (black solid line) obtained from a simulated spectrum for a rotating gas disk with various values of the central BH and stellar mass. In gray we plot the isophotes of the simulated line spectrum. Upper left panel: case of BH without mass contribution from stars. Upper right panel: stars only, with . Lower panels: ( left) and ( right) with included the same stellar mass distribution as in model in the upper right panel. Note that the differences in the kinematics observed in the left panels (both BH dominated, but with and without stellar contributions) are marginal, and only visible i.e. beyond from the nucleus. 

Open with DEXTER  
In the text 
Figure 3: Left panel: spectroastrometric curves from spectra differing only in spatial resolution. Solid line: spatial resolution of . Dotted line: spatial resolution of . Dashed line: spatial resolution of . The vertical longdashed lines denote the ``high velocities''. Right panel: corresponding rotation curves. 

Open with DEXTER  
In the text 
Figure 4: Left panel: spectroastrometric curves for spectra differing in the value BH mass. Solid line: . Dotted line: . Right panel: rotation curves for the same models. 

Open with DEXTER  
In the text 
Figure 5: Left panel: simulated spectroastrometric curves from spectra differing only in slit position angle (referred to the disk line of nodes ). Right panel: measured amplitude of the spectroastrometric curve (distance between the maximum and the minimum) as a function of the slit position angle referred to the disk line of nodes; the dashed line denotes a cosine function opportunely rescaled to match the value of the amplitude of . 

Open with DEXTER  
In the text 
Figure 6: Left panel: isovelocity contour map of the line of sight velocity field on the plane of sky for a rotating gas disk with the of line of nodes along the x axis, with superposed a wide slit. Right panel: comparison of spectroastrometric curves for simulated spectra differing only in slit width (we remind that we assumed a spatial resolution of 0 ). Solid line: wide slit. Dotted line: wide slit. Dashed line: wide slit. 

Open with DEXTER  
In the text 
Figure 7: Comparison of spectroastrometric curves for spectra differing only in spectral resolution. Solid line: spectral resolution of 10 km s^{1}. Dotted line: spectral resolution of 50 km s^{1}. Dashed line: spectral resolution of 100 km s^{1}. Dotdashed line: spectral resolution of 150 km s^{1}. 

Open with DEXTER  
In the text 
Figure 8: Comparison of spectroastrometric curves for spectra differing only in intrinsic velocity dispersion. Solid line: intrinsic velocity dispersion of 10 km s^{1}. Dotted line: intrinsic velocity dispersion of 50 km s^{1}. Dashed line: s intrinsic velocity dispersion of 100 km s^{1}. Dotdashed line: intrinsic velocity dispersion of 150 km s^{1}. 

Open with DEXTER  
In the text 
Figure 9: Comparison of spectroastrometric curves for spectra differing only in intrinsic flux distribution. Solid line: base model flux distribution (explained in the text). Dotted line: model with . Dashed line: model with . Dotdashed line: base model with a central hole of radius . 

Open with DEXTER  
In the text 
Figure 10: Comparison of spectroastrometric curves for spectra differing only in intrinsic flux distribution. Solid line: base model flux distribution (explained in the text). Dotted line: base model centered on position. Dashed line: base model centered on position. Dotdashed line: base model centered on position. 

Open with DEXTER  
In the text 
Figure 11: Spectroastrometric curves for three simulated spectra of the same model with different PA of the slit. Solid line: slit with PA (referred to the disk line of nodes position angle) . Dotted line: slit with . Dashed line: slit with . 

Open with DEXTER  
In the text 
Figure 12: Spectroastrometric 2D map derived from the three spectroastrometric curves of Fig. 11. Left panel: derived photocenter positions on the sky plane, the black open diamonds are the points actually used for the minimization. Right panel: the 3D plot of the map, where the z axis is velocity. 

Open with DEXTER  
In the text 
Figure 13: Left row: spectroastrometric curves for the model of Fig. 11 with added noise for a S/N of 20, 50, and 100 ( from top to bottom). Solid line: slit at PA = . Dotted line: slit at PA . Dashed line: slit at PA = . For simplicity we plot the error bars of the points only for the PA curve. Right row: spectroastrometric 2D map derived from the three spectroastrometric curves of the left row. The red filled circle marks the derived position of the BH, while the red solid line represents the position angle of the line of nodes ( ) obtained by fitting a straight line to the XY position of the high velocity points. 

Open with DEXTER  
In the text 
Figure 14: Fit results for models of Fig. 13 with S/N=20 ( top panel), S/N=50 ( middle) and S/N=100 ( bottom). Points with error bars are the result of ``observations''. The red line is the best fit model and the dotted line is the expected contribution from the stellar mass. 

Open with DEXTER  
In the text 
Figure 15: Observed and model spectroastrometric curves at S/N=100 (e.g. vs. ). The horizontal dotted line is the BH ``S'' coordinate S_{0}. The vertical dotted line is the systemic velocity . 

Open with DEXTER  
In the text 
Figure 16: Effect of varying the intrinsic flux distribution at subresolution scales: results for various intrinsic flux distributions. Top left panel: exponential flux distribution centered on with . Top right panel: exponential flux distribution centered on with . Bottom left panel: exponential flux distribution centered on with and a central hole of radius . Bottom right panel: exponential flux distribution centered on with . 

Open with DEXTER  
In the text 
Figure 17: Application of the method to low BH masses: results from simulations with S/N=100, ( left) and ( right). Top: spectroastrometric curves. Middle: 2D spectroastrometric maps. Bottom: fit results. 

Open with DEXTER  
In the text 
Figure 18: Application of the method to low BH masses: results from the simulations with BH mass of and velocity bins of 10 km s^{1}. Top: spectroastrometric curves (for simplicity we plot the error bars of the points only for the PA curve). Middle: spectroastrometric map. Lower panel: the result of the fit. 

Open with DEXTER  
In the text 
Figure A.1: Example of a light profile along the slit for a given velocity bin for real data. Three different position ranges are indicated as horizontal lines and the corresponding centroid positions are indicated with vertical lines with the same style. 

Open with DEXTER  
In the text 
Figure A.2: Example of centroid determination on a real spectrum. The dashed line denotes the fitted Gausshermite functions and the relative centroid position (vertical dashed line). The solid lines denotes the Gaussian function which has been fitted in a window around the principal peak as described in the text. The vertical solid lines represent the edges of the window and the Gaussian center respectively. 

Open with DEXTER  
In the text 
Figure A.3: FWHM of the Gaussian fitted to the principal peak of the light profile for spectroastrometric curves of Fig. 11. Solid line: PA1 curve. Dotted line: PA2 curve. Dotdashed lines: PA3 curve. The horizontal dashed line denotes 1.1 times the PSF FWHM. The gray dashed lines are the estimated limits of the ``high velocities''. 

Open with DEXTER  
In the text 
Copyright ESO 2010