Detecting supermassive black holes in ultracompact dwarf galaxies

One fundamental question in astronomy remains that of how galaxies form and evolve. Cosmological simulations and large scale surveys of the past two decades have started to test formation mechanism theories. While the combination of these tools have confirmed predictions of the large scale structure of the Universe, they have exposed a number of issues on smaller scales, especially in the low-mass galaxy regime. It is unclear whether these issues are associated with the simulations or observational constraints. Therefore, detailed studies of low-mass stellar systems are of vital importance in understanding this conflict. It is the goal of this dissertation to provide such information from objects known as ultracompact dwarf galaxies (UCDs). UCDs are a class of compact stellar systems intermediate to that of globular clusters and dwarf galaxies. Recent analyses have shown that the dynamical mass appears to be systematically elevated when compared to the stellar population estimates, which can be explained by central black holes (BH) making up 10-15% of the total mass. In my dissertation, I test this hypothesis in two projects searching for central BHs in three UCDs with the Jeans Anisotropic Modeling (JAM) and Schwarzschild dynamical modeling methods. My first project presents the JAM modeling of two Virgo UCDs, VUCD3 and M59cO. Assuming isotropy, I detect central BHs in both objects with masses of 4:4+2:5 -3:0 X 106M. in VUCD3 and 5:8+2:5 - 2:8 x 106M. in M59cO. These BHs make up an astonishing ~13% and ~18% of the total mass for VUCD3 and M59cO, respectively. My second project presents dynamical modeling using both models, mentioned above, of the most massive UCD, M59-UCD3. The best fit models in the JAM and an axisymmetric Schwarzschild model have BHs ranging from 2:5 􀀀 5:9  106 M., while a triaxial Schwarzschild model prefers a solution with no BH. However, models with a BH in each technique provide a better fit to the central VRMS profile, and thus I estimate the BH mass to be 4:2+2:1 1:7 x 106 M. These BHs, combined with structural information of each UCD, suggest all three objects are the tidally stripped remnants of ~ 109 - 1010M. galaxies.

DETECTING SUPERMASSIVE BLACK HOLES IN ULTRACOMPACT DWARF GALAXIES by Christopher Ahn A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Physics Department of Physics and Astronomy The University of Utah August 2018 Copyright c Christopher Ahn 2018 All Rights Reserved The University of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL This dissertation of Christopher Ahn has been approved by the following supervisory committee members: Anil C. Seth , Chair May 2, 2018 Date Approved Kyle Dawson , Member May 2, 2018 Date Approved , Member Zheng Zheng May 2, 2018 Date Approved , Member Stephan LeBohec May 2, 2018 Date Approved Miriah Meyer , Member May 2, 2018 Date Approved and by of Peter Trapa the Department of Physics and Astronomy and by David B. Kieda of the Graduate School. , Chair ABSTRACT One fundamental question in astronomy remains that of how galaxies form and evolve. Cosmological simulations and large scale surveys of the past two decades have started to test formation mechanism theories. While the combination of these tools have confirmed predictions of the large scale structure of the Universe, they have exposed a number of issues on smaller scales, especially in the low-mass galaxy regime. It is unclear whether these issues are associated with the simulations or observational contraints. Therefore, detailed studies of low-mass stellar systems are of vital importance in understanding this conflict. It is the goal of this dissertation to provide such information from objects known as ultracompact dwarf galaxies (UCDs). UCDs are a class of compact stellar systems intermediate to that of globular clusters and dwarf galaxies. Recent analyses have shown that the dynamical mass appears to be systematically elevated when compared to the stellar population estimates, which can be explained by central black holes (BH) making up ∼10-15% of the total mass. In my dissertation, I test this hypothesis in two projects searching for central BHs in three UCDs with the Jeans Anisotropic Modeling (JAM) and Schwarzschild dynamical modeling methods. My first project presents the JAM modeling of two Virgo UCDs, VUCD3 and M59cO. 6 Assuming isotropy, I detect central BHs in both objects with masses of 4.4+2.5 −3.0 × 10 M 6 in VUCD3 and 5.8+2.5 −2.8 × 10 M in M59cO. These BHs make up an astonishing ∼13% and ∼18% of the total mass for VUCD3 and M59cO, respectively. My second project presents dynamical modeling using both models, mentioned above, of the most massive UCD, M59-UCD3. The best fit models in the JAM and an axisymmetric Schwarzschild model have BHs ranging from 2.5 − 5.9 × 106 M , while a triaxial Schwarzschild model prefers a solution with no BH. However, models with a BH in each technique provide a better fit to the central VRM S profile, and thus I estimate the BH mass 6 to be 4.2+2.1 −1.7 × 10 M . These BHs, combined with structural information of each UCD, suggest all three objects are the tidally stripped remnants of ∼ 109 − 1010 M galaxies. CONTENTS ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii CHAPTERS 1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Overview of Galaxy Formation and Evolution . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Black Holes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1 Detection Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Black Hole Scaling Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3 Dynamical Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3.1 Ingredients for Dynamical Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3.2 JAM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.3.3 Schwarzschild Dynamical Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.4 Outline of the Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2. DETECTION OF SUPERMASSIVE BLACK HOLES IN TWO VIRGO ULTRACOMPACT DWARF GALAXIES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Observations and Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 HST Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Gemini NIFS data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3 Kinematic Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Creating a Mass Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Dynamical Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 Isotropic Model Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.2 Impact of Anisotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3. 15 18 18 19 20 24 30 33 34 39 40 THE BLACK HOLE IN THE MOST MASSIVE ULTRACOMPACT DWARF GALAXY M59-UCD3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Data and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Imaging Data and Deriving a Mass Model . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Kinematic Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Dynamical Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Jeans Anisotropic Models (JAM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 47 47 55 58 60 3.3.2 Axisymmetric Schwarzschild Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Triaxial Schwarzschild Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.4 Summary of Dynamical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Radio and X-ray Observations of UCDs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Radio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 X-ray . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.3 Fundamental Plane of BH Accretion . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 Summary of Main Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Implications of a Central Massive BH . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.3 Progenitor Galaxy Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.4 Further Evidence and Future Prospects for SMBHs in UCDs . . . . . . . . 4. 67 68 71 73 73 74 77 77 78 79 81 82 CONCLUSIONS AND FUTURE PROSPECTS . . . . . . . . . . . . . . . . . . . . 84 4.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.2 Future: Galaxy Cluster Assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.3 Future: Scaling Relations and UCD Formation . . . . . . . . . . . . . . . . . . . . . . . 85 REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 v LIST OF FIGURES 1.1 Examples of the four major galaxy types with the name and image courtesy provided below each image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Correlation between the mass of a central black hole and the line-of-sight velocity dispersion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3 Example of one-dimensional MGE fit to R−2 profile. . . . . . . . . . . . . . . . . . . . . 10 2.1 The two galaxy-UCD systems discussed in this paper. . . . . . . . . . . . . . . . . . . . 17 2.2 Integrated spectra of VUCD3 (left) and M59cO (right). . . . . . . . . . . . . . . . . . . 22 2.3 Surface brightness profile of VUCD3 (left) and M59cO (right) in HST filters used for dynamical modeling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.4 The color profiles of VUCD3 (left) and M59cO (right) . . . . . . . . . . . . . . . . . . . 27 2.5 Cumulative likelihood, assuming isotropy, of the BH mass in VUCD3 (left) and M59cO (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.6 Dispersion profiles of VUCD3 (left) and M59cO (right) . . . . . . . . . . . . . . . . . . 34 2.7 Contour plots showing the degeneracy between βz , MBH , and Γ . . . . . . . . . . . 36 2.8 The dynamical-to-stellar mass ratio Γ vs. total dynamical mass. . . . . . . . . . . . 38 2.9 Dispersion profiles of VUCD3 (left) and M59cO (right) . . . . . . . . . . . . . . . . . . 40 3.1 The M59/M59-UCD3 system discussed in this paper. . . . . . . . . . . . . . . . . . . . . 46 3.2 Surface brightness profile of M59-UCD3 in HST/F814W . . . . . . . . . . . . . . . . . 48 3.3 Color profile of M59-UCD3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.4 Integrated spectrum of M59-UCD3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.5 Full kinematic measurements of M59-UCD3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.6 MCMC post burn-in phase distributions for our default model. . . . . . . . . . . . . 62 3.7 Contour plots showing all three dynamical modeling techniques. . . . . . . . . . . . 63 3.8 VRM S comparison between smoothed data (left) and best-fit JAM model (right). 63 3.9 Cumulative likelihood of MBH from the JAM modeling . . . . . . . . . . . . . . . . . . 64 3.10 VRM S comparison between the data, best-fit default JAM model . . . . . . . . . . . 66 3.11 Black points show a rectangular aperture along the semi major axis. . . . . . . . . 72 3.12 VLA images of three massive UCDs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 3.13 The fundamental plane for UCDs with known BHs . . . . . . . . . . . . . . . . . . . . . 78 3.14 Dynamical-to-stellar mass ratio Γ vs. total dynamical mass. . . . . . . . . . . . . . . 80 LIST OF TABLES 2.1 Resolved Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.2 Best-fit Sérsic Parameters of VUCD3 and M59cO . . . . . . . . . . . . . . . . . . . . . . 25 2.3 Multi-Gaussian Expansions (MGEs) Used as the Default Mass and Luminosity Models in Dynamical Modeling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.1 Best-fit Sérsic Parameters of M59-UCD3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.2 Default Multi-Gaussian Expansion (MGE) Used in the Dynamical Modeling. . 54 3.3 Gemini/NIFS Kinematic Measurements of M59-UCD3 . . . . . . . . . . . . . . . . . . . 59 3.4 Calculations of Triaxial Schwarzschild Model Reduced χ2 Independently . . . . . 70 3.5 Chandra X-ray Constraints and Measurements . . . . . . . . . . . . . . . . . . . . . . . . . 75 CHAPTER 1 INTRODUCTION When we look up to the night sky we see, on average, a few thousand stars. However, all of these stars only make up a small fraction of the total number of stars within our galaxy, the Milky Way. Furthermore, the Milky Way is only one of billions of galaxies in the Universe. Most of these galaxies are one of four major galaxy types. Based on their appearance, these galaxy types are elliptical, spiral, lenticular, and irregular, as shown in Figure 1.1. Galaxies are far more complex than just their appearance. These galaxies can be subdivided based on a wide range of characteristics such as their mass, size, and luminosity. Furthermore, the boundary between what is classified as a galaxy versus a star cluster is unclear. For example, some globular clusters (GCs) appear to be the central remnants of once much larger objects. One of the biggest questions in astronomy to this day remains how galaxies and their constituents formed and evolved from one type to another. The time scale on which this occurs is so long that it is only possible for us to see a snapshot, which is effectively frozen in time. This gives astronomers two general approaches to attempt to understand the Universe as a whole: 1.) observe as many nearby objects as possible with the widest range of characteristics, or 2.) observe objects that are farther away (at “high redshift”), as we can look back in time, due to the finite speed of light, which has had to travel to reach us. Nearby objects are generally brighter and spatially resolved, which allows us to study them in great detail. For example, modern telescopes and instruments enable us to observe the properties of their stars, gas, and dust in multiple positions at once. This information can be used to infer the interplay between various internal components building up to the global properties of galaxies. The results of these detailed studies also provide us with information about objects at high redshift, at a time when these objects were just beginning to form and evolve, and where we are limited to the global properties of the system. 2 M87: J.-C. Cuilandre and G. Anselmi– Canada-France-Hawaii Telescope (CFHT)/Coelum M101: CFHT/NOAO/AURA/NSF/NASA/ESA/STScI NGC 4753: ARC and the SDSS collaboraPon, www.sdss.org NGC 4449: NASA/STScI Figure 1.1: Examples of the four major galaxy types with the name and image courtesy provided below each image. From top left to bottom right these are elliptical (M87), spiral (M101), lenticular (NGC 4753), and irregular (NGC 4449). 3 1.1 Overview of Galaxy Formation and Evolution There are two main theories of galaxy formation: the top-down theory and bottom-up theory. Both of these theories attempt to explain the formation of spiral type galaxies, while lenticular and elliptical galaxies have been shown to be a product of major mergers of spiral galaxies (Toomre & Toomre 1972). In the top-down theory of galaxy formation, galaxies form through direct collapse of large gas clouds (Eggen et al. 1962). The matter content of the early Universe consisted of mostly dark matter concentrated in large clumps, which interacted gravitationally. These interactions created tidal torques that pulled the clumps into a slow rotation. As the Universe expanded, the baryonic matter began to cool and the gas clouds, within the dark matter clumps, began to interact through collisions. These collisions acted to dissipate some of the energy of the gas cloud, which contracted and fell toward the center of the dark matter clumps. Through conservation of angular momentum, these gas clouds began to rotate faster as they contracted, which caused matter to form in tight disks. In the bottom-up theory of galaxy formation, matter started out in smaller dark matter clumps, which merged to form galaxies (Searle & Zinn 1978). Since the baryonic matter started out in small dark matter clumps and interacted in the same way as the large dark matter clumps, this theory still produces disk-like structure formation. This theory, commonly referred to as hierarchical merging, matches observations of a large number of small galaxies relative to large ones. Therefore, hierarchical merging has become widely accepted as the dominant galaxy formation mechanism. Over the past two decades, galaxy simulations and large scale galaxy surveys have started to test these formation mechanism theories (e.g., Colless et al. 2001; Springel et al. 2005; Vogelsberger et al. 2014). While these simulations and large scale surveys have confirmed many hierarchical merging predictions of the large scale structure of the Universe, they have also exposed a number of issues on smaller scales. Some of these issues include the following: 1.) The dwarf galaxy problem that refers to the number of observed dwarf galaxies, relative to normal galaxies, is an order of magnitude too low (Klypin et al. 1999; Moore et al. 1999). One proposed solution to this problem is that there may exist a large number of dwarf galaxies that astronomers haven’t been able to observe because their stars have been stripped via tidal interactions. However, this solution creates a new problem, the so-called “too big to fail” problem, which states that many of the satellites predicted by cosmological simulations are so massive that there’s no way they could not have visible stars. 2.) The cusp-core problem that refers to the discrepancy between the centers of galaxies in simulated 4 dark matter haloes showing the dark matter density rising steeply near the center (cuspy), and the rotation curves of observed dwarf galaxies showing flat central dark matter densities (cored; Moore 1994). Solutions to this problem have been suggested in the form of feedback mechanisms from supernovae and active galactic nuclei (AGN), which act to smooth out the central dark matter density distribution. Thus far, it is unclear whether these issues are associated with the simulations or observational constraints. Therefore, it is important that observational astronomers continue to provide detailed, quantitative information on the intrinsic structure of as many galaxies as possible, especially in the low-mass regime where simulations seem to break down. It is the goal of this dissertation to provide these types of measurements from a special class of objects known as ultracompact dwarf galaxies (UCDs). These UCDs are compact stellar systems with masses and radii of M > 2 × 106 M and r > 10 pc. In the mass/luminosity vs. size plane, UCDs fall directly between GCs and dwarf galaxies. Prior to their discovery in the late 1990s, this region in parameter space was bare and GCs are relatively dim, which restricted detailed studies in this mass regime to the nearest galaxies. With the discovery of UCDs, which are generally larger and brighter than GCs, astronomers have only recently started to investigate low-mass stellar systems at the distance of the nearest galaxy clusters. Early interpretations of UCD formation suggested that they could either be the high mass end of the GC distribution (e.g., Mieske et al. 2002; Fellhauer & Kroupa 2002, 2005; Kissler-Patig et al. 2006; Murray 2009), or the tidally stripped remnants of dwarf galaxies with progenitor masses of the order of 109 − 1010 M (e.g., Bekki et al. 2001, 2003; Drinkwater et al. 2003; Strader et al. 2013; Forbes et al. 2014; Seth et al. 2014). However, since there is a smooth transition between these three classes of objects, UCDs are likely a combination of both (e.g., Brodie et al. 2011; Da Rocha et al. 2011; Misgeld & Hilker 2011; Norris et al. 2014; Janz et al. 2016). Therefore, UCDs provide an important test of the formation and evolution of low-mass galaxies. 1.2 Black Holes Black holes (BH)s are ubiquitous in galaxies. Current measurements have placed BHs into three broad categories based on their mass: stellar mass BHs, intermediate mass black holes (IMBHs), and supermassive black holes (SMBHs). The latter of which has been shown to correlate with a number of host galaxy properties. This suggests SMBHs play an important role in galaxy evolution. The combination of stellar observations and evolutionary models provides us with a 5 pretty good idea of how stellar mass BHs are formed. Main sequence stars with M & 8 M quickly exhaust their fuel source for fusing elements together and collapse in violent supernova explosions. At this end stage of stellar evolution, if the mass of remnant core exceeds the maximum mass sustainable with neutron degeneracy pressure, the star collapses to form a BH with a mass ranging from ∼ 5 − 100 M . IMBHs are black holes with masses ranging from ∼ 100−1×106 M . Their formation is thought to occur through the merging of stellar mass BHs (e.g., Abbott et al. 2016b) or from the collision of massive stars in dense cluster environments. Although most astronomers believe they must exist, the search to confirm this has been plagued with controversy. For example, IMBH detections have been claimed in ω Cen (e.g., Noyola et al. 2010; Baumgardt 2017), M54 (Ibata et al. 2009), and the Andromeda GC G1 (Gebhardt et al. 2005), but the mass measurements have been inconsistent and the search for accretion evidence has been unsuccessful (e.g., van der Marel & Anderson 2010; Miller-Jones et al. 2012; Haggard et al. 2013). Finally, SMBHs are black holes with masses M > 106 M . The formation of SMBHs is an active field of research in astronomy, in large part, due to the apparent existence of SMBHs in nearly all massive galaxies. For example, the Milky Way hosts a SMBH with M ∼ 4 × 106 M (Genzel et al. 1997; Ghez et al. 1998). Furthermore, observations of distant quasars, which are SMBHs that are actively accreting matter, have implied that most massive galaxies should host SMBHs in the present day (Soltan 1982; Marconi et al. 2004). Astronomers generally agree that once a BH settles into the center of a galaxy, it can grow through the accretion of matter. The location at the centers of galaxies is not surprising due to dynamical friction. Dynamical friction is the loss of momentum and energy of moving objects through gravitational interactions. As a massive object moves through a cloud of smaller objects, the gravitational interactions accelerate the smaller objects, while the massive object loses momentum and kinetic energy, allowing it to sink toward the center of the cloud. However, the origin of the initial “seeds” of SMBHs remains an open question. Possible explanations have ranged from the direct collapse of large gas clouds (Begelman et al. 2006) to rapidly accreting remnants of the first massive stars (Whalen & Fryer 2012). To distinguish between the growing number of explanations, astronomers need to continue to measure and study these exotic objects. 1.2.1 Detection Techniques Detecting and measuring the mass of BHs in galaxies is often very difficult. By definition, BHs do not emit any light. Therefore, to find them, astronomers rely on their influence on 6 the objects around them. There are three main techniques used to detect BHs: modeling the gravitational influence of BHs on stars and gas around them, from the radiation from matter accreting onto a BH, and from the emission of gravitational waves. Modeling the gravitational influence of BHs on the stars and gas surrounding them is, by far, the most common technique for accurately estimating the mass of BHs. In this method, the measured motion of the stars or gas is used to constrain dynamical models. Meanwhile, the light emitted by stars and gas is used to infer the gravitational potential, which can separate the contribution to the potential from a BH and the rest of the galaxy, thereby providing the mass of the BH. However, to constrain the mass of the BH, astronomers need to probe the region of the galaxy where the BH dominates the potential. This region is known as the sphere of influence and is defined by: rinf l = GMBH /σ 2 (1.1) where G is Newton’s gravitational constant, MBH is the mass of the black hole, and σ is the galaxy velocity dispersion at the effective radius. Early on, the spectroscopic capabilities of the Hubble Space Telescope (HST) provided the only means to probe the sphere of influence region in galaxies outside our local group (Kormendy & Ho 2013). This created a plateau in the number of BH mass measurements for nearly 15 years. In the mid 2000s, a new technology emerged capable of resolving the sphere of influence region from ground-based measurements. This technology, known as adaptive optics (AO), is used to improve photometric and spectroscopic observations by deforming a mirror to correct for distortions caused by atmospheric turbulence. These distortions are measured using either a natural guide star (a bright nearby star) or a laser guide star (an artificial star created by a laser on the telescope). AO paved the way to extend the types and number of galaxies with BH mass measurements. Another approach to detecting BHs relies on the radiation emitted from material falling onto a BH. Any material that falls into a BH is both super-heated and accelerated. Due to conservation of angular momentum, this accreted matter generally forms a disk-like structure around the BH. Within the disk, internal friction causes the angular momentum to be transported radially outward, allowing the material to fall farther inward, further increasing the temperature. This activity leads to the emission of X-rays. In many of these cases, the emitted material ejected at high speeds can also form relativistic jets, which can be detected with telescopes. As such, some of the most energetic events in the Universe have been attributed to the accretion of matter onto BHs. In particular, AGN and quasars are believed to be the accretion of matter onto SMBHs at the center of galaxies (Celotti et al. 7 1999). By modeling the dynamics of the material falling onto the BH, through processes such as reverberation mapping, the mass of the BH can be determined. The final approach to detecting BHs, through the emission of gravitational waves, has only recently emerged as observationally possible. Albert Einstein predicted the existence of gravitational waves in 1916 in his general theory of relativity. Einstein showed that massive, dense objects orbiting each other, such as neutron stars or black holes, would disrupt space-time in such a way that waves would radiate outward from the source at the speed of light. However, even he doubted his mathematical conclusions when it was pointed out that changing the coordinate system could isolate the unwanted singularities associated with gravitational waves. In 1976, observations confirmed the existence of a binary system composed of a pulsar in orbit around a neutron star (Taylor et al. 1976). A few years later, in 1989, it was confirmed that the orbit of this binary system was slowly shrinking due to the emission of energy in the form of gravitational waves (Taylor & Weisberg 1989). Nearly 100 years after Einstein’s initial prediction, the Laser Interferometer Gravitational-Wave Observatory (LIGO) detected the distortions in space-time for the first time, confirming the existence of gravitational waves (Abbott et al. 2016a). The properties of these waves provide us with information about the orbiting objects such as their mass, size, and orbital time scales. 1.2.2 Black Hole Scaling Relations The combination of these detection techniques have provided astronomers with a large sample of BH mass determinations. Interestingly, these measurements have revealed a striking correlation between the mass of the BH and a number of the host galaxies’ properties (e.g., Ferrarese & Ford 2005; Kormendy & Ho 2013; Georgiev et al. 2016; Graham 2016; Saglia et al. 2016). The most well known correlation is with the mass of the central BH and the stellar velocity dispersion of the galaxy, as shown in Figure 1.2 (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002; Saglia et al. 2016). Here, the data are taken from Saglia et al. (2016). These correlations are largely based on black hole mass estimates that were derived with dynamical models. Amazingly, these scaling relations span several orders of magnitude, from the biggest galaxies (M > 1012 M ) with black holes of ∼ 109 M to black holes of ∼ 105 M in dwarf galaxies (M ∼ 109 M ). However, very few BH masses have been measured in the galaxy mass regime below ∼ 109 M . It is the goal of this dissertation to measure the mass and understand the implications of central massive black holes in this low-mass regime using dynamical modeling techniques. 8 Figure 1.2: Correlation between the mass of a central black hole and the line-of-sight velocity dispersion. Here, I reproduced the data from Figure 13 of Saglia et al. (2016) with the 98 galaxies in their sample shown as black points. The red line indicates the best-fit from their parameters. 9 1.3 Dynamical Models Collisionless models are often a very good approximation of the stellar dynamics of galaxies because gravity is the dominant force. The most widely used tools to construct dynamical models are the Jean Anisotropic Modeling method (JAM) and the Schwarzschild dynamical modeling method. In these techniques, current methodologies require the luminosity and mass models of the galaxies to be expressed as a series of Gaussians (i.e., multi-Gaussian expansion (MGE); Cappellari 2002). Both of these methods are described briefly in Chapters 2 and 3, but I provide the general overview of each modeling method here. 1.3.1 Ingredients for Dynamical Modeling In order to create dynamical models that can accurately reproduce the measured motion of the stars or gas in a galaxy, one needs to determine the luminosity and mass density. Furthermore, to determine the mass of a central BH, kinematic measurements are needed in the region where the BH dominates the gravitational potential (within the sphere of influence). These kinematic measurements are generally determined from AO-assisted integral field unit (IFU) data mounted on 8-10m telescopes, which provide three-dimensional spectral data cubes where two dimensions represent a recreated image of the galaxy, and the third dimension represents the wavelength. To determine the luminosity and mass density, imaging data are needed, preferably, with better spatial resolution than the spectroscopy. Currently, the only source of such data is the HST. In the dynamical modeling techniques discussed in this dissertation, I use HST imaging to determine the spatial distribution of the stars and determine the threedimensional distribution of mass using the multi-Gaussian expansion (MGE) formalism (Cappellari 2002). To summarize, the MGE method consists of a series expansion of galaxy images using two-dimensional Gaussians, which, when summed together, reproduce the overall surface brightness profile of the galaxy. An example of a one-dimensional MGE fit is shown in Figure 1.3 where the black points show a r−2 profile, the black line shows the summed total of the Gaussians, and the red lines show individual Gaussians. This approach has significant advantages over the conventional technique of fitting ellipses of increasing semimajor axis to the galaxy image. For example, strong deviations from ellipses are hard to model with a small number of Fourier terms. Since the MGE method is not restricted to a certain number of Gaussians, it can easily capture multiple structural components of a galaxy without deviating in shape. Furthermore, deprojecting 10 Figure 1.3: Example of one-dimensional MGE fit to R−2 profile is shown as black diamonds. The black line shows the summed total of the nine Gaussian functions. The red lines show the individual Gaussian contribution to the fit. 11 Gaussians to form a three-dimensional gravitational potential for use in dynamical models is significantly easier than deprojecting ellipses. In this dissertation, the luminosity and mass density profiles are constructed using the MGE formalism. For every galaxy presented here, I fitted the galaxy image with Sérsic function profiles. The Sérsic profiles are then parameterized by MGEs. These MGEs are manipulated using different scale factors, discussed below, to construct the luminosity and mass density for use in the dynamical models. 1.3.2 JAM The JAM method is a generalization of the solutions to the collisionless Boltzmann equation, which describes the motion of stars in a gravitational potential and is defined as: 3 X ∂Φ ∂f ∂f (vi − )=0 (1.2) ∂xi ∂xi ∂vi i=1 Here, x represents the position, v is the velocity, Φ is the gravitational potential, and f is the distribution function, which is a function of six variables; and therefore, an infinite family of solutions exists. However, the collisionless Boltzmann equation can be simplified by assuming a cylindrical coordinate system, and axisymmetry as: vφ2 vR vφ ∂f ∂f ∂f ∂Φ ∂f ∂Φ ∂f vR + vz +( − ) − − = 0. (1.3) ∂R ∂z R ∂R ∂vR ∂z ∂vz R ∂vφ Multiplying this equation by vR and vz , and integrating over all velocities yields the two Jeans equations: 2 − νv 2 νvR φ 2) ∂(νvR ∂(νvR vz ) ∂Φ + = −ν R ∂R ∂z ∂R νvR vz ∂(νvz2 ) ∂(νvR vz ) ∂Φ + + = −ν R ∂z ∂R ∂z + (1.4) (1.5) where: Z νvk vj = vk vj f d3 v. (1.6) Here, ν represents the density of stars. The JAM method attempts to solve these equations by assuming a distribution function under the influence of a gravitational potential defined by the MGEs discussed in the previous section. The full implementation details of this method are described in Cappellari (2008). To summarize, this method makes two general assumptions: 1.) The mass-to-light ratio (M/L) is constant; and 2.) The velocity ellipsoid is aligned with cylindrical coordinates (R, z) and characterized by the anisotropy parameter βz = 1 − vz2 2 , vR where vz is the velocity parallel to the rotation axis and vR is the velocity in the radial direction in the plane of the galaxy. From here, this method is applied in a series of steps. First, a gravitational potential is constructed from a mass model of the galaxy. This potential can also contain the mass of a theoretical central BH. Next, 12 the potential is applied to find solutions to the Jeans equations (Jeans 1915). Finally, the intrinsic quantities are convolved with the point spread function (PSF) of the kinematic data to generate observable quantities that can be compared to kinematic measurements. The general caveat of this method is that the existence of a solution to the Jeans equations does not guarantee the existence of a corresponding physical distribution function, which can only be verified using a different modeling technique (Cappellari 2008). 1.3.3 Schwarzschild Dynamical Models The Schwarzschild dynamical modeling method is a three-integral dynamical modeling technique based on Schwarzschild’s numerical orbit superposition method (Schwarzschild 1979). Here, the three integrals are defined by the conservation laws under the influence of a gravitational potential, which are conservation of energy, conservation of angular momentum along the symmetry axis Lz , and I3 (the analytically unspecified integral of motion). In theory, this approach is much simpler than the JAM method. Starting from a predefined number of stars, Schwarzschid’s method samples all possible orbits in a gravitational potential. The orbits are then weighted and added together to reproduce the observed surface brightness and kinematic measurements of a galaxy. However, in practice, the computational resources needed to sample all possible orbits is expensive. Furthermore, the systematic effects, such as the stellar mass distribution within the black hole sphere of influence, have a strong impact on the results of the Schwarzschild models. Therefore, comparisons between the results of this method and the JAM method are of vital importance for black hole mass measurements as a way to interpret and possibly fix the problems associated with each method individually. For example, as discussed above, a solution to the Jeans equations does not guarantee the existence of a physical distribution function. On the other hand, the increased generality of the Schwarzschild method leaves it susceptible to numerical artifacts, which may not be physical. Part of the goal of this dissertation is to present a direct comparison of the two techniques as a way to disentangle possible advantages and disadvantages of each method. This implementaion details of the Schwarzschild models used in this dissertation are discussed in detail in van den Bosch et al. (2008). To briefly summarize, Schwarzschild’s method is also applied in a series of steps, with the first step, to create a gravitational potential, being identical to the JAM method. Next, the equations of motion are integrated numerically for a representative library of orbits. Finally, the weighted superposition of the orbits is determined for which the surface brightness and the full line of sight velocity distribution best fit the observed surface brightness and kinematics. In this way, the 13 Schwarzschild dyncamical modeling method not only provides the best fit parameters (e.g., BH mass, M/L, viewing angles, etc.), but also allows investigation of the distribution function through orbital mass weights. 1.4 Outline of the Dissertation This dissertation consists of two separate projects addressing galaxy formation and evolution of low mass galaxies. In particular, I present the search for supermassive black holes in three of the most massive Virgo cluster UCDs using dynamical models. These include a search for SMBHs using the JAM method in the UCDs, VUCD3, and M59cO (Chapter 2), and a comparison between the BH mass measurement in M59-UCD3 with the JAM and Schwarzschild methods (Chapter 3). In Chapter 2, I present the detection of SMBHs in the two Virgo UCDs, VUCD3, and M59cO. In this study, we use adaptive optics assisted data from the Gemini/NIFS instrument to derive radial velocity dispersion profiles for both objects. Mass models for the two UCDs are created using multiband (HST) imaging, including the modeling of mild color gradients seen in both objects. We then find a best-fit stellar M/L and BH mass by combining the kinematic data and the deprojected stellar mass profile using JAM. Assuming axisymmetric isotropic Jeans models, we detect BHs in both objects with masses +2.5 6 6 of 4.4+2.5 −3.0 × 10 M in VUCD3 and 5.8−2.8 × 10 M in M59cO (3σ uncertainties). The BH mass is degenerate with the anisotropy parameter, βz ; for the data to be consistent with no BH requires βz = 0.4 and βz = 0.6 for VUCD3 and M59cO, respectively. Comparing these values with nuclear star clusters shows that while it is possible that these UCDs are highly radially anisotropic, it seems unlikely. These detections constitute the second and third UCDs known to host supermassive BHs. They both have a high fraction of their total mass in their BH; ∼13% for VUCD3 and ∼18% for M59cO. They also have low best-fit stellar M/Ls, supporting the proposed scenario that most massive UCDs host high mass fraction BHs. The properties of the BHs and UCDs are consistent with both objects being the tidally stripped remnants of ∼109 M galaxies. In Chapter 3, I present an examination of the internal properties of the most massive UCD, M59-UCD3, by combining adaptive optics assisted near-IR integral field spectroscopy from Gemini/NIFS, and (HST) imaging. We use the multiband HST imaging to create a mass model that suggests and accounts for the presence of multiple stellar populations and structural components. We combine these mass models with kinematics measurements from Gemini/NIFS to find a best-fit stellar M/L and BH mass using JAM, axisymmetric 14 Schwarzschild models, and triaxial Schwarzschild models. The best fit parameters in the JAM and axisymmetric Schwarzschild models have black holes between 2.5 and 5.9 million solar masses. The triaxial Schwarzschild models point toward a similar BH mass, but show a minimum χ2 at a BH mass of ∼ 0. Models with a BH in all three techniques provide better fits to the central Vrms profiles, and thus we estimate the BH mass to be 6 4.2+2.1 −1.7 × 10 M (estimated 1σ uncertainties). We also present deep radio imaging of M59-UCD3 and two other UCDs in Virgo with dynamical BH mass measurements, and compare these to X-ray measurements to check for consistency with the fundamental plane of BH accretion. We detect faint radio emission in M59cO, but find only upper limits for M60-UCD1 and M59-UCD3 despite X-ray detections in both these sources. M59-UCD3 8 is the most massive UCD discovered to date (M? = 2.9+0.2 −0.3 × 10 M ), so the detection of a central massive black hole making up ∼2% of the total dynamical mass serves as an important test of the formation mechanisms of UCDs. The BH mass and nuclear light profile of M59-UCD3 suggests it is the tidally stripped remnant of a ∼109−10 M galaxy. CHAPTER 2 DETECTION OF SUPERMASSIVE BLACK HOLES IN TWO VIRGO ULTRACOMPACT DWARF GALAXIES This chapter originally appeared in the 2017 Astrophysical Journal (Ahn et al. 2017). In that format it is c 2017 by the American Astronomical Society. It is reprinted here with their permission. Sections, tables and figures have been renumbered and reformatted to match the overall flow of this thesis. The coauthors of the paper are Anil C. Seth, Mark den Brok, Jay Strader, Holger Baumgardt, Remco van den Bosch, Igor Chilingarian, Matthias Frank, Michael Hilker, Richard McDermid, Steffen Mieske, Aaron J. Romanowsky, Lee Spitler, Jean Brodie, Nadine Neumayer, and Jonelle L. Walsh. 2.1 Introduction Ultracompact dwarf galaxies (UCDs) are stellar systems discovered in the late 1990s through spectroscopic surveys of the Fornax cluster (Hilker et al. 1999; Drinkwater et al. 2000). With masses ranging from a few million to a hundred million solar masses and sizes . 100 pc, UCDs are among the densest stellar systems in the Universe. In the luminosity-size plane, UCDs occupy the region between globular clusters (GCs) and compact ellipticals (cEs) (e.g., Brodie et al. 2011; Misgeld & Hilker 2011; Norris et al. 2014; Janz et al. 2016). The smooth transition between these three classes of objects has led to significant debate as to how UCDs were formed. Explanations have ranged from UCDs being the most massive GCs (e.g., Fellhauer & Kroupa 2002, 2005; Mieske et al. 2002; Kissler-Patig et al. 2006) to UCDs being the tidally stripped nuclear remnants of dwarf galaxies (Bekki et al. 2001, 2003; Pfeffer & Baumgardt 2013; Strader et al. 2013; Forbes et al. 2014). Recently, analyses of the integrated dispersions of UCDs revealed an interesting property; the dynamical mass appears to be elevated ∼50% for almost all UCDs above 107 M when compared to the mass attributed to stars alone (e.g., Haşegan et al. 2005; Mieske et al. 2013). These dynamical mass estimates have been made combining structural information from Hubble Space Telescope (HST) imaging with ground-based, global velocity dispersion 16 measurements. These models assume that mass traces light, stars are on isotropic orbits, and are formed from a Kroupa-like initial mass function (IMF) (Haşegan et al. 2005; Mieske et al. 2008, 2013). Possible explanations for this unique phenomenon have included ongoing tidal stripping scenarios (Forbes et al. 2014; Janz et al. 2015) and central massive black holes (BHs) making up ∼10-15% of the total mass (Mieske et al. 2013). Alternatively, the elevated dynamical-to-stellar mass ratios can be explained by a change in the stellar IMF in these dense environments. For example, a bottom-heavy IMF would imply an overabundance of low-mass stars that contribute mass but very little light (Mieske & Kroupa 2008), and a top-heavy IMF would allow for an overabundance of stellar remnants contributing mass but virtually no light. The former case has been suggested in giant ellipticals (e.g., van Dokkum & Conroy 2010; Conroy & van Dokkum 2012), while Dabringhausen et al. (2012) argued that the relative abundance of X-ray binaries in UCDs favored a top-heavy IMF, although an increased X-ray luminosity in UCDs was not found in subsequent work (Phillipps et al. 2013; Pandya et al. 2016). In the context of the tidal stripping scenario, the elevated dynamical-to-stellar mass ratios could potentially be explained if UCDs still reside within progenitor dark matter haloes. However, to have a measurable effect on the kinematics of compact objects such as UCDs, the central density of the dark matter halo would need to be orders of magnitude higher than expected for dark matter halos of the stripped galaxies (Tollerud et al. 2011; Seth et al. 2014). In addition, the search for an extended dark matter halo in Fornax UCD3, based on its velocity dispersion profile, yielded a nondetection (Frank et al. 2011). In this paper, we follow up on the idea that the elevated values of the dynamical-to-stellar mass ratios, which we denote in this paper as Γ (≡ (M/L)dyn /(M/L)∗ ), can be explained by the presence of a supermassive BH (Mieske et al. 2013). This scenario was confirmed in one case, M60-UCD1, which hosts a BH that makes up 15% of the total dynamical mass of the system and a best-fit Γ of 0.7 ± 0.2 (Seth et al. 2014). As M60-UCD1 is one of the highest density UCDs, its low stellar mass-to-light ratio (M/L) suggests that a systematic variation of the IMF with density is not the cause for high M/L estimates found in most massive UCDs, and strengthens the case that these may be due to high mass fraction BHs. As part of our ongoing adaptive optics kinematics survey of UCDs, we investigate internal kinematics of two Virgo UCDs, VUCD3 and M59cO, with the goal of constraining the mass of putative central massive BHs. While we resolve the stellar kinematics of these two UCDs, they are fainter than M60-UCD1, and therefore have lower S/N data. This forces us to make more assumptions in our modeling. However, making reasonable assumptions, 17 we clearly detect BHs in both objects. Images of VUCD3, M59cO, and their host galaxies are shown in Figure 2.1. VUCD3 is located 14 kpc in projection from the center of M87 and has MV = −12.75 (Mieske et al. 2013). The metallicity of VUCD3 has been estimated to be between -0.28 and 0.35 in several studies (Evstigneeva et al. 2007; Firth et al. 2009; Francis et al. 2012), and it has an [α/Fe]∼0.5 (Francis et al. 2012). M59cO is located 10 kpc in projection from the center of M59 and has MV = −13.26 (Mieske et al. 2013). Its metallicity has been measured in several studies with [Z/H] between 0.0 and 0.2, with [α/Fe]∼0.2 (Chilingarian & Mamon 2008; Sandoval et al. 2015; Janz et al. 2016). We assume a distance of 16.5 Mpc for both objects. All magnitudes are reported in the AB magnitude system unless otherwise noted. All magnitudes and colors have been corrected for extinction; in VUCD3, we use AF 606W = 0.061 and AF 814W = 0.034, while for M59cO we used AF 475W = 0.107 and AF 850LP = 0.041 (Schlafly & Finkbeiner 2011). This paper is organized as follows: In Section 2.2, we discuss the data used for analysis and how the kinematics were modeled. In Section 2.3, we present our methods for determining the density profile of our UCDs. We present our dynamical modeling methods in Section 2.4. Our results for the best-fit BH are presented in Section 2.5, and we conclude in Section 2.6. VUCD3 M59cO M87 M59 Figure 2.1: The two galaxy-UCD systems discussed in this paper. The left panel shows the M87-VUCD3 system, and the right panel shows the M59-M59cO system. The main images show 2M ASSLGA images of both galaxies (Jarrett et al. 2003). VUCD3 and M59cO are the point-like images outlined in the yellow boxes. The inset images are zoom-in HST archival images of each UCD. The red line connecting the UCD to the host galaxy shows the projected distance, assuming each object is at a distance of 16.5 Mpc. 18 2.2 Observations and Kinematics In this section, we discuss the data and reduction techniques used for our analysis. Section 2.2.1 discusses the HST archival images, and Section 2.2.2 explains the reduction of our Gemini NIFS integral field spectroscopy. The derivation of kinematics is discussed in Section 2.2.3. 2.2.1 HST Data We obtained images of VUCD3 and M59cO from the Hubble Legacy Archive. VUCD3 was originally imaged as part of HST snapshot program 10137 (PI: Drinkwater). These data were taken using the High Resolution Channel (HRC) on the Advanced Camera for Survey (ACS), through the F606W and F814W filters. The ACS/HRC pixel scale is 0.02500 pixel−1 . The exposure times were 870s in F606W and 1050 s in F814W. M59cO was originally imaged as part of GO Cycle 11 program 9401 (PI: Côté). These data were taken using the ACS, Wide Field Camera (WFC), through the F475W and F850LP filters. The ACS/WFC pixel scale is 0.0500 pixel−1 . The exposure times were 750s in F475W and 1210s in F850LP. We synthesized a point spread function (PSF) for the HST images in each filter using the TinyTim software,1 MultiDrizzle,2 and following the procedure described in Evstigneeva et al. (2007). First, we generated the distorted PSFs with TinyTim, such that they represent the PSF in the raw images. We then place the PSFs in an empty copy of the raw HST flat fielded images at the location of the observed target. Finally, we pass these images through MultiDrizzle using the same parameters as were used for the data. This procedure produces model PSFs that are processed in the same way as the original HST data. The size of the PSFs was chosen to be 1000 × 1000 , which is much larger than the minimum size recommended by TinyTim in an attempt to minimize the effects of the ACS/HRC PSF halo problem (see §4.3.1 of Evstigneeva et al. (2008)). Although the ACS/WFC PSF does not suffer the same problems, we modeled the PSFs the same for consistency. The background (sky) level was initially subtracted from the HST images by MultiDrizzle. Instead of following the conventional way of estimating the sky level from empty portions of the image, we elected to add the MultiDrizzle sky level back in and model the background ourselves. This choice was motivated by the fact that these UCDs fall within the envelope of their host galaxy, and thus the background is not uniform across the image. To accomplish this, we first masked our original object and all foreground/background objects 1 http://www.stsci.edu/software/tinytim/ 2 http://stsdas.stsci.edu/multidrizzle/ 19 in the image. Next, we weighted the remaining good pixels (from the DQ extension of the image) by their corresponding error. Finally, we fitted a plane to the image to determine the sky level in each individual pixel; this was then subtracted from the image. 2.2.2 Gemini NIFS data Our spectroscopic data were obtained on 2015 May 2nd and 3rd and 2014 May 20th using the Gemini North telescope with the Near-Infrared Integral Field Spectrometer (NIFS) instrument (McGregor et al. 2003). The observations were taken using Altair laser guide star adaptive optics (Herriot et al. 2000; Boccas et al. 2006). The Gemini/NIFS instrument supplies infrared spectroscopy with a 300 field of view in 0.100 × 0.0400 pixels with spectral resolution λ δλ ∼ 5700 (σinst = 22 km s−1 ). The observations were taken in the K band, covering wavelengths from 2.0 to 2.4µm. The NIFS data were reduced following a similar procedure to the previous work done by Seth et al. (2010). The data were reduced using the Gemini version 1.13 IRAF package. Arc lamp and Ronchi mask images were used to determine the spatial and spectral geometry of the images. For M59cO, the sky images were subtracted from their closest neighboring onsource exposure. The images were then flat-fielded, had bad pixels removed, and split into long-slit slices. The spectra were then corrected for atmospheric absorption using an A0V telluric star taken on the same night at similar airmass with the NFTELLURIC procedure. Minor alterations were made to the Gemini pipeline to enable error propagation of the variance spectrum. This required the creation of our own IDL versions of NIFCUBE and NSCOMBINE programs. Each dithered cube was rebinned using our version of NIFCUBE to a 0.0500 × 0.0500 pixel scale from the original 0.100 × 0.0400 . These cubes were aligned by centroiding the nucleus, and combined using our own IDL program, which includes bad pixel rejection based on the nearest neighbor pixels to enhance its robustness. For VUCD3, sixteen 900s on-source exposures were taken with a wide range of image quality. We selected eight of the images with the best image quality and highest peak fluxes (two taken on May 2nd, and six taken on May 3rd, 2015) for a total on-source exposure time of two hours to create our final data cube. The data were dithered on chip in a diagonal pattern with separations of ∼100 . Due to the very compact nature of VUCD3, the surface brightness of the source in the sky region of the exposure that was subtracted from the neighboring exposure had very little signal (S/N of the sky portion of each individual image was always < 1). The final data cube for M59cO was created using twelve 900s on-source exposures (six taken on May 20th, 2014, two taken on May 2nd 2015, and four taken on May 3rd, 2015; 20 an additional five exposures were not used due to poor image quality) for a total on-source exposure time of three hours. We used an object-sky-object exposure sequence, and sky images were subtracted from both of their neighboring exposures. The object exposures were dithered to ensure the object did not fall on the same pixels in the two exposures with the same sky frame subtracted: this gives independent sky measurements for each exposure, improving the S/N relative to undithered exposures. The PSF for the kinematic data was derived by convolving a HST model image to match the continuum emission in the kinematic data cube. The HST model image was derived in K -band to best match the kinematic observations as follows. First, we generated 2D model images in each band using our best-fit Sérsic profiles described in Section 2.3. We then determined the color in each pixel. Next, we generated simple stellar population (SSP) models from Bruzual & Charlot (2003) using their PADOVA 1994 models at solar metallicities and assuming a Chabrier IMF. This is a reasonable assumption as both objects have near solar metallicities (see Section 2.1). Using our derived color and a color-color diagram from the SSPs, we determined the K -band luminosity for each pixel. The resulting image was convolved with a double Gaussian and Gauss+Moffat function and fitted to the NIFS image using the MPFIT2DFUN code3 (Markwardt 2009). In each case, we assessed which function provided the best fit to the PSF. For VUCD3, a Gauss+Moffat function was determined to provide the best fit; the Gaussian had a FWHM of 0.13800 and contained 29% of the light. The Moffat contained the remaining 71% of the light with a FWHM of 1.0800 that was parameterized by a series of Gaussians using the MGE-FIT-1D code developed by Cappellari (2002). For M59cO, a double Gaussian model was determined to provide the best fit with the inner component having a FWHM of 0.21600 and containing 53% of the light and the outer component having a FWHM of 0.93100 and containing 47% of the light. We note that in order to quantify the systematic effects on our choice of model PSF, we also report the results with the best-fit counterpart function (i.e., a double Gaussian for VUCD3, and a Gauss+Moffat function for M59cO). Furthermore, we verified the reliability of our PSF determination by comparing the results of fits to the HST model image to those from Lucy-Richardson deconvolved images. 2.2.3 Kinematic Derivation The kinematics were measured from the NIFS data by fitting the wavelength region between 2.29 µm and 2.37 µm, which contains the CO bandheads. Due to the low S/N 3 http://www.physics.wisc.edu/ craigm/idl/fitting.html 21 of the data (central pixel median S/N = 10 per 2.13 Å pixel for VUCD3 and S/N = 11 per 2.13 Å pixel for M59cO), it was not possible to make kinematic maps as were made for M60-UCD1 (Seth et al. 2014). We therefore constructed a radial dispersion profile for dynamical modeling. To create our dispersion profile, the data were binned radially such that the median S/N was ≈ 25 in each bin. The line spread function (LSF) was determined in each bin by combining sky exposures in the same dither pattern as the science exposures; λ we convolved the high-resolution Wallace & Hinkle (1996) stellar templates ( ∆λ = 45000) by the LSF in each radial bin before fitting. We fitted the radial velocity and dispersion to the data using the penalized pixel fitting algorithm pPXF (Cappellari & Emsellem 2004); a fourth order additive polynomial was also included in the fit. The integrated spectra and their corresponding fits are shown in Figure 2.2. The one sigma random uncertainties on the determined kinematics were estimated using Monte Carlo simulations. Gaussian random noise was added to each spectral pixel, and then we fitted kinematics and took the standard devation as the uncertainty. For VUCD3, a portion of the fit was left out due to bad sky subtraction. We found the integrated (r < 0.600 ) barycentric-correct radial velocity to be 707.3 ± 1.4 km s−1 , and the integrated dispersion to be 39.7 ± 1.2 km s−1 . This dispersion value is in agreement with the measurement in Evstigneeva et al. (2007) of 41.2±1.5 km s−1 using Keck/ESI data and a 1.500 aperture. For M59cO, we found the integrated (r < 1.100 ) velocity to be 721.2 ± 0.5 km s−1 , and the integrated dispersion to be 31.3 ± 0.5 km s−1 . The dispersion of M59cO is significantly lower than the measurement in Chilingarian & Mamon (2008) of 48±5 km s−1 ; however, this study was based on relatively low resolution SDSS spectra. Our values are consistent with the higher resolution Keck/DEIMOS observations presented in Norris et al. (2014) that find σ = 29.0 ± 2.5 km s−1 . Our binned, resolved values for the radial kinematics are given in Table 2.1. For each bin, we give the light weighted radius, S/N in each bin, line-of-sight velocity, dispersion, the corresponding uncertainty for both the velocity and dispersion, as well as the rotational speed and angle (see below). We also tested our velocity and dispersion values by fitting to the Phoenix stellar templates (Husser et al. 2013). Our results are consistent within one sigma for all dispersion values, while we see velocity discrepancies of up to 6 km s−1 , suggesting the velocity measurement uncertainties may be underestimated due to template mismatch or sky subtraction issues (especially in VUCD3 due to low S/N ). To calculate the rotational speed, we first split the integrated velocity bin in half and rotated the line separating the two halves through 360◦ in increments of 5◦ , fitting the 22 Figure 2.2: Integrated spectra of VUCD3 (left) and M59cO (right) are shown here in black. In both spectra, the red lines indicate the best kinematic fit and the residuals are shown in green. For visibility, the zero points of the residuals are given as the green lines at 2510 counts and 4406 counts for VUCD3 and M59cO, respectively. A portion of the VUCD3 fit had to be masked due to bad sky subtraction. The integrated dispersion was determined to be σ = 39.7 ± 1.2 km s−1 with a median S/N = 42 per pixel and σ = 31.3 ± 0.5 km s−1 with a median S/N = 69 per pixel for VUCD3 and M59cO, respectively. velocity on each side. Next, we fitted a sinusoidal curve to the difference between the two halves as a function of angle. Using the angle where the sinusoidal curve was either maximum or minimum, we created a line to split the radial bins in half and quote the difference in Table 2.1. We note that the rotational speed quoted is of order unity of the true rotation value (the velocity difference is a factor of 4/πvrot assuming smooth azimuthal variation). Neither object is rotation dominated, but VUCD3 shows substantial rotation oriented roughly along the major axis with v/σ ∼ 0.5; this is similar to M60-UCD1 (Seth et al. 2014). We note that only the dispersion values for the full annuli are used for the Vrms profile for fits to our dynamical models. We note two important characteristics of the M59cO kinematics. First, from the low rotational velocity combined with the nearly circular Sérsic profile fits, discussed below, M59cO appears to be nearly face-on or simply nonrotating. Second, there appears to be an upturn in the velocity dispersion for the last radial bin. This upturn could be due to several possibilities. First, it may be due to a problem in the characterization of the PSF on large scales resulting in scattered light from the center of the UCD to large radii; however, tests Radius [00 ] 0.033 0.069 0.118 0.217 0.437 0.025 0.065 0.105 0.145 0.191 0.245 0.317 0.475 S/N 38.27 41.66 40.10 25.96 15.72 33.07 45.84 47.87 49.65 47.92 42.31 43.49 33.42 v [km s−1 ] 709.6 713.7 715.6 704.0 706.7 720.8 719.5 719.4 718.5 720.1 720.3 724.5 725.2 verr [km s−1 ] 3.0 3.0 2.8 1.7 1.8 2.0 1.7 1.7 1.6 1.4 1.3 1.1 0.9 σ [km s−1 ] 52.9 51.2 49.0 40.3 33.1 40.2 39.9 37.6 34.9 33.6 31.8 28.4 29.6 σerr [km s−1 ] 2.5 2.2 2.1 1.6 2.1 1.6 1.4 1.3 1.2 1.1 1.0 1.0 0.7 Rotation1 [km s−1 ] 4.7 7.7 18.7 22.5 4.3 5.3 7.5 8.3 7.0 3.7 4.7 5.9 3.8 PA2 [◦ ] 77 77 82 112 127 82 82 47 52 22 7 -23 -28 Rotation is the maximum difference (amplitude) between the two halves of the radial bin split by a line at the PA. This value is on the order of unity of the true rotational velocity (amplitude is a factor 4/π vrot ). 2 The PA orientation is N=0◦ and E=90◦ . 1 M59cO Object VUCD3 Table 2.1: Resolved Kinematics 23 24 with alternative PSF models with more power in the wings did not create the observed upturn in our dynamical models. Alternatively, it may be due to poor sky subtraction at large radii, tidal stripping at the edges of the UCD, or background light contamination from the host galaxy; assuming one of these is the cause, we exclude this data point from our dynamical modeling. 2.3 Creating a Mass Model To create dynamical models predicting the kinematics of our UCDs, we first needed to create a model for the luminosity and mass distribution in each object. While typically, the mass is assumed to trace the light (e.g., Seth et al. 2014), in both UCDs considered here previous works have detected color gradients, suggesting a nonconstant M/L (Chilingarian & Mamon 2008; Evstigneeva et al. 2008). Fortunately, we have two filter data available for both UCDs, and we make use of the surface brightness profile fits in both filters to estimate the luminosity and mass profiles of the UCDs. We consider the uncertainties from different combinations of luminosity and mass models in our best-fit dynamical models in Section 2.4, and find that they do not create a significant uncertainty in our BH mass determinations. Neither source is well fit using a single Sérsic profile, and both appear to have two components (Evstigneeva et al. 2007; Chilingarian & Mamon 2008). We therefore determine the surface brightness profile by fitting the data in each filter to a PSF-convolved, two component Sérsic profile using the two-dimensional fitting algorithm, GALFIT (Peng et al. 2002). The parameters in our fits are shown in Table 2.2 and include, for each Sérsic profile: the total magnitude (mtot ), effective radius (Re ), Sérsic exponent (n), position angle (P A), and axis ratio (q). The fitting was done in two ways. First, we allowed all of the above free parameters to vary in both filters; these fits are henceforth referred to as the “free” fits. Next, we fitted the data again fixing the shape parameters of one filter to the best-fit model from the other filter; specifically, we fixed the effective radius, Sérsic exponent, position angle and axis ratio, allowing only the total magnitude to vary. For example, in VUCD3, our fixed fit in F814W was done by fixing the shape parameters to the best-fit free model in F606W. By using the same shape parameters, these “fixed” fits provide a well-defined color for the inner and outer Sérsic profiles. These Sérsic profile fits, shown in Figure 2.3 for the filters used to create the luminosity and mass models, were performed on a 500 × 500 image centered on each UCD with a 100 × 100 pixel convolution box. We note that M59cO also has data in the F850LP filter available, which we use to model the color gradients as discussed below. However, due to the lack of χ2 14.16 6.342 1.065 0.974 mtot 18.99 18.55 19.30 18.13 Inner Sérsic Re [00 ] n 0.08 3.51 0.66 0.07 3.25 0.62 0.16 1.06 0.97 0.15 1.02 0.99 PA [◦ ] 19.0 18.0 -65.2 34.1 mtot 18.83 17.96 18.27 16.68 Outer Sérsic Re [00 ] n 0.61 1.28 0.91 0.62 1.74 0.89 0.64 1.09 0.98 0.61 1.21 0.98 PA [◦ ] 18.4 20.7 88.4 17.7 min,tot Fixed 19.13 18.40 19.45 17.94 mout,tot Fixed 18.66 18.11 18.21 16.74 Notes. The last two columns show the total magnitude when all shape parameters of the Sérsic profiles are held fixed to the other filter. 1 Object VUCD3 (F606W) VUCD3 (F814W) M59cO (F475W) M59cO (F850LP) Table 2.2: Best-fit Sérsic Parameters of VUCD3 and M59cO 25 26 Figure 2.3: Surface brightness profiles of VUCD3 (left) and M59cO (right) in HST filters used for dynamical modeling. Black stars are data, cyan lines are convolved models, red lines are the double Sérsic reconstructed profile, and green and blue lines are the individual Sérsic components. The residuals between the data and convolved models are shown in the bottom panel of each figure. a red cutoff in the filter, the PSF is difficult to characterize 4 ; therefore, we chose to use the Sérsic fits to the F475W filter as the basis for our luminosity and mass models. The outputs for the best fitting Sérsic profiles in each filter are shown in Table 2.2. For VUCD3, the total luminosity and effective radius calculated from the double Sérsic profile was found to be LF 814W = 17.8 × 106 L and Re = 18 pc or 0.22500 , with Sérsic indices of 3.25 for the inner component and 1.74 for the outer component. These indices are similar to what was found for M60-UCD1 (Strader et al. 2013). For M59cO, the total luminosity and effective radius was found to be LF 475W = 20.3 × 106 L and Re = 32 pc or 0.400 , with Sérsic indices of 1.06 and 1.21 for the inner and outer components, respectively. These values are comparable to the n ∼ 1 used to fit the system in previous work (Chilingarian & Mamon 2008). The best-fit Sérsic profiles were then parameterized by a series of Gaussians or MGEs for use in our dynamical models (Emsellem et al. 1994; Cappellari 2002). For a uniform stellar population, the luminosity profile in any band can be used to obtain an accurate mass model. However, if the stellar populations vary with position, we need to take this into account in our dynamical modeling. We tested for stellar population variations by creating color profiles as shown in Figure 2.4. It is clear that both objects show 4 http://www.stsci.edu/hst/HSToverview/documents/calworkshop/workshop2002/CW2002Papers/gilliland.pdf 27 Figure 2.4: The color profiles of VUCD3 (left) and M59cO (right) are shown here in black diamonds based on HST data. The solid lines indicate the two-component Sérsic model fits that have been convolved with the PSF, while dashed lines indicate model fits that are unconvolved. The colors represent whether the parameters of the Sérsic fits were left independent (black) or fixed to the other band (red and blue). Blue lines indicate that the shape parameters of the Sérsic profile in the bluest filter were fixed to that of the redder filter, while red lines are vice versa. The fixed unconvolved models provide a well-defined color for the inner and outer Sérsic profiles. For VUCD3, the inner color is 0.59 mag and the outer color is 0.72 mag. For M59cO, the inner and outer colors are 1.36 mag and 1.53 mag, respectively. The effects of convolution with the PSF make the color differences more extreme. a trend toward redder colors as a function of radius, confirming the color gradients shown in previous works (Chilingarian & Mamon 2008; Evstigneeva et al. 2008). Therefore, we could not assume a mass-follows-light model as was done with M60-UCD1 (Seth et al. 2014). To incorporate the stellar population variations into our dynamical models, we needed two ingredients: (1) a mass profile of the object needed for dynamical modeling, and (2) a K band luminosity profile to enable comparison of our dynamical models with our observed data. The unconvolved fixed models (dashed blue and red lines in Figure 2.4) provide a welldefined color for the inner and outer components of the Sérsic fits since the shape parameters are held fixed. For VUCD3, the inner component color is F 606W − F 814W = 0.59 mag, while the outer component color is F 606W − F 814W = 0.72 mag. In M59cO, we found the color to be F 475W − F 850LP = 1.36 and 1.53 mag for the inner and outer components, 28 respectively. These colors were then used to find the mass-to-light ratio, assuming solar metallicity, using the Bruzual & Charlot (2003) Padova 1994 SSP models. To evaluate the errors on our M/L, we assumed an error of ±0.02 mags in our color determinations. For VUCD3, we found the inner component M/LF 814W to be 1.4±0.4, with a corresponding age of 9.6 Gyrs, and the outer component to be 2.7±0.4, with a corresponding age of 11 Gyrs. For M59cO, we found 2.8+0.3 −0.2 and 5.5±0.5 for the inner and outer M/LF 475W , with corresponding ages of 5.5 and 11.5 Gyrs, respectively. To determine a mass density profile in M pc−2 to be used in the dynamical models, we multiplied the luminosities of each MGE subcomponent by these M/Ls. These M/L values can also be used to estimate total masses of the inner and outer Sérsic components. For VUCD3, we found the inner component to contain 11 ± 3 × 106 M and the outer component 27 ± 4 × 106 M . For M59cO, we found 6 14+2 −1 and 84 ± 8 × 10 M for the inner and outer components, respectively. We note that we also computed a mass profile for the free Sérsic fits where we determined the M/L from the color at the FWHM of each Gaussian component in the MGE light profile (discussed in Section 2.4). Our best fit mass model MGE for each UCD is given in the 1st column of Table 2.3. The color profiles and SSP models were also used to calculate the K-band luminosity MGEs that are used to compare our dynamical models to the kinematic data. Since the unconvolved fixed models provide an accurate determination of the color profile of each UCD, we used these colors, described above, to create a K-band luminosity profile for the dynamical models. In both cases, we use the BC03 models to infer the colors between our best-fit model and K-band. For VUCD3, we find that the inner component has F 814W − K = 2.14 and outer component has F 814W − K = 2.26. This leads to a scale factor in luminosity surface density of 0.44 and 0.40 for the inner and outer components, respectively. For M59cO, we found the inner component F 475W −K = 3.35 with a scale factor of 0.24 and the outer component F 475W −K = 3.58 with a scale factor of 0.20. These scale factors were multiplied by the inner and outer component luminosity profiles to make K-band MGEs for use in the dynamical models. Our best fit K-band luminosity model MGE for each UCD is given in the second column of Table 2.3. Finally, the color profiles and SSP models were used to calculate the total stellar population M/L. This was accomplished by first calculating the flux within the central 2.500 from model images of the inner and outer Sérsic profiles. Next, we used the M/L calculated from the color profiles to find a flux weighted total M/L. For VUCD3, we found M/LF 814W,∗ = 2.1 ± 0.6, which, assuming V − I = 1.27 based on observations (Evstigneeva 29 Table 2.3: Multi-Gaussian Expansions (MGEs) Used as the Default Mass and Luminosity Models in Dynamical Modeling. Object VUCD3 M59cO 1 M ass (M pc−2 )1 1687276. 1435008. 1145834. 798986.4 489960.9 265194.5 123165.2 48162.50 15761.48 4225.199 948.0368 173.9253 25.42110 2.945087 0.192683 938.0369 1724.193 2272.561 2018.410 1101.361 348.8148 57.96370 3.441089 3820.00 8870.65 13691.4 12547.7 6057.40 1376.19 95.0511 2082.22 4051.69 5736.21 5457.66 3284.28 1137.78 206.300 12.7642 IK (L pc−2 )2 537235. 456912. 364838. 254400. 156005. 84439.0 39216.4 15335.1 5018.52 1345.32 301.859 55.3785 8.09415 0.93773 0.06135 135.776 249.400 328.719 291.958 159.309 50.4551 8.38429 0.49774 328.197 762.126 1176.31 1078.04 520.425 118.236 8.16636 75.4811 146.874 207.938 197.841 119.056 41.2449 7.47872 0.46271 σ(00 ) .0001 .0003 .0008 0.002 0.005 0.009 0.019 0.038 0.072 0.131 0.228 0.387 0.640 1.046 1.836 0.012 0.044 0.121 0.262 0.484 0.795 1.210 1.808 0.005 0.019 0.047 0.092 0.151 0.223 0.315 0.013 0.047 0.123 0.262 0.473 0.761 1.138 1.670 q 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.91 0.91 0.91 0.91 0.91 0.91 0.91 0.91 0.99 0.99 0.99 0.99 0.99 0.99 0.99 0.98 0.98 0.98 0.98 0.98 0.98 0.98 0.98 PA (◦ ) 19.04 19.04 19.04 19.04 19.04 19.04 19.04 19.04 19.04 19.04 19.04 19.04 19.04 19.04 19.04 18.43 18.43 18.43 18.43 18.43 18.43 18.43 18.43 34.06 34.06 34.06 34.06 34.06 34.06 34.06 17.69 17.69 17.69 17.69 17.69 17.69 17.69 17.69 The M/L in F814W (VUCD3) and F475W (M59cO) were used to determine the mass profiles for dynamical modeling. These methods are described in detail in Section 2.3 2 Creation of MGEs required an assumption on the absolute magnitude of the sun. These values were assumed to be 4.53 mag in F 814W , 5.11 mag in F 475W , and 3.28 mag in K taken from http://www.ucolick.org/∼ cnaw/sun.html. See Section 2.3 for a discussion on how the K-band luminosity was determined. 30 et al. 2007) corresponds to M/LV,∗ = 5.2±1.5. We found M/LF 475W,∗ = 4.8+0.6 −0.5 for M59cO. For the overall object we estimate a g − V ∼ 0.47, yielding a M/LV,∗ = 4.1+0.5 −0.4 . Both values of M/LV are consistent with the 13 Gyr population estimates in Mieske et al. (2013). 2.4 Dynamical Modeling In this section we describe the technical details of the dynamical modeling, while the results of the modeling are presented in the next section. We fit the radial dispersion profiles of each UCD to dynamical models using the Jeans Anisotropic Models (JAM) method with the corresponding code discussed in detail in Cappellari (2008). To briefly summarize, the dynamical models are made in a series of steps making two general assumptions: (1) The velocity ellipsoid is aligned with the cylindrical coordinate system (R, z, φ); (2) The anisotropy is constant. Here, the anisotropy is defined as βz = 1 − (σz /σR )2 where σz is the velocity dispersion parallel to the rotation axis and σr is the velocity dispersion in the radial direction in the plane of the galaxy. The first step in the dynamical modeling process is to construct a three-dimensional mass model by deprojecting the two-dimensional mass model MGEs discussed in the previous section. In the self-consistent case, the luminosity and mass profile are the same. However, in our case, we used the mass profile to construct the potential, and we used the light profile to calculate the observable properties of the model, both described below. The choice to parameterize the light profile with MGEs is motivated by the ease of deprojecting Gaussians and the accuracy in reproducing the surface brightness profiles (Emsellem et al. 1994; Cappellari 2002). The second step in the dynamical modeling process is to construct a gravitational potential using our mass model. This potential also contains a Gaussian to represent a supermassive BH with the axis ratio, q = 1, and width, σ . rmin /3, where rmin is the smallest distance from the BH that needs to be accurately modeled. Although a supermassive BH can be modelled by adding a Keplerian potential, it is much simpler to model the BH as this small Gaussian (Emsellem et al. 1994). Next, the MGE formalism is applied to the solution of the axisymmetric anisotropic Jeans equations (see Section 3.1.1 of Cappellari 2008). Finally, the intrinsic quantities are integrated along the line-of-sight (LOS) and convolved with the PSF from the kinematic data to generate observables that can be compared with the radially binned dispersion profiles. Supermassive BH masses are frequently measured with dynamical models that allow for fully general distribution functions (e.g., Schwarzschild), which is important to include because of the BH mass-anisotropy degeneracy in explaining central dispersion peaks in galaxies. Since 31 plunging radial orbits have an average radius that is far from the center of the galaxy, these orbits can raise the central dispersion without significantly enhancing the central mass density. Similarly, a supermassive BH also raises the dispersion near the center of the galaxy. Other dynamical modeling techniques break this degeneracy by fitting the full orbital distribution without assumptions about the anisotropy. However, given the quality of our kinematic data, a more sophisticated dynamical modeling technique is not feasible; we further discuss the assumptions and limitations of our modeling at the beginning of Section 2.5.1. For our dynamical models, we created a grid of the four free parameters: Γ, BH mass, inclination angle, and anisotropy. For VUCD3, our initial run consisted of the following: • 40 values of Γ ranging from 0.05 to 2.0 in increments of 0.05. Note that the best-fit dynamical M/LF 814W = 2.1Γ. We use this translation when reporting the dynamical M/Ls for VUCD3. • 16 values for the BH mass from 0.0 to 7.0 × 106 M in increments of 5 × 105 M plus one point at 1 × 105 M • 11 values for an anisotropy parameter in increments of 0.1 from −0.2 to 0.8 • 4 values for the inclination angle from 60◦ to 90◦ in 10◦ increments M59cO was fitted using the following: • 35 values of Γ from 0.1 to 3.5 in increments of 0.1. Note that the best-fit dynamical M/LF 475W = 4.8Γ. We use this translation when reporting the dynamical M/Ls for M59cO. • 19 values for the BH mass from 0.0 to 8.5 × 106 M in increments of 5 × 105 M plus one point at 1 × 105 M • 11 values for an anisotropy parameter from −0.2 to 0.8 in 0.1 increments • 9 values for the inclination angle from 14.5◦ (the lowest possible) to 90◦ in 10◦ increments The grids for Γ, BH mass and anisotropy values are explained in further detail below. To determine the best-fit BH mass, we assumed isotropy (motivation explained in Section 2.5) and marginalized over Γ and the inclination angle. Next, we computed the cumulative likelihood, shown in Figure 2.5. 32 Figure 2.5: Cumulative likelihood, assuming isotropy, of the BH mass in VUCD3 (left) and M59cO (right). In both figures, the black solid line represents the best-fit model. The red, blue, and cyan lines indicate PSF, mass, and luminosity model variations, repectively. The grey vertical lines indicate the best-fit, one σ, and three σ BH mass estimates. See Section 2.4 for further explanation of the individual red, blue, and cyan lines. Here, the different linestyles and colors represent different models and variations in the kinematic PSF. Unless explicitly stated, all of our dynamical models make use of the K-band luminosity MGEs; we note that the Γ values are scalings of our mass model, and the luminosity model is used only to calculate the model dispersion values. Our default dynamical models (shown in black) are as follows: • For VUCD3 the default mass model was obtained by fixing the best-fit double Sérsic model from the F814W data to the F606W data allowing only the Sérsic amplitudes to vary. The best-fit PSF was a Gauss+Moffat profile. • For M59cO the default mass model was obtained by fixing the best-fit double Sérsic model from the F475W data to the F850LP data allowing only the amplitude to vary. The best-fit PSF was a double Gaussian profile. To explore the systematic errors created by our choices of mass modeling and the fitting of the kinematic (NIFS) PSF, we also ran JAM models varying the mass model and PSF. We used our default mass model and varied only the PSF (shown in red) as follows: 33 • Solid: the best-fit PSF from the function that did not best match the continuum (i.e., a double Gaussian for VUCD3, and a Gauss+Moffat profile for M59cO). • Dotted: the PSF created using the HST model image in the reddest filter available for convolution. We also ran three separate JAM models with various mass models (shown in blue) as follows: • Solid: mass model using the best-fit free double Sérsic profile with the mass profile determined from the color at the FWHM of the individual MGEs. • Dashed: model using the best-fit free double Sérsic profile assuming mass follows light. • Dotted: model where the shape parameters of the double Sérsic profile were fixed assuming mass follows light. Finally, we tested the effects of our choice of the luminosity model by running one dynamical model with the default mass model and PSF, but using the luminosity model from the original filter (F814W for VUCD3 and F475W for M59cO; shown in cyan). The default model was chosen based on the accuracy of reproducing the surface brightness profiles, as well as the ease and accuracy of determining the luminosity and mass profiles. The systematic effects of our model and PSF variations were taken into account when reporting the uncertainties on our final results (see Section 2.5). We also ran a finer grid of models for our default isotropic models to better sample and obtain a best-fit value for the cumulative distributions (Figure 2.5) and predicted Vrms profiles (Figure 2.6). This smaller grid sampled the BH mass in 18 linear steps of 100,000 M ranging from: 3.4 million to 5 million M for VUCD3, and 5.2 million to 6.8 million M for M59cO. For comparison with the dispersion profile, we also included a zero mass BH for both objects. For VUCD3, Γ was modeled in 40 linear steps of 0.05, ranging from 0.05 to 2.0. For M59cO, we sampled the Γ in 35 linear steps of 0.1 ranging from 0.1 to 3.5. This final grid was modeled at the best-fit inclination angle from the first grid, which was 60◦ for VUCD3, and 14.5◦ for M59cO. However, the inclination angle has a negligible effect on the best-fit BH mass and M/L. The best-fit model from this grid was used to determine the final results (see Section 2.5.1). 2.5 Results In this section we report the results of our dynamical modeling, including the best-fit BH mass, and stellar M/L, assuming isotropy. We also report the impact of including anisotropic orbits in our dynamical modeling. 34 Figure 2.6: Dispersion profiles of VUCD3 (left) and M59cO (right) where black points are the measured velocity dispersions. The blue and red lines represent the best-fit isotropic models to the dispersion profile with and without a BH, respectively. The grey line represents the best-fit dynamical model without a BH, but including anisotropy. The grey point in the M59cO dispersion profile was not fitted. 2.5.1 Isotropic Model Results We start by considering the best-fit results for isotropic Jeans models. We are forced to adopt simple dynamical models given the quality of our data sets, with just a small number of radially integrated kinematic measurements. The assumption of isotropy is a resonable one; in M60-UCD1, for which higher fidelity data were available, the results from isotropic Jeans models were fully consistent with the more sophisticated Schwarzschild model results (Seth et al. 2014), but with somewhat smaller error bars due to the lack of orbital freedom in the Jeans models. Nearby nuclei, including the Milky Way, have also been found to be nearly isotropic (Verolme et al. 2002; Cappellari et al. 2009; Schödel et al. 2009; Hartmann et al. 2011; Nguyen et al. 2016) and their transformation into UCDs by tidal stripping is not expected to affect the mass distribution near the center of the galaxy (e.g., Pfeffer & Baumgardt 2013). We also note that other works have also shown JAM models consistent with Schwarzschild model and maser BH mass estimates (Cappellari et al. 2009, 2010; Drehmer et al. 2015). Given all of these factors, we present our results assuming isotropic Jeans models, and then consider the effects of anisotropy in the following section. +0.9 6 6 For the BH mass we found the best fit to be 4.4+0.6 −0.7 × 10 M and 5.8−0.5 × 10 M for 35 VUCD3 and M59cO, respectively. We found the best-fit M/LF 814W to be 1.8 ± 0.3 for VUCD3 and M/LF 475W = 1.6+0.3 −0.4 for M59cO. Here, the uncertainties are quoted as the one-sigma deviations calculated from the cumulative likelihood. Due to the lack of orbital freedom in the JAM models, we also quote the three-sigma deviations for both objects, which also encompass the systematic effects of the model/PSF variations. For VUCD3, we 6 found the best-fit BH mass and M/LF 814W to be 4.4+2.5 −3.0 ×10 M and 1.8±1.2, respectively. 6 For M59cO, we found 5.8+2.5 −2.8 × 10 M for the BH mass and 1.6+1.2 −1.1 for the M/LF 475W . Using the color information from Section 2.3, we found M/LV = 3.0 and 1.4 for VUCD3 and M59cO, respectively. Figure 2.6 shows the comparison of our kinematic data (black points) with the bestfitting dynamical model, using the values stated above for the mass of the BH and M/L. The red line represents the best-fit dynamical model without a BH, and the blue line represents the best-fit dynamical model with a BH. The grey line indicates the best-fit dynamical model without a BH, but including anisotropy (discussed in Section 2.5.2). Changing the mass of the BH affects the overall shape of the dispersion profile, while the M/L merely scales the model dispersion vertically. In both objects, it is clear that when isotropy is assumed a central massive BH better reproduces the kinematic dispersion profile. From Figure 2.7, we see that adding a central massive BH to the dynamical modeling has the effect of reducing the best-fit M/L (as well as increasing the anisotropy as discussed below). Therefore, we determined the total dynamical mass and Γ, assuming isotropy, as a check to our original hypothesis; the addition of a central massive BH reduces Γ to values comparable to globular clusters and compact elliptical galaxies. The total dynamical mass was calculated by multiplying the dynamical M/L with the total luminosity. For VUCD3, we found that with a BH mass of 4.4 × 106 M , Γ with three-sigma error bars was 0.8 ± 0.6 resulting in a total dynamical mass of 32 ± 21 × 106 M . For comparison, without a BH component Γ with three-sigma error bars was 1.7 ± 0.2 resulting in a total dynamical mass of 66 ± 8 × 106 M , which is consistent with previous results (Evstigneeva et al. 2007; Mieske et al. 2013). For M59cO, we found that with a BH mass of 5.8 × 106 M , Γ with 6 three-sigma error bars was 0.3 ± 0.2 resulting in a total dynamical mass of 32+24 −22 × 10 M . Without a BH component Γ with three-sigma error bars was 0.9 ± 0.1 resulting in a total 6 dynamical mass of 83+5 −12 × 10 M . These results for Γ and the total dynamical mass without a BH are inconsistent with previous works; it is lower than those based on a low resolution dispersion determination (Chilingarian & Mamon 2008; Mieske et al. 2013), and higher than the measurement by Forbes et al. (2014) based on a lower integrated dispersion 36 Figure 2.7: Contour plots showing the degeneracy between βz , MBH , and Γ, where VUCD3 is shown in the top panel figures and M59cO in the bottom. Data in the left panels have been marginalized over Γ (≡ (M/L)dyn /(M/L)∗ ), while data in the right panels are marginalized over βz . The blue points represent the extent of our grid over these parameters, and the green point represents the best fit determined for all free parameters. The black, blue and red contours represent the one, two, and three σ confidence levels, respectively, corresponding to ∆χ2 values of 2.3, 6.2, and 11.8 (assuming two degrees of freedom). The green, orange, and yellow lines correspond to the best fit Γ and βz assuming the BH mass makes up 1%, 5%, and 10% of the total dynamical mass, respectively. We note that the M/L∗ used corresponds to a M/LV of 5.2 for VUCD3 and 4.1 for M59cO. 37 measurement. Taking the ratio of the best-fit BH mass and the total dynamical mass (including both stars and BH), we found a central massive BH making up ∼13% and ∼18% of the total mass for VUCD3 and M59cO, respectively. These large mass fractions suggest a large black hole sphere of influence, which quanitifies the ability to detect a BH given a set of observations. Using the conventional definition of the black hole sphere of influence (rinf l = GM/σ 2 ), we find rinf l = 0.1500 for VUCD3 and rinf l = 0.3100 for M59cO (assuming the integrated dispersion values). This large sphere of influence is the reason for the large uncertainty in our stellar masses. The BH mass fractions are comparable to the mass fraction found in M60-UCD1 (Seth et al. 2014), and consistent with the estimates made by Mieske et al. (2013). They are also similar to the mass fractions of BHs within nuclear star clusters in the Milky Way, M32, and NGC 4395 (Graham & Spitler 2009; den Brok et al. 2015). Furthermore, these BH mass fractions reduce Γ to values comparable to those in many globular clusters and compact ellipticals (Strader et al. 2011; Forbes et al. 2014). Figure 2.8 illustrates this effect. Here the grey points represent globular clusters and UCDs. The stars represent objects which our group will analyze and test for the presence of central massive BHs. The blue stars are VUCD3 and M59cO while the red star shows M60-UCD1. These colored stars show Γ after accounting for a central massive BH. The colored arrows illustrate the effect that the central massive BH has on Γ and the dynamical estimate of the stellar mass. Figure 2.8 also illustrates the fact that all three UCDs with central massive BHs have stellar components with lower than expected dynamical masses (i.e., their Γ value is below one) assuming a Kroupa/Chabrier IMF. In both objects, we assumed solar metallicities. However, if the metallicity were significantly below solar, this could lead to an overestimate of the population mass estimates; this seems possible in VUCD3 where the existing measurements span a wide range from [Z/H]=-0.28 to 0.35 (Evstigneeva et al. 2007; Firth et al. 2009; Francis et al. 2012). Globular cluster dynamical mass estimates also seem to be lower than expected (Strader et al. 2011; Kimmig et al. 2015), although this appears to be in part because of mass segregation within the clusters combined with the assumption of mass-traces-light models (Shanahan & Gieles 2015; Baumgardt 2017). However, we note that no mass segregation is expected in any of the UCDs with BHs due to their long relaxation times (the half-mass relaxation times are 203 Gyr, 624 Gyr, and 350 Gyr in VUCD3, M59cO and M60-UCD1 respectively). The most massive globular clusters also have long relaxation times, and these clusters also seem to have lower than expected M/Ls 38 Figure 2.8: The dynamical-to-stellar mass ratio Γ vs. total dynamical mass. Grey points represent globular clusters (GCs) and UCDs with mass estimates based on integrated dispersions and assuming mass-traces-light models from Mieske et al. (2013) and references therein (with the exception of UCD3; Frank et al. 2011). Stars represent seven UCDs and two GCs for which we have AO-assisted stellar kinematic data in hand. The colored stars represent the new stellar mass measurements after accounting for a central massive BH. The arrows show the change caused in the stellar mass estimates by including the BH. 39 for clusters with [Fe/H]>-1 (Strader et al. 2011). Both these clusters and the less massive metal-rich clusters in the Milky Way with dynamical mass estimates based on N-body models by Baumgardt (2017) have masses 70-80% of the values expected for a Kroupa IMF, consistent with all three UCDs we have measured so far. We also note that Γ (assuming a Chabrier IMF) was recently found to be significantly below unity in the nucleus of NGC 404 (Nguyen et al. 2016) and in many compact elliptical galaxies (Forbes et al. 2014). 2.5.2 Impact of Anisotropy Due to the intrinsic degeneracy between the BH mass, stellar M/L, and anisotropy parameter, we also tested the impact of including anisotropic orbits in our JAM models. These degeneracies are represented by the contour plots shown in Figure 2.7. From left to right, top to bottom, the panels represent VUCD3 anisotropy vs. BH mass, VUCD3 Γ vs. BH mass, M59cO anisotropy vs. BH mass, and M59cO Γ vs. BH mass. The blue points represent our grid sample, and the green point is the best-fit determined over the entire grid. The colored lines represent the best-fit anisotropy and Γ assuming the BH mass makes up 1% (green), 5% (orange), and 10% (yellow) of the total dynamical mass of the system. The contours were calculated by determining the minimum chi-squared value between the four free parameters for each pair of grid points shown in the plot (i.e., for each pair of grid points shown in Figure 2.7, we marginalised out the two parameters not shown). Here, it is clear that the BH mass scales inversely with both the anisotropy and M/L. We note that the green points in Figure 2.7 show the best-fit BH mass and Γ determined over the entire grid are consistent with the results we obtained when we assumed isotropy. For the kinematic data to be consistent with no BH, the anisotropy parameter needs to be as high as βz = 0.4 and βz = 0.6 for VUCD3 and M59cO, respectively (shown as a grey line in Figure 2.6). This would require both of these objects to have a high degree of radial anisotropy. However, we recognize that a lower value for the mass of the BH could lead to a more reasonable value for βz . Therefore, we also tested what the best-fit values for βz and Γ would be assuming the mass of the BH was 1%, 5% and 10% of the total dynamical mass. The best-fit values are represented by the green (1%), orange (5%) and yellow (10%) colored lines in Figure 2.7, and shown again as dispersion profile fits in Figure 2.9. In each case, the dynamical models fit the dispersion profile well, but βz at each BH mass remains high at 0.3 − 0.4 for VUCD3 and 0.5 − 0.6 for M59cO, similar to the no BH case. As discussed at the beginning of Section 2.5.1, high βz would be at odds with existing nuclei and the one previous UCD where this measurement has been made. Furthermore, Γ remains elevated in the case of VUCD3. The best-fit no BH mass stellar 40 Figure 2.9: Dispersion profiles of VUCD3 (left) and M59cO (right) where black points are the measured velocity dispersions. The green, orange, and yellow lines represent the best-fit anisotropic models to the dispersion profile assuming the mass of the BH is 1%, 5%, and 10% of the total dynamical mass, respectively. For VUCD3, the best-fit βz and Γ values are 0.4 and 1.8 (green), 0.4 and 1.55 (orange), 0.3 and 1.25 (yellow). For M59cO, the best-fit βz and Γ values are 0.6 and 0.8 (green), 0.6 and 0.7 (orange), 0.5 and 0.6 (yellow). The grey point in the M59cO dispersion profile was not fitted. M/Ls would be a factor of 1.9 above the population estimate. This is inconsistent with the stellar population results in M60-UCD1 and local group globular clusters which fall below the stellar population estimates. We compared our anisotropy values to similar objects tested for anisotropic orbits. M60UCD1, the only other UCD with a known value, was shown to be nearly isotropic (Seth et al. 2014). Other compact objects such as the nuclear star clusters in the Milky Way, CenA, NGC 404, NGC 4244 and compact object M32 have also been shown to have nearly isotropic orbits (Verolme et al. 2002; Cappellari et al. 2009; Schödel et al. 2009; Hartmann et al. 2011; Nguyen et al. 2016). Therefore, we conclude that while it is possible to have highly anisotropic orbits, it seems unlikely. 2.6 Discussion and Conclusions In this paper we have tested the hypothesis that the existence of central massive BHs making up ∼ 10% of the total mass of UCDs can explain the elevated dynamical-to-stellar mass ratios observed in almost all UCDs above 107 M (Haşegan et al. 2005; Mieske et al. 41 2008, 2013; Seth et al. 2014). For our analysis, we observed two Virgo UCDs, VUCD3 and M59cO, using adaptive optics assisted kinematics data from the Gemini/NIFS instrument combined with multiband HST archival imaging. The Gemini/NIFS data were used to determine radial dispersion profiles for each object. We found integrated dispersion values of 39.7±1.2 km s−1 for VUCD3 and 31.3±0.5 km s−1 for M59cO with central dispersion values peaking at 52.9 km s−1 and 40.2 km s−1 for each object, respectively. The HST archival images were fitted with a double Sérsic profile to model the mass density, total luminosity, and to test for the presence of stellar population variations. We found a total luminosity of LF 814W = 17.8 × 106 L and LF 475W = 20.3 × 106 L for VUCD3 and M59cO, respectively. Both objects showed a mild positive color gradient as a function of radius, implying multiple stellar populations. These effects were accounted for in our mass models by multiplying the luminosity by the M/L determined from SSP models. Combining our mass models and velocity dispersion profiles, we created dynamical models using JAM. We found that the best-fit dynamical models contained central massive +2.5 6 BHs with masses and three sigma uncertainties of 4.4+2.5 −3.0 and 5.8−2.8 × 10 M for VUCD3 and M59cO, respectively, assuming isotropy. These BHs make up an astonishing ∼13% of VUCD3’s and ∼18% of M59cO’s total dynamical mass. The addition of a central massive BH has the effect of reducing Γ, as illustrated by the red and blue arrows in Figure 2.8. For comparison, the best-fit dynamical model, assuming isotropy, without a central BH returns a Γ value of 1.7 with a total dynamical mass of 66 × 106 M and 0.9 with a total dynamical mass of 83 × 106 M for VUCD3 and M59cO, respectively. The best-fit dynamical models reduce Γ to 0.8 with a total dynamical mass of 32 × 106 M for VUCD3 and 0.3 with total dynamical mass 32 × 106 M for M59cO. Due to the intrinsic degeneracy between the BH mass and anisotropy parameter, βz , in the JAM models, we also tested the impact of including anisotropic orbits. We found that βz values of 0.4 for VUCD3 and 0.6 for M59cO allow the kinematic data to be consistent with no BH. Furthermore, a central massive BH making up 1%, 5%, and 10% of the total dynamical mass also match the kinematic data, but leave βz relatively unchanged at ∼ 0.4 for VUCD3 and ∼ 0.6 for M59cO. Comparing these values with other nuclear star clusters and UCDs shows that highly radially anisotropic orbits in UCDs are improbable. We conclude that both VUCD3 and M59cO host supermassive BHs and are likely tidally stripped remnants of once more massive galaxies. We can estimate the progenitor mass 42 assuming that these UCDs follow the same scaling relations between BH mass and bulge or galaxy mass as unstripped galaxies (e.g., Kormendy & Ho 2013; van den Bosch 2016), as well as similar scaling relations for NSCs (e.g., Scott & Graham 2013; Georgiev et al. 2016). Mieske et al. (2013) used BH scaling relations to show that today’s UCDs are consistent with typically having ∼ 1% of the luminosity of their progenitor galaxy, suggesting progenitor galaxies for VUCD3 and M59cO of roughly ∼ 109 M . With the measured BH masses, we can use scaling relations to estimate more precise progenitor masses; using the Saglia et al. (2016) relations for all galaxies, we estimate dispersions of ∼100 km s−1 , and bulge masses of 1.2×109 and 1.7×109 M for VUCD3 and M59cO, respectively, with the scatter in the latter relationship suggesting about an order of magnitude uncertainty. We note the high inferred galaxy σ values are not necessarily expected to be observed in the nucleus (e.g., Koleva et al. 2011; Feldmeier et al. 2014). Assuming the inner components of the UCDs represent the nuclear star clusters of the progenitor (Pfeffer & Baumgardt 2013), the apparent magnitudes of these components are ∼19.5-20 in g with inferred masses of 4-10×106 M . Nuclei with these magnitudes and similar effective radii are seen in Virgo Cluster galaxies ranging from Mg of −16 to −18 (with galaxy masses of 1–8×109 M assuming M/Lg similar to that of M59cO; Côté et al. 2006). Both inner components are slightly bluer; if interpreted as being due to a younger age population (rather than lower metallicity), this suggests that they formed stars as recently as 6-8 Gyr (Section 2.3), which may suggest that they were stripped more recently than this. Measurements for both VUCD3 and M59cO suggest they have near-solar metallicity and are α enhanced (see Section 2.1). The metallicity seems consistent with present day nuclei in the luminosity range expected for the progenitors (Geha et al. 2003; Paudel et al. 2011). However, only a small fraction of Virgo NSCs seem to be significantly α enhanced (Paudel et al. 2011). Overall, the BH mass and NSC luminosity and metallicity suggest that the host galaxies for both objects were of order 109 M . This is about an order of magnitude lower than the likely progenitor of M60-UCD1 (Seth et al. 2014), and is in a galaxy mass regime where very few BH masses have been measured (Verolme et al. 2002; Seth et al. 2010; Reines et al. 2013; den Brok et al. 2015; Nguyen et al. 2016). In systems with measured BH masses, the ratio of BH to NSC masses ranges from 10−4 to 104 (Georgiev et al. 2016), thus it is consistent with the measurements here of roughly equal NSC and BH masses. These UCDs constitute the second and third UCDs known to host supermassive BHs. All UCDs with adaptive optics kinematic data available thus far have been shown to host central massive BHs. After taking these BHs into account, the stellar mass of UCDs is no 43 longer higher than expected, suggesting other UCDs with high Γ may host BHs (Figure 2.8). Nondetection of BHs in two UCDs based on ground-based data have been published. In NGC4546-UCD1 (M? ∼ 3 × 107 M ), Norris et al. (2015) suggests that any BH is .3% of the stellar mass despite finding evidence that this UCD is in fact a stripped nucleus. This result depends on the assumption of a stellar M/L based on the age estimate of the stellar population; a lower stellar M/L such as those we find (Figure 2.8) would result in a higher possible BH mass. Another BH nondetection was reported by Frank et al. (2011) using isotropic models of the bright and extended UCD3; a 5% mass fraction BH is consistent with their data within 1σ, while a 20% BH mass fraction is excluded at 96% confidence. Our high resolution results reinforce the hypothesis that UCD BHs could represent a large increase in the number density of massive BHs (Seth et al. 2014). Simulations of tidal stripping from cosmological simulations suggest that all high-mass UCDs (>107.3 M ) are consistent with being stripped nuclei (Pfeffer et al. 2014, 2016), with a mix of globular clusters and stripped nuclei at lower masses. Depending on how common stripped nuclei are, these objects may represent the best way of studying the population of BHs in lower mass galaxies, a critical measurement for understanding the origin of supermassive BHs (Volonteri 2010). This emphasizes the value in making similar studies of nearer, lower mass UCDs. For example, some local group globular clusters are also thought to be tidally stripped remnants (e.g., ω Cen, G1 Norris et al. 1997; Meylan et al. 2001). BH detections have been claimed in ω Cen (e.g., Noyola et al. 2010; Baumgardt 2017), M54 (Ibata et al. 2009), and 47 Tucanae (Kızıltan et al. 2017), but these remain controversial (van der Marel & Anderson 2010; Haggard et al. 2013). A BH has also been claimed in the Andromeda globular cluster G1 (Gebhardt et al. 2005), but accretion evidence for this BH has been elusive (Miller-Jones et al. 2012). In all these cases, the mass fraction of the black hole is certainly lower than the mass fractions of >10% that we find here. The lack of knowledge of BH demographics in low-mass host galaxies prevents easy comparison with nonstripped systems. Nonetheless, nearby UCDs represent the best place to push towards lower masses; we have ongoing observing programs for six additional UCDs, including objects in M31 and NGC 5128. CHAPTER 3 THE BLACK HOLE IN THE MOST MASSIVE ULTRACOMPACT DWARF GALAXY M59-UCD3 This chapter originally appeared in the 2018 Astrophysical Journal (Ahn et al. 2018). In that format it is c 2018 by the American Astronomical Society. It is reprinted here with their permission. Sections, tables and figures have been renumbered and reformatted to match the overall flow of this thesis. The coauthors of the paper are Anil C. Seth, Michele Cappellari, Davor Krajnović, Jay Strader, Karina T. Voggel, Jonelle L. Walsh, Arash Bahramian, Holger Baumgardt, Jean Brodie, Igor Chilingarian, Laura Chomiuk, Mark den Brok, Matthias Frank, Michael Hilker, Richard M. McDermid, Steffen Mieske, Nadine Neumayer, Dieu D. Nguyen, Renuka Pechetti, Aaron J Romanowsky, and Lee Spitler. 3.1 Introduction The past few decades have seen the emergence of compact stellar systems that have blurred the conventional lines, based on properties such as mass, luminosity and size, between star clusters and galaxies (Hilker et al. 1999; Drinkwater et al. 2000). One such classification of objects, the ultracompact dwarf galaxies (UCDs), occupy the region between globular clusters (GCs) and compact ellipticals (cEs) (e.g., Brodie et al. 2011; Misgeld & Hilker 2011; Pfeffer & Baumgardt 2013; Norris et al. 2014; Janz et al. 2016). With masses and radii of M > 2 × 106 M and r > 10 pc, UCDs are among the densest stellar systems in the Universe. However, the nature and origin of these dense objects is still widely debated. Early interpretations suggested that UCDs could be the most massive GCs (e.g., Mieske et al. 2002; Fellhauer & Kroupa 2002, 2005; Kissler-Patig et al. 2006; Murray 2009) or possibly the tidally stripped remnants of dwarf galaxies (Bekki et al. 2001, 2003; Drinkwater et al. 2003; Strader et al. 2013; Pfeffer & Baumgardt 2013; Forbes et al. 2014). However, there is evidence that both formation mechanisms could contribute to the UCD population we observe (Brodie et al. 2011; Da Rocha et al. 2011; Norris & Kannappan 2011; Janz et al. 45 2016). In the last decade, observational results based on structural information from Hubble Space Telescope (HST) imaging combined with integrated velocity dispersion measurements have shown an interesting trend: UCDs dynamical mass-to-light ratios (M/Ldyn ) appear to be systematically elevated when compared to the canonical stellar population estimates (Haşegan et al. 2005; Dabringhausen et al. 2008; Mieske et al. 2008; Taylor et al. 2010; Frank et al. 2011; Strader et al. 2013). These results prompted suggestions of variations in the stellar initial mass function (IMF) of UCDs (top-heavy: Murray 2009; Dabringhausen et al. 2009, 2010; bottom-heavy: Mieske & Kroupa 2008). Further explanations have suggested that these elevated M/Ls could be explained by ongoing tidal stripping (Forbes et al. 2014; Janz et al. 2016), or, as a relic of a massive progenitor galaxy in the tidal stripping scenario, a central massive black hole (BH) making up ∼ 10−15% of the total mass (Mieske et al. 2013). Supermassive black holes (SMBHs) have been confirmed in four UCDs with masses M > 107 M ; three in the Virgo cluster (Seth et al. 2014; Ahn et al. 2017), and one in the Fornax cluster (Afanasiev et al. 2018). A search for SMBHs in two lower mass (M < 107 M ) UCDs in Centaurus A yielded a non-detection (Voggel et al. 2018). However, Voggel et al. (2018) also showed that the dynamical-to-stellar M/Ls were overestimated in previous studies. The combination of this evidence still supports the idea that most UCDs with apparently high dynamical-to-stellar mass ratios (including a vast majority of UCDs above 107 M ) host SMBHs. Lower mass UCDs may be the high mass end of the globular cluster distribution. This view would be consistent with the analysis of the stripped nuclei contribution to UCDs in ΛCDM simulations by Pfeffer et al. (2014) and Pfeffer et al. (2016). Despite the fact that all of the detected SMBHs are found in massive (M > 107 M ) UCDs, the most massive UCD discovered to date, M59-UCD3 (M∗ ∼ 2 × 108 M re ∼ 25 pc), has been left out. This is in part due to its recent discovery (Liu et al. 2015; Sandoval et al. 2015), and the lack of high resolution imaging data needed for dynamical modeling. Thus, M59-UCD3 serves as an important test of the idea that the most massive UCDs host SMBHs. In this paper, we present the dynamical modeling techniques and results for M59-UCD3. An image of M59-UCD3 and its host galaxy (M59=NGC 4621) is shown in Figure 3.1. M59-UCD3 is located 10.2 kpc in projection from the center of M59, assuming an average distance of 16.5 Mpc to the Virgo Cluster based on surface brightness fluctations (Mei et al. 2007). We note that the individual distance of M59 has been measured to be 14.9 ± 0.4 Mpc (Mei et al. 2007), and this distance has been used in previous luminosity and mass 46 M59cO 10.2 kpc M59 (Host galaxy) M59-UCD3 (Ultracompact Dwarf) Figure 3.1: The M59/M59-UCD3 system discussed in this paper. Here, the main image shows the 2M ASS LGA image (Jarrett et al. 2003). M59-UCD3 is outlined in the yellow box. The inset image is a zoom-in HST image taken through the F814W filter on the WFC3 instrument. We also outline another UCD in this system, M59cO, in blue. The red line connecting the UCD to the host galaxy shows the projected distance assuming M59-UCD3 is at a distance of 16.5 Mpc. estimates; our assumed 16.5 Mpc will yield a 10% higher dynamical mass estimate relative to the previous mass determination (Liu et al. 2015), while at 16.5 Mpc M59-UCD3 has a measured MV = −14.8 (Sandoval et al. 2015). We adopt the conventional definition of Γ ≡ (M/L)dyn /(M/L)∗ , which is the ratio between the dynamically determined total M/L and the stellar M/L inferred from stellar population modeling. Throughout this paper, we assume a Chabrier IMF for the stellar population models. The metallicity of M59-UCD3 has been estimated to be near solar, with [Fe/H]∼ −0.01 and [α/Fe]∼0.21 (Liu et al. 2015; Sandoval et al. 2015; Janz et al. 2016; Villaume et al. 2017). These values of near solar metallicity and moderate alpha-element enhancement are consistent with previously measured high mass UCDs (Evstigneeva et al. 2007; Chilingarian & Mamon 2008; Firth et al. 2009; Francis et al. 2012; Strader et al. 2013; Janz et al. 2016). All magnitudes are reported in the AB magnitude system. Fur- 47 thermore, all magnitudes and colors have been extinction corrected using AF 475W = 0.100 and AF 814W = 0.052 (Schlafly & Finkbeiner 2011). This paper is organized as follows: Section 3.2 discusses the data used for analysis, how we determined a density profile, and how the kinematics were modeled. In Section 3.3, we present our three dynamical modeling techniques and the results from each. Section 3.4 discusses the radio/X-ray observations of UCDs and whether these observations can be used to infer the presence or not of an accreting SMBH. In Section 3.5, we discuss the implications of the results and present our conclusions. 3.2 Data and Methods In this section, we present the data and our reduction techniques. Section 3.2.1 discusses the HST images and our methods for deriving a mass model. Section 3.2.2 explains the reduction of our Gemini/NIFS integral field spectroscopy and the derivation of the kinematics. 3.2.1 Imaging Data and Deriving a Mass Model We obtained images of M59-UCD3 from the HST GO Cycle 23 program 14067 (PI: Ahn) with the Wide Field Camera 3 (WFC3) instrument, which has a pixel scale of 0.0400 pixel−1 . Our data were taken through the F475W and F814W filters. The exposure times in each filter were 1470s and 747s for F475W and F814W, respectively. We derived a point spread function (PSF) for each filter following the procedure outlined in previous studies (Evstigneeva et al. 2007; Ahn et al. 2017). To briefly summarize, we generated the distorted PSF with TinyTim and placed these PSFs in an empty copy of the raw HST flat fielded image at the location of our observed target. The distorted PSFs were then passed through MultiDrizzle using the same parameters as were used for the data. This produces model PSFs that are processed in the same way as the original data. The background (sky) level is traditionally determined from empty portions of the image. However, UCDs generally fall within the stellar halo of their host galaxy. Therefore, the sky level is not uniform across the image. To account for this, we added the MultiDrizzle level, subtracted by the HST reduction pipeline, back in and modeled the sky as a tilted plane. This was accomplished by masking all foreground/background objects including our UCD in the image. The good pixels (determined from the DQ extension of the image) were then weighted by their corresponding errors. Finally, a plane was fitted to the image to represent the sky level. The formal uncertainties on the sky level determination in this method are negligible. However, a clear systematic effect is seen in that the mean value of 48 the data minus sky model is offset from zero. We regard this as indicative of the systematic uncertainty, which are 0.86 counts in F475W and 1.38 counts in F814W. We use these uncertainties for plotting purposes only in the surface brightness profile (Figure 3.2 as grey bands), and color profile (Figure 3.3 as our error bars on the data). To enable dynamical modeling of M59-UCD3, we needed to create a model to represent the luminosity and mass distribution. Typically, in compact objects such as UCDs, the mass is assumed to trace the light (e.g., Mieske et al. 2013; Seth et al. 2014). However, previous studies have found significant color gradients in UCDs, suggesting multiple stellar populations (Chilingarian & Mamon 2008; Evstigneeva et al. 2008; Ahn et al. 2017). Figure 3.2: Surface brightness profile of M59-UCD3 in HST/F814W, which was used for dynamical modeling. Black stars are data, cyan lines are convolved double Sérsic profile models, yellow lines are convolved triple Sérsic profile models, the red line is the triple Sérsic reconstructed profile, and green, blue, and purple lines are the individual Sérsic components. The gray bands represent the uncertainty in our background sky determination. The residuals between the data and convolved models are shown in the bottom panel. 49 Figure 3.3: Color profile of M59-UCD3 is shown as black diamonds. The error bars are calculated from the uncertainty in our background (sky) level determinations. The solid lines indicate the triple-component Sérsic model fits that have been convolved with the HST PSF. Dashed lines show models that are unconvolved. The colors represent whether the shape parameters of the Sérsic profiles were independent (black) or fixed (red and blue). Blue lines indicate that the shape parameters of the F475W filter were held fixed to the best-fit F814W Sérsic models, and red lines are vice-versa. The unconvolved fixed models (red and blue) provide a well defined color for each Sérsic component. Our default model is shown in red. Here, the inner, middle, and outer colors are 1.26, 1.32, and 1.06 mag for our default model, respectively. See Section 3.2.1 for a discussion on our choice of the default model. Therefore, two filter data are essential for determining the most accurate luminosity and mass profiles of UCDs. The uncertainties in the luminosity and mass profile combinations are discussed in Section 3.3.1. For now, we discuss the general procedure for determining our luminosity and mass distributions. The surface brightness profile was determined by fitting the data in each filter to a PSFconvolved, multiple component Sérsic profile using the two-dimensional fitting algorithm, GALFIT (Peng et al. 2002). The parameters of the individual Sérsic profiles that were fitted include: the total magnitude (mtot ), effective radius (Re ), Sérsic exponent (n), position angle (P A), and the axis ratio (q). The fitting was done in two ways, similar to our previous study 50 (Ahn et al. 2017). In short, we fitted while allowing all of the above parameters to vary, henceforth referred to as the “free” fit. The initial fits showed an isophotal twist between the individual Sérsic profiles. However, the axis ratios of the outer profiles were nearly circular (q ∼ 0.99). Furthermore, two of the three dynamical modeling techniques are restricted to axisymmetric potentials and thus, do not allow for isophotal twists. To enable comparison between all three techniques, we fixed the axis ratio of the outer profiles to be perfectly circular, and fitted the data again. Next, we fitted the data while fixing Re , n, P A, and q to the best-fit model from the other filter, which we call the “fixed” fit. For example, the fixed F814W fit contains all of the shape parameters from the best-fit F475W model, where only the total magnitude is varied. Since the only free parameter is the total magnitude, these fits provide a well defined color for each Sérsic profile. The Sérsic profiles used to create the default luminosity and mass models are shown in Figure 3.2, and the parameters of the best-fit models are shown in Table 3.1. We chose the default model to be the fixed F814W fit (outlined in bold in Table 3.1) due to 1.) the ability to accurately reproduce the surface brightness profile, 2.) it clearly provides the best fit to the color profile (discussed below), and 3.) it provides a well defined color for each Sérsic component. However, as discussed in Section 3.3.1, the choice of the luminosity and mass model produces a minor effect on the results of the dynamical models. The most massive (M > 107 M ) UCDs have been found to consist of two components: a dense central component, and a more diffuse extended component, as shown by the two-component profile fits in previous studies (Evstigneeva et al. 2007, 2008; Chilingarian & Mamon 2008; Strader et al. 2013; Ahn et al. 2017; Afanasiev et al. 2018; Voggel et al. 2018). Shown in cyan in Figure 3.2, a two-component Sérsic profile provides a suitable fit to the data within the central 2.500 . However, a much better fit out to 300 was obtained by allowing a three-component Sérsic profile, shown in Figure 3.2 as yellow. The outer component is rounder, and as we will see below appears to be bluer than the inner components. The inner Sérsic components can also be replaced by a single King profile; a two-component King + Sérsic profile provides an equally suitable fit as the three-component Sérsic out to 300 . The central King model component in this fit has a core radius Rc = 0.0500 (4 pc) and a concentration parameter c = 1.42. We stick with our three-Sérsic model fit for the remainder of the paper due to the ease of transforming this model into our mass models. The simulations of Pfeffer & Baumgardt (2013) show that when galaxies are stripped into UCDs, they show a two-component structure with the inner component consisting of the galaxies’ nuclear star cluster (NSC) surrounded by the remains of the tidally stripped galaxy. 51 Table 3.1: Best-fit Sérsic Parameters of M59-UCD3 Component Inner Middle Outer 1 Parameter χ2 mtot mtot Fixed1 Re [00 ] Re [pc] n q P A [◦ ]2 mtot mtot Fixed1 Re [00 ] Re [pc] n q mtot mtot Fixed1 Re [00 ] Re [pc] n q F 475W 1.045 17.27 17.58 0.21 16.8 2.04 0.74 -6.41 18.02 17.86 0.50 40.0 0.80 1.00 18.84 18.20 1.24 99.2 0.90 1.00 F 814W 1.356 16.24 16.01 0.26 20.8 1.72 0.73 -6.36 16.70 16.70 0.35 28.0 5.23 1.00 16.95 17.78 0.56 44.8 0.90 1.00 The “Fixed” magnitudes show the total magnitude when the shape parameters of the Sérsic profiles are held fixed to the other filter. The PA orientation is N=0◦ and E=90◦ . The default model (in bold) was constructed by fitting F475W first, allowing all parameters to be free. Then the shape parameters from the previous fit were held fixed, and only mtot is fit for F814W. Next, we tested fitting F814W first, allowing all parameters to be free. Then the shape parameters from the previous fit were fixed, and only mtot for F475W is fit. The numbers from the second approach are not bolded. 52 Therefore, the multiple component model that is required to fit the surface brightness profile (regardless of the model choice) provides evidence that M59-UCD3 is likely a tidally stripped remnant. The total luminosity in each band (within the central 300 ) and effective radius calculated from the unconvolved Sérsic profiles are found to be LF 814W = 1.18 − 1.19 × 108 L , LF 475W = 6.1 − 6.2 × 107 L , and Re = 26 − 29 pc or 0.32 − 0.3600 , respectively. Here the ranges are quoted based on the total luminosity and effective radius calculated across all free and fixed models. Our effective radii are consistent with previous measurements. However, our calculated luminosity is slightly lower than the previously estimated Lg = 9.5 × 107 L (Liu et al. 2015; Sandoval et al. 2015). We attribute these deviations in the luminosity to our deeper HST imaging data, where previous studies were limited to ground based data only. Each of the best-fit Sérsic profiles were then parameterized by a multi-Gaussian expansion (MGE: Emsellem et al. 1994; Cappellari 2002), using the MGE FIT 1D fitting method and code1 of Cappellari (2002) for use in the dynamical modeling. Creating dynamical models that reproduce the observed kinematics requires a luminosity and mass profile. To determine the mass profile, we make use of our dual filter HST data to test for the presence of stellar population variations. M59-UCD3 shows a complicated color profile that varies by 0.25 mags within the central 300 , as shown in Figure 3.3. Here, the diamonds represent the data, where the error bars are calculated from our estimated systematic background effects, discussed above. Solid lines represent the convolved model and dashed lines represent the unconvolved model. The colored lines show whether the parameters of the Sérsic fits were independent (black) or fixed (blue and red). Blue lines indicate that the shape parameters from the F475W filter were held fixed to the F814W filter, and red lines are vice versa. It is clear, from the unconvolved models, that the bluest colors near the center are due to PSF effects and there is a gradient toward bluer colors at larger radii. This is unique when compared to other UCDs, which generally show either no overall trend or a gradient toward redder colors are larger radii (Chilingarian & Mamon 2008; Evstigneeva et al. 2008; Seth et al. 2014; Janz et al. 2015; Ahn et al. 2017). The unconvolved fixed models (dashed blue and red lines in Figure 3.3) provide a well defined color for the three Sérsic components. However, the fixed F814W Sérsic shape parameters (red lines) provide a better fit to the color profile, which motivated our choice of the default model. Here, the inner component color is F475W − F814W = 1.26, the middle component is F475W − F814W = 1.32, and the outer component color is F475W − F814W 1 http://purl.org/cappellari/software 53 = 1.06. We used these colors combined with the Bruzual & Charlot (2003) Padova 1994 simple stellar population (SSP) models and corresponding code2 , assuming solar metallicity and a Chabrier IMF, to determine the mass-to-light ratio (M/L). The default Bruzual & Charlot (2003) color tables and corresponding M/Ls do not include our HST filter set. Therefore, we downloaded the filter transmission curves and reran the composite stellar population model code. We found the inner, middle and outer M/LF 814W,∗ to be 2.5 ± 0.4, 2.9 ± 0.4, and 1.0 ± 0.1, with corresponding ages of 9.9, 13.7, and 2.8 Gyrs, respectively. Here, the error bars are calculated assuming an uncertainty of ±0.05 mags in our color determinations. We also determined the total mass by multiplying the total luminosity with the corresponding M/L for each Sérsic component. We found the inner component total stellar mass to be 17.3 ± 2.8 × 107 M , middle component to be 10.6 ± 1.5 × 107 M , and the outer component to be 1.2 ± 0.1 × 107 M . The mass density profile was then determined by multiplying the luminosity model MGEs by their corresponding M/Ls. We note that to test the systematic effects of our choice of the default mass profile, we also determined two additional mass density profiles by 1.) computing the M/L from the color of the free Sérsic profile fits (black lines in Figure 3.3) at the FWHM of each Gaussian component in the MGE luminosity profile and 2.) by assuming a mass-follows-light model, which is equivalent to creating a mass density profile where all of the Sérsic components have been scaled by the flux weighted M/L described below. The latter test was motivated by the bluer outer most component, which could also be interpreted as a metal-poor old stellar population. As discussed in Section 3.3.1, this results in a M/LF 814W,∗ = 1.9, which is intermediate to our mass-follows-light model. The luminosity is used to calculate the center of each kinematic bin and create the observed kinematic field. Since our kinematic data were taken in the K-band, we determined a luminosity MGE in that band using the color profiles and SSPs. This was accomplished by creating a color-color diagram of F814W −K vs. F475W − F814W from the SSP models. We then used our derived color to infer the F814W −K color. For each Sérsic component, we found these colors to be F814W −K = 2.28, 2.32, and 2.12 for the inner, middle, and outer components, respectively. These colors lead to a scale factor in the luminosity surface density for each component of 2.59 (inner), 2.69 (middle), and 2.24 (outer). These scale factors were multiplied by the luminosity profile for each component to make our K-band MGEs used in the dynamical models. The best-fit model MGE, for our default model, is shown in Table 3.2. 2 http://software.astrogrid.org/p/cea/latest/cec/config/galaxev/GALAXEV.html 54 Table 3.2: Default Multi-Gaussian Expansion (MGE) Used in the Dynamical Modeling. M ass (M pc−2 )1 294397. 361060. 370914. 308330. 199887. 99196.3 36071.5 9451.14 1737.86 210.794 11.8194 2547.19 6801.10 8380.02 3882.25 427.301 44.2053 116.224 173.015 128.835 40.5159 3.46218 1 IK (L pc−2 )2 304328. 373239. 383426. 318732. 206630. 102542. 37288.3 9769.96 1796.49 217.905 12.2182 2355.68 6289.76 7749.98 3590.37 395.175 101.759 267.543 398.277 296.574 93.2663 7.96984 σ(00 ) 0.001 0.004 0.011 0.027 0.059 0.115 0.209 0.356 0.576 0.903 1.459 0.055 0.180 0.355 0.541 0.738 0.073 0.263 0.601 1.050 1.566 2.175 q 0.74 0.74 0.74 0.74 0.74 0.74 0.74 0.74 0.74 0.74 0.74 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 PA (◦ ) -0.11 -0.11 -0.11 -0.11 -0.11 -0.11 -0.11 -0.11 -0.11 -0.11 -0.11 -56.78 -56.78 -56.78 -56.78 -56.78 -70.25 -70.25 -70.25 -70.25 -70.25 -70.25 The M/L was used to determine the mass profiles. These methods are described in Section 3.2.1. 2 The luminosity MGEs were created in the K-band, which required an assumption on the absolute magnitude of the sun. We assumed these values to be 4.523 in F 814W and 3.29 in K taken from http://www.baryons.org/ezgal/filters.php. See Section 3.2.1 for a more detailed explanation of how the luminosity profile was derived. The PA adopted for all of the dynamical models was -6.41◦ where N=0◦ and E=90◦ . The horizontal lines separate the individual Sérsic models 55 Finally, the SSP models and color profiles were used to determine a flux weighted average M/L. We computed this by determining the flux within the central 300 from model images of the inner, middle, and outer Sérsic profiles, and then weighting the individual M/Ls calculated above by their corresponding flux. We found the average M/LF 814W,∗ to be 2.47 ± 0.25 (assuming ∼10% uncertainties similar to individual M/Ls). We calculated the overall F475W − F814W color to be 1.26. Using the SSP models, described above, we estimate V − F814W = 0.86 (Bruzual & Charlot 2003). This corresponds to M/LV,∗ = 4.2 ± 0.4, which is consistent with previous results (Liu et al. 2015). We note that the Mieske et al. (2013) eq. 4 estimate of M/LV,∗ for [Fe/H]=-0.01 is 4.07 based on the mean between the Bruzual & Charlot (2003) and Maraston (2005) model predictions. 3.2.2 Kinematic Derivation Our spectroscopic data were obtained with the Near-Infrared Integral Field Spectrometer (NIFS) mounted on the Gemini North telescope using Altair laser guide star adaptive optics (Herriot et al. 2000; McGregor et al. 2003; Boccas et al. 2006). The observations were taken on the nights of 2015 May 4th, 5th, and 6th, in the K band at wavelengths from 2.0 µm to 2.4 µm. Gemini/NIFS data are taken in 0.100 × 0.0400 pixels over a 300 field of view with a spectral resolution of λ δλ ∼ 5700 (σinst = 22 km s−1 ). Our data were reduced following the same procedure as our previous studies, using the Gemini version 1.13 IRAF package (Seth et al. 2014; Ahn et al. 2017). To summarize, first, arc lamp and Ronchi mask images were used to determine the spatial and spectral geometry of the images. Next, the spectra were sky subtracted, flat-fielded, had bad pixels removed, and split into long-slit slices. Finally, the spectra were corrected for telluric absorption with an A0V telluric star taken on the same night using the NFTELLURIC procedure. The final data cube consisted of eight (five on May 4th, one on May 5th, and two on May 6th) 900s on-source exposures taken in either an object-sky-object or object-sky sequence. The object exposures were dithered to give independent sky measurements for each exposure and to improve the signal-to-noise ratio (S/N ). We created our own version of NIFCUBE and NSCOMBINE to rebin the spectra to a 0.0500 × 0.0500 pixel scale, and enable error propagation of the variance spectrum. Finally, the spectra were combined using our IDL program, which centroids the nucleus and rejects bad pixels based on the nearest neighboring pixels. The kinematic PSF was determined by convolving a Gauss+Moffat (1969) function with an HST model K band image to match the surface brightness in the kinematic data cube. The HST model K band image was derived following the same procedure outlined in 56 Section 3.2.1, where we inferred the F814W −K color using our derived color in each pixel combined with the SSP models. This model image was then convolved with a Gauss+Moffat function and fitted to the NIFS continuum image using the MPFIT2DFUN IDL code3 (Markwardt 2009). We note that in order to quantify the systematic effects of our PSF determination on the dynamical models, we also convoled the HST model image with a double Gaussian function. In this model PSF, the central Gaussian contained 58% of the light with a FWHM of 0.20100 , and the outer Gaussian contained 42% of the light with a FWHM of 0.89400 . Furthermore, we also determined a PSF where we fixed the central Gaussian FWHM in our Gauss+Moffat model to be the diffraction limit (F W HM = 0.0700 ) of the NIFS instrument. In this diffraction limited PSF, the central Gaussian contained 24% of the light, while the Moffat function contained the remaining 76% with a FWHM of 1.0800 . The effects of these PSFs are discussed in Section 3.3.1. In our best-fit model function the central Gaussian was found to contain 35% of the light with a FWHM of 0.16500 . The Moffat function contained the remaining 65% of the light with a FWHM of 1.0800 . The one dimensional analytic Gaussian+Moffat PSFs were again fitted by Gaussians using the MGE FIT 1D procedure of Cappellari (2002). The kinematics were measured by fitting the CO bandhead region (2.28 µm - 2.395 µm) to stellar templates using the IDL version of the penalized pixel algorithm pPXF, version 5.0.11 (Cappellari & Emsellem 2004; Cappellari 2017). We convolved the high-resolution λ ( ∆λ = 45000) Wallace & Hinkle (1996) stellar templates with the line spread function (LSF), determined in each bin using sky exposures, before fitting. The bins were created using the Voronoi binning method and IDL code1 to achieve S/N = 25 per resolution element (Cappellari & Copin 2003). However, beyond 0.500 the bins were remade using sectors to further increase the S/N at the largest radii and enable identification of rotation signatures. These outer bins have S/N ∼ 20, with the outermost four bins having S/N ∼ 12. We fit the radial velocity (V), dispersion (σ), skewness (h3 ), and kurtosis (h4 ) to the data. The uncertainties associated with the determined kinematics were estimated using a Monte Carlo simulation. Gaussian random noise was added to each spectral pixel in each bin, and the kinematics were fitted again. The standard deviations of the fits were taken as the one sigma uncertainties. An example of the integrated (r < 0.7500 ) kinematic fits is shown in Figure 3.4. The barycentric systemic velocity and integrated (r < 0.7500 ) dispersion of M59-UCD3 were found to be 434.5±0.6 km s−1 and 65.7±0.6 km s−1 , respectively. Previous measurements of the radial velocity have ranged from 373 ± 18 km s−1 to 447 ± 3 km s−1 3 http://purl.com/net/mpfit 57 Figure 3.4: Integrated (r < 0.7500 ) spectrum of M59-UCD3 is shown in black. The red line indicates the best kinematic fit, and residuals are shown in green. For visibility, the residuals are offset by 1.04 × 104 counts. The integrated dispersion was found to be σ = 65.7 ± 0.6 km s−1 with a median S/N = 57 per pixel. (Liu et al. 2015; Sandoval et al. 2015). Our measured radial velocity is within the range of previous measurements, but outside the quoted uncertainties. Our velocity dispersion is significantly lower than previous measurements, which have ranged from ∼ 70 km s−1 (no quoted uncertainty) to 77.8 ± 1.6 km s−1 (Liu et al. 2015; Janz et al. 2016). We attempted to simulate the Liu et al. (2015) data by creating a Jeans Anisotropic Model (JAM) with the best-fit parameters (described below), a pixel size equivalent to their slit width, and seeing FWHM= 1.8500 . Using these values, we found an integrated velocity dispersion of 67 km s−1 , which is still significantly lower than their measurement. The kinematic maps for the four measured velocity moments are shown in Figure 3.5 and their corresponding values are shown in Table 3.3. The velocity map shows clear rotation, and h3 shows the common anticorrelation with the velocity. We also see a very definitive 58 -0.5 0.5 -24 -1.0 -1.0 -0.5 -0.5 0.0 0.0 RA Offset ["] 0.5 1.0 0.5 1.0 1.0 -40 85 77 0.5 70 0.0 0.0 62 -0.5 -0.5 55 -1.0 1.0 -1.0 -1.0 -0.5 0.0 RA Offset ["] 0.5 1.0 0.5 0.0 -0.5 -1.0 RA Offset ["] 48 0.5 0.15 1.0 0.09 0.5 h3 0.03 0.0 0.0 -0.03 -0.5 -1.0 1.0 0.5 -0.5 -0.09 -1.0 -1.0 -0.5 -0.5 0.0 0.0 RA Offset ["] 0.5 0.5 1.0 1.0 1.0 -0.15 0.15 0.09 0.5 0.03 0.0 h4 -8 -0.5 Dec Offset ["] 0.0 1.0 Dec Offset ["] 8 0.0 -1.0 1.0 Dec Offset ["] 24 0.5 Dispersion [km s-1] Dec Offset ["] 0.5 RA Offset [pc] 50 0 -50 40 1.0 Relative Velocity [km s-1] RA Offset [pc] 50 0 -50 1.0 0.0 -0.03 -0.5 -0.5 -0.09 -1.0 1.0 -1.0 -1.0 -0.5 0.0 RA Offset ["] 0.5 1.0 -0.15 0.5 0.0 -0.5 -1.0 RA Offset ["] Figure 3.5: Full kinematic measurements of M59-UCD3, which include the radial velocity (top-left), dispersion (bottom-left), the skewness h3 (top-right), and kurtosis h4 (bottomright). Black contours show the K band continuum at 1 mag/arcsecond2 intervals. The median 1σ uncertainties are 3.9 km s−1 for the velocity, 3.7 km s−1 for dispersion, 0.04 for h3 , and 0.05 for h4 . Here, h3 clearly shows the common anticorrelation with the velocity. peak near the center of the dispersion map. The dispersion map has some significant asymmetries, with higher values at larger radii to the East of the nucleus. Using kinemetric fits4 (Krajnović et al. 2006), we find the rotation amplitude peaks at ∼30 km s−1 at a radius of 0.200 . The V /σ reaches ∼0.4 at this radius and stays constant at larger radii with both the dispersion and rotation decreasing outwards. This is similar to the V /σ of 0.3 − 0.5 seen in other UCDs including M60-UCD1 and VUCD3 (Seth et al. 2014; Ahn et al. 2017; Voggel et al. 2018). This level of rotation is higher than what is seen in any Milky Way globular clusters (e.g., Kimmig et al. 2015; Kamann et al. 2018), and within the range of nearby galaxy nuclei with resolved observations (e.g., Feldmeier et al. 2014; Nguyen et al. 2017). 3.3 Dynamical Modeling In this section, we present our dynamical modeling techniques and their results. Here, we present three techniques: JAM, axisymmetric Schwarzschild modeling, and triaxial Schwarzschild modeling (but using an axisymmetric shape and mass model). Previous 4 http://davor.krajnovic.org/idl/ X rad [pixels] ... -2.000 -1.000 0.000 ... Y rad [pixels] ... -4.866 -4.440 -4.000 ... Number of Pixels ... 3 2 1 ... χ2 ... 0.603 0.541 0.343 ... S/N ... 35.02 36.09 36.13 ... v [km s−1 ] ... 473.92 475.05 478.13 ... verr [km s−1 ] ... 2.93 3.36 4.53 ... σ [km s−1 ] ... 65.77 71.81 71.03 ... σerr [km s−1 ] ... 3.04 3.73 4.84 ... ... -0.04 -0.05 -0.05 ... h3 ... 0.04 0.04 0.05 ... h3,err ... -0.04 -0.03 -0.03 ... h4 ... 0.03 0.04 0.05 ... h4,err Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available. Bin Number ... 15 16 17 ... Table 3.3: Gemini/NIFS Kinematic Measurements of M59-UCD3 59 60 studies have shown that these techniques typically provide consistent results in estimating BH masses (e.g., Verolme et al. 2002; Cappellari et al. 2010; van den Bosch & de Zeeuw 2010; Seth et al. 2014; Drehmer et al. 2015; Feldmeier-Krause et al. 2017; Krajnović et al. 2018). We note that our dynamical modeling effort and interpretation for M59-UCD3 is more complicated than for any previous results on BHs in UCDs (Seth et al. 2014; Ahn et al. 2017; Afanasiev et al. 2018; Voggel et al. 2018). This is in part due to the high quality of the data, which enabled us to run triaxial Schwarzschild models in addition to the JAM models (as we also did in M60-UCD1; Seth et al. 2014), and also due to the lower BH mass fraction in M59-UCD3 relative to previous UCD detections. Our initial results showed the triaxial Schwarzschild and JAM codes gave inconsistent results for our BH mass, which led us to also run additional axisymmetric Schwarzschild models. We were unable to fully resolve the differences here, but do our best to weigh the evidence from these models, and in Section 3.3.4 conclude that there is evidence for a ∼2% mass fraction BH in M59-UCD3. We discuss each modeling method individually below before reaching these conclusions. 3.3.1 Jeans Anisotropic Models (JAM) We use the JAM method to fit the VRM S = √ V 2 + σ 2 using the code1 described in Cappellari (2008). This technique solves the anisotropic Jeans equations under two general assumptions: (1) The velocity ellipsoid is aligned with the cylindrical coordinate system (R, z, φ); (2) The anisotropy is constant and defined as βz = 1 − (σz /σR )2 , where σz is the velocity dispersion parallel to the rotation axis and σR is the velocity dispersion in the radial direction (Cappellari 2008). The primary input ingredients for the JAM models are the mass and luminosity density model (parameterized by MGEs as described in Section 3.2.1). These define the gravitational potential and the distribution of the tracer population, which produces the observed kinematics. We use the JAM code to fit the following parameters: (1) √ the intrinsic axial ratio, which acts as a parametrization of the q 02 −cos2 i inclination angle (q = , where q is the parameter we sample over and q 0 is the sin i axis ratio of the flattest MGE), (2) the mass of a point-like black hole MBH , (3) Γ, which parameterizes the M/L relative to the best-fit stellar population estimate, and (4) the anisotropy parameter βz (see §3.1 of Cappellari 2008 or §4 of Ahn et al. 2017 for a more detailed explanation of how these parameters are used). For a given set of these four parameters, the JAM model generates observable model kinematic measurements that can be compared with the measured values. For our JAM dynamical models we created a Markov Chain Monte Carlo (MCMC) 61 simulation to fully sample parameter space. The parameters sampled include: the axial ratio, anisotropy parameter (βz ), Γ, and the mass of the BH. We used the emcee python package5 developed by Foreman-Mackey et al. (2013), which is an implementation of the affine-invariant MCMC ensemble (Goodman & Weare 2010). This algorithm uses a set of walkers to explore the parameter space. At each step, the result of the likelihood for each walker informs the next choice of model parameters to be evaluated. We ran all of our models with 200 steps per walker. 3.3.1.1 Default Model For our default model, we used the K band luminosity model to represent the tracer population, mass model determined from the fixed F814W surface brightness profile fit to represent mass distribution, and the best-fit Gauss+Moffat function PSF (all described in §2). We ran our MCMC chain for this default model with 100 walkers for a total of 20000 steps. We consider the first 50 steps of each walker as the burn-in phase. Figure 3.6 shows the post burn-in phase distributions for this model. Here, the scatter plots show the two-dimensional distributions for each parameter with points colored according to their likelihood (white high and blue/black low). The histograms show the one-dimensional distributions for each parameter. We used the one-dimensional distributions to calculate the best-fit values and their corresponding uncertainties. The axis ratio is nearly unconstrained. The best-fit βz = 0.1 ± 0.1. We also see a slight degeneracy between βz and both Γ and the mass of the BH (MBH ). The degeneracy seen here is less significant than we found in our previous study (Ahn et al. 2017), which we attribute to our higher quality kinematic data. The best-fit Γ = 0.64 ± 0.02 (which corresponds to M/Ldyn,F 814W = 1.58 ± 0.05 or M/Ldyn,V = 2.69 ± 0.08), and MBH = (5.9 ± 1.1) × 106 M . The degeneracy between MBH and Γ is expected, well known, evident across all of our dynamical modeling techniques (shown in Figure 3.7, where blue contours represent the JAM models), and similar to what has been seen in previous studies (Seth et al. 2014; Ahn et al. 2017; Voggel et al. 2018). Here we quote the 1σ uncertainties calculated from the 16th and 84th percentiles of the one-dimensional MCMC distributions. Due to the lack of orbital freedom in the JAM models, we also quote the 3σ uncertainties (0.2 and 99.8 percentiles), which encompass the 0.3 , Γ = 0.64 ± 0.06 systematic effects, discussed below. Quoting 3σ errors we find βz = 0.1−0.2 (M/Ldyn,F 814W = 1.58 ± 0.15, M/Ldyn,V = 2.69 ± 0.25), and MBH = (5.9 ± 3.1) × 106 M . Figure 3.8 shows the comparison of the smoothed data with the model VRM S calculated 5 https://github.com/dfm/emcee 0. 2 0. 65 1e7 0. 6 0. 4 Γ 0. 2 0. 4 0. 6 0. 8 1. 0 0. 75 0. 70 0. 65 0. 60 0. 55 βz 0. 4 0. 2 0. 0 −0 .2 q 0. 60 0. 45 0. 30 0. 15 0. 2 MBH 0. 8 10. .5 05 0. 60 Γ 0. 70 0.−0 75.2 0. 0 βz 0. 4 62 MBH 1e7 Figure 3.6: MCMC post burn-in phase distributions for our default model. The scatter plots show the projected two-dimensional distributions for each parameter. The histograms show the projected one-dimensional distribution. From top-left to bottom-right the panels show the axis ratio, anisotropy parameter βz , Γ, and MBH . 63 Figure 3.7: Contour plots showing all three modeling techniques. Here, the blue contours represent the JAM models, red are the axisymmetric Schwarzschild models, and green are the triaxial Schwarzschild models. The black crosses denote the BH mass and Γ values used to make the VRM S comparison discussed below. 0.0 62 -0.5 -0.5 55 -1.0 1.0 -1.0 -1.0 -0.5 0.0 RA Offset ["] 0.5 0.5 0.0 -0.5 RA Offset ["] 1.0 -1.0 48 0.5 85 1.0 77 0.5 70 0.0 0.0 62 -0.5 -0.5 55 -1.0 1.0 -1.0 -1.0 -0.5 0.0 RA Offset ["] 0.5 0.5 0.0 -0.5 RA Offset ["] 1.0 VRMS [km s-1] Best-fit BH 77 0.5 70 0.0 1.0 Dec Offset ["] 0.5 85 1.0 VRMS [km s-1] Dec Offset ["] 1.0 48 -1.0 Figure 3.8: VRM S comparison between the smoothed data (left) and best-fit JAM model(right). Black contours show isophotes in the K band stellar continuum. We note that the we smoothed the kinematic data in this figure for visual comparison only. 64 with the best-fit parameters. 3.3.1.2 Quantifying the Systematic Effects To quantify the systematic effects of our choice of the default model (discussed above), we also ran MCMC chains with 100 walkers and 20000 total steps. To demonstrate the effects of each systematic effect, we varied the PSF, mass model, and luminosity model one by one, while holding the other parameters fixed to our default model. The cumulative likelihood for all model variations is shown in Figure 3.9. The black line represents the default model described above, with grey lines indicating the 1 and 3σ uncertainties. The colored lines show variations in the PSF (red), mass model (blue), and luminosity model Figure 3.9: Cumulative likelihood of MBH from the JAM modeling. The black line represents the best-fit of the default model. The colored lines represent variations in the PSF (red), mass model (blue), and luminosity model (cyan). The grey vertical lines indicate the best-fit, 1, and 3σ BH mass estimates based on the default model. See section 3.3.1 for a detailed explanation the individual red, blue, and cyan lines. 65 (cyan), which are explained in more detail below. The PSF is the largest source of systematic errors, seen in Figure 3.9 as the red lines. Our default PSF was the best-fit model to the NIFS continuum image from the K band model. The solid line represents a double Gaussian model PSF, and the dashed line represents a Gauss+Moffat function model PSF where the central Gaussian FWHM is assumed to have a FWHM = 0.0700 , which is the diffraction limit of the NIFS instrument. In this model PSF, the Moffat function FWHM was left unchanged from the default PSF, but the corresponding weights of the Gauss and Moffat functions were recalculated to match the NIFS continuum. In all of the dynamical modeling techniques, the model predictions are convolved with the kinematic PSF before comparison with the data (see Appendix A of Cappellari 2008). Therefore, our determination of the kinematic PSF is crucial in generating model observables that match the data. We tested the double Gaussian PSF as it is a common way to represent the adaptive optics PSF (but clearly did not fit the radial profile of the NIFS continuum as well in our case). We also tested a PSF created by convolving of the F814W image as opposed to our K band model and found very similar results (line not shown). The diffraction limited PSF was tested due to our inability to match the central few pixel VRM S values. As shown in Figure 3.10, the best-fit parameters with the diffraction limited PSF allow the model VRM S values to better match the data near the center, but have minimal effect on the BH mass. Here we took an elongated rectangular aperture along the semimajor (red) and semiminor (blue) axes. The solid lines represent the best-fit model VRM S with the default model, and the dashed lines represent the best-fit model with the diffraction limited PSF. Despite the better fit to the central data, we prefer our default model as the diffraction limited NIFS PSF, convolved with a model of the HST photometry, provides a worse fit to the surface brightness inferred from the NIFS datacube. The blue lines in Figure 3.9 represent mass model variations with the default luminosity model and PSF. In this case, the solid line shows a mass model that was determined from the color of the free Sérsic profile fits at the FWHM of each Gaussian in the MGE. This variation was motivated by the uncertainty in determining our mass profile from the fixed Sérsic models. The dashed line shows our mass-follows-light model, which is equivalent to scaling all of the Sérsic components by the flux weighted M/L (M/LF 814W,∗ = 2.47). This test was motivated by the bluer color of the outer component, which could also be due to an older, more metal-poor population at larger radii. The de-extincted color of this component is redder than the Bruzual & Charlot (2003) models for metallicities below Z=0.004 ([Fe/H]∼ −0.7). If we assume a model at that metallicity, we get a M/LF 814W,∗ =1.9; this is closer 66 Figure 3.10: VRM S comparison between the data, best-fit default JAM model (solid with MBH = 5.9 × 106 M , Γ = 0.64), and best-fit diffraction limited PSF JAM model (dashed with MBH = 6.3 × 106 M , Γ = 0.63). Here, we take an elongated rectangular aperture one pixel wide along the semimajor axis (red) and the semiminor axis (blue). The semimajor axis has been offset by 10 km s−1 for visibility. (but still lower than) the inner component M/Ls, and thus the resulting mass profile is intermediate between our constant M/L and default model M/Ls. Figure 3.9 shows that these mass model variations provide MBH constraints within the 1σ deviations from the default model, and therefore our results do not depend critically on the stellar population variations. Finally, the cyan line represents a luminosity model variation where we used the default mass model and PSF, but the original fixed F814W luminosity MGE. The luminosity model variation makes the least difference. This is expected since the luminosity model is only used to determine the center of each kinematic bin and generate the observed kinematic field. Figure 3.9 shows that our choice for the default model is reasonable given that all of the model variations fall within the 3σ error bars calculated from the default model 67 likelihood. 3.3.2 Axisymmetric Schwarzschild Models We fit the full line-of-sight velocity distribution (LOSVD) using an axisymmetric Schwarzschild orbit superposition model described in detail in Cappellari et al. (2006). This three-integral dynamical modeling technique is based on Schwarzschild’s numerical orbit superposition method (Schwarzschild 1979), which has been shown to reproduce kinematic observations (Richstone & Tremaine 1988; Rix et al. 1997; van der Marel et al. 1998). This method assumes axisymmetry, which also requires the potential to not vary on the time scale required to sample the density distribution of an orbit. Since the orbital time scale within M59-UCD3 is ∼ 106 yrs (assuming our effective radius and integrated dispersion), while the relaxation time is ∼ 1012 yrs and orbital time scale of M59-UCD3 around M59 is ∼ 108 yrs, the potential is unlikely to vary during the orbital sampling period. The generality of this method has allowed it to become the standard for determining the mass of central BHs when high resolution kinematic data are available (e.g., Cappellari et al. 2002; Verolme et al. 2002; Gebhardt et al. 2003; Valluri et al. 2005; Shapiro et al. 2006; van den Bosch et al. 2006; Nowak et al. 2007, 2008; Cappellari et al. 2009; Krajnović et al. 2009). However, the more general approach, which allows for triaxial systems, is described in van den Bosch et al. (2008), and discussed in Section 3.3.3. The full details of this method are described in Cappellari et al. (2006). In short, this method consists of four steps. First, as with the JAM models, a stellar potential is created by deprojecting the mass model MGEs assuming an axisymmetric shape and stellar M/L. Second, a representative, dithered orbit library is constructed with even sampling across the observable sampling space (based on the three integrals of motion and the stellar potential). Next, the orbits are projected onto the observable space using sky positions, taking into account the kinematic PSF and apertures (Voronoi bins, discussed in Section 3.2.2). Finally, the weights of each orbit are determined using a nonnegative least squares fit (Lawson & Hanson 1974), which are coadded to reproduce the observed kinematics. For our models, we follow the approach outlined in Krajnović et al. (2009), with the only exception being that for M59-UCD3; we do not assume mass follows light. We used our default model, described above, to construct the mass and luminosity profiles. Here, the mass MGE is used to calculate the orbit libraries. For these models, we fixed the inclinations angle to be 85◦ . This choice was arbitrary as the inclination angle has virtually no effect on the mass of the BH and Γ (see Section 3.3.1). We created a grid of the two free parameters: MBH and Γ. The orbit libraries are constructed for each MBH at an expected Γ and consist 68 of 21x8x7x2 orbital bundles, which are composed of 63 dithers (see Cappellari et al. 2006) This means there are 508032 total orbits of which 2352 are free to vary to optimize the fit. It is not necessary to compute an orbit library for every (MBH , Γ) combination because the orbit libraries can be scaled to match different Γs. For our grid, we sampled 32 MBH s between 2 × 104 M and 1.2 × 107 M , and 47 Γs between 0.43 and 0.89. The red contours in Figure 3.7 show the 1, 2 and 3σ contour results for Γ and MBH . These contours were calculated from the LOESS smoothed χ2 distribution. Likewise, the best fits determined 6 from the likelihood distribution are MBH = 2.5+1.8 −1.3 × 10 M , and Γ = 0.67 ± 0.03 (1σ uncertainties from the 16th and 84th percentiles). We note that these models are consistent with no BH within the 3σ contour. Figure 3.7 also shows the similar size and overlap between the 1σ uncertainties from the axisymmetric Schwarzschild models (red contours) and the JAM models (blue contours). However, we note that these are formal errors and are smaller for JAM due to the reduced freedom of the model. 3.3.3 Triaxial Schwarzschild Models Finally, we also fit the full LOSVD using the more general triaxial Schwarzschild models and corresponding code described in detail in van den Bosch et al. (2008). This model is also based on Schwarzschild’s numerical orbit superposition method (Schwarzschild 1979), but is not restricted to axisymmetry as described above. This method is implemented in a series of steps, similar to those described in Section 3.3.2. First, the stellar potential is created by deprojecting the mass model MGE, as described above. However, in the triaxial case, the viewing angles must be provided, which parameterize the intrinsic shape of the galaxy (see §3 of van den Bosch et al. 2008). Second, the initial conditions for each orbit library are found. These orbits must include all possible types of orbits that the potential can support (Thomas et al. 2004; van den Bosch et al. 2008). Next, the orbits are integrated for a fixed period of time, while storing the projected properties on a grid for comparison with the data. These properties are convolved with the same PSF as the kinematic observations. Finally, the orbital weights are determined using a sparse quadratic programming solver from the GALAHAD library, which is capable of fitting the kinematics in a least-squares sense while also satisfying the mass constraints (Gould et al. 2003). For M59-UCD3, we provided an oblate axisymmetric shape by specifying the viewing angles (θ,φ,ψ) = (85, -49.99, 89.99) degrees. In the oblate limit, the model does not depend on φ, while ψ must be 90◦ , and θ represents the inclination angle. Therefore, our model is nearly axisymmetric and seen at an inclination angle of 85◦ . The choice of the inclination 69 angle was again arbitrary but matches the axisymmetric Schwarzschild models. We sampled over all possible inclination angles and found consistent results. With this set up, the triaxial Schwarzcshild models are sampled in the axisymmetric limit but still allow for all possible orbits in a triaxial potential. For the other two parameters, Γ, and MBH , we created a grid similar to the axisymmetric Schwarzschild models described above. We sampled 47 MBH s ranging from 5.5 × 103 -2 × 107 M , and 62 Γs ranging from 0.43 to 1.04. The main results are shown as green contours in Figure 3.7. Here, the best-fit Γ = 0.75 ± 0.06, and we find MBH is consistent with no black hole. There is a clear disagreement on MBH between the JAM models/axisymmetric Schwarzschild models and these triaxial Schwarzschild models. To attempt to resolve these differences, we explored a wide range of tests for our triaxial models including: fitting only the inner higher S/N region, fitting sectors of the data, symmetrizing the kinematics, fitting only the radial velocity and velocity dispersion, adding various amounts of regularization, changing the total number of integrated orbits, and varying the input models and PSFs. In every test, the fitting results remained consistent. However, we note two interesting observations: (1) The green contours shown in Figure 3.7 show significant χ2 differences in the model in the region MBH < 2 × 105 M . At these masses, the BH sphere of influence is <0.00200 , which is well below the diffraction limit of our NIFS data. This is clearly unphysical as the data cannot possibly constrain BH masses in this low-mass regime (i.e., the green contours are closed well below the diffraction limit of our instrument). We note that if we ignore models with MBH . 3 × 105 the triaxial model results become fully consistent with the JAM/axisymmetric Schwarzschild models. (2) We calculated the χ2 value for each of the model kinematic moments and VRM S , independently. These values for two model BH masses are shown in Table 3.4, which shows the even kinematic moments and the VRM S favor a high-mass BH. However, the overall fit is clearly being driven by the odd velocity moments, especially the radial velocity. This is also unphysical, as the odd moments are supposed to provide virtually no constraints on the gravitational potential, as they have large freedom to vary, at fixed potential, to fit the data. As discussed in Section 3.3.4, comparing the VRM S profiles of the best-fit no black hole with a best-fit MBH ∼ 4 × 106 M shows a significantly better fit to the central pixels in the latter case. These observations lead us to speculate that the minimum χ2 at zero BH may be a numerical artifact and to favor the results from the JAM and axisymmetric Schwarzschild models of a detectable SMBH. MBH [M ] 104 4 × 106 Γ 0.75 0.65 χ2 Original 0.765 0.793 χ2 Vel Only 1.005 1.072 χ2 σ Only 0.875 0.819 χ2 h3 Only 0.791 0.805 χ2 h4 Only 0.498 0.494 χ2 VRM S 0.799 0.753 Table 3.4: Calculations of Triaxial Schwarzschild Model Reduced χ2 Independently 70 71 3.3.4 Summary of Dynamical Results In summary, we detect a central massive black hole with the JAM dynamical models where the best-fit MBH and Γ are 5.9 ± 1.1 × 106 M and 0.64 ± 0.02, respectively. With 6 the axisymmetric Schwarzschild models, we find the best-fit MBH = 2.5+1.8 −1.3 × 10 M and Γ = 0.67 ± 0.03 (1σ uncertainties). Finally, with the triaxial Schwarzschild models we find the results are consistent with no black hole and Γ = 0.75. However, the triaxial models show a small region that overlaps with the JAM/axisymmetric Schwarzschild models at the 3σ level. Despite the variations in the dynamical modeling results, all of the models provide better fits to the VRM S data with a BH mass between 2 − 6 × 106 M . This is particularly true in the central pixels, as shown in Figure 3.11. Here, we show a VRM S model comparison for all of the dynamical modeling techniques along the semimajor axis. The colored lines show the JAM (blue), axisymmetric Schwarzschild models (red), and triaxial Schwarzschild models (green) best-fit parameters for two hypothetical MBH , Γ combinations, shown as X’s in Figure 3.7. In this case, we show a MBH ∼ 4 × 106 M with Γ = 0.67 as solid lines and MBH ∼ 104 M with Γ = 0.74 as dashed lines. It is clear from this comparison plot that the high mass BH is favored in the VRM S profile for all of the dynamical modeling techniques, especially near the center where we expect the effects of a central massive BH are the most significant. The results of the dynamical modeling techniques show that we cannot constrain the lower limit of the mass of a central massive BH. However, the better fits to the central VRM S profiles provides evidence in favor of a detectable BH mass. Furthermore, the JAM and axisymmetric Schwarzschild models are nearly consistent at the 1σ level. By combining the 1σ confidence levels of the JAM and axisymmteric Schwarzschild models, we suggest the BH 6 mass in M59-UCD3 is MBH = 4.2+2.1 −1.7 ×10 M . This estimate is based on the average of the best-fit JAM and axisymmetric Schwarzschild models, where the uncertainties from each model were added in quadrature. We do the same for the best Γ value to find Γ = 0.65±0.04, which corresponds to an average M/LF 814W,dyn = 1.61 ± 0.10 and M/LV,dyn = 2.73 ± 0.17. Finally, we note that this study is, to our knowledge, the first time that a direct comparison has been made between these three dynamical modeling codes. As noted at the beginning of this section, in general, comparisons of JAM and Schwarzschild modeling have found consistent results. One interesting recent study by Leung et al. (2018) has compared both Schwarzschild and JAM models against circular velocities derived from molecular gas for 54 galaxies with CALIFA integral-field stellar kinematics. The study found that JAM 72 Figure 3.11: Black points show a rectangular aperture along the semi major axis. The solid lines represent a model ∼ 4 × 106 M BH with Γ = 0.67 for the JAM (blue), axisymmetric Schwarzschild model (red), and triaxial Schwarzschild model (green). The dashed line represents a ∼ 104 M BH with Γ = 0.74 using the same colored convention described above. Note that these MBH , Γ combinations are not the best-fit model from any of the dynamical models. This choice is arbitrary and is for visual comparison between a low and high mass BH only. and Schwarzschild recover consistent mass profiles, without evidence for systematic biases (their Fig. D1). However, it also found that JAM recovers more reliable circular velocities than the Schwarzschild models in the large-radii regime, where the gas velocities are more reliable (their Fig. 8). Although the study was not specific to SMBHs, it shows that the reduced generality of the JAM method, with respect to Schwarzschild’s, can lead to a more robust mass-profiles recovery from real observations. The lack of flexibility could be leading to a more robust result here too, especially if the kinematic data includes any outliers that are not well described by their error bars. Finally, we note that despite the disagreement in the BH mass, the overall the agreement between the models is quite good; apart from the triaxial Schwarzschild model χ2 minimum at zero BH mass, the confidence regions of 73 all three models overlap in both Γ and MBH . 3.4 Radio and X-ray Observations of UCDs An alternative method for inferring the presence of a SMBH in UCDs is via accretion, which produces X-ray and radio emission. X-ray emission alone is only suggestive, as low-mass X-ray binaries are common in dense stellar systems and can mimic the X-ray emission from a low-luminosity AGN. However, radio emission from low-mass X-ray binaries is not detectable at the distance of the Virgo Cluster and hence is a more secure indication of a SMBH. Here we consider the radio and X-ray emission from three massive UCDs around the Virgo galaxies M59 and M60: M59cO, M59-UCD3, and M60-UCD1, which all have dynamical evidence for SMBHs. We note that no deep radio data exist for the other UCDs with evidence of SMBHs. 3.4.1 Radio We obtained deep radio continuum data for M59 and M60 with the Karl G. Jansky Very Large Array (VLA) as part of program 15A-091 (PI: Strader) in February and March 2015. All data were taken in B configuration and with C band receivers in 3-bit mode, split into subbands centered at 5 and 7 GHz, each with 2 GHz of bandwidth. Four 1.75 hr long blocks were observed, and in each block observations alternated between the two targets, giving 3.5 hr of observations (2.6 hr on source) per galaxy. The data were flagged and calibrated in AIPS using standard methods, then imaged with Briggs robust weighting (Briggs 1995). The subband data were imaged separately (at central frequencies of 4.6 and 7.1 GHz after flagging) and together, at a mean frequency of 5.8 GHz. The beam in the combined images is 1.3300 × 1.1400 . M59-UCD3 is not detected in the individual subbands or in the combined image. The local rms noise in the region of M59-UCD3 is 2.6 µJy bm−1 . Hence, we set a 3σ upper limit of < 7.8µJy bm−1 (L < 1.27 × 1034 ) at a mean frequency of 5.8 GHz. M60-UCD1 is also undetected, with a local rms of 2.4 µJy bm−1 and a corresponding upper limit of < 7.2µJy bm−1 (L < 1.17 × 1034 ) at 5.8 GHz. In contrast, we do detect M59cO in the 4.6 GHz subband at a flux density of 10.8 ± 3.8 µJy bm−1 (L = 1.75 × 1034 ). It is not detected in the 7.1 GHz image. The UCD is detected in Gaia with a J2000 position of (R.A., Dec) = (12:41:55.334, +11:40:03.79), only 0.100 from the VLA position of the radio source in the 4.6 GHz image (R.A., Dec) = (12:41:55.331, +11:40:03.69). The astrometric match suggests that the radio emission, while faint, is indeed real and associated with M59cO. Here, the 74 luminosity is the flux density in µJy ×10−29 × 4πR2 × 5 × 109 . These are all given at 5 GHz (i.e., assuming a flux density slope of α = 0 (flat)). VLA mosaic images for these three UCDs is shown in Figure 3.12. 3.4.2 X-ray These UCDs have been studied in the X-rays using Chandra by several previous authors (Luo et al. 2013; Strader et al. 2013; Hou & Li 2016; Pandya et al. 2016), but we revisit this analysis to ensure consistency. All our results are consistent with these past studies. As noted in these previous studies, the X-ray emission from UCDs can be explained by low-mass X-ray binaries (LMXBs). In fact the number of X-ray sources falls short of expectations based on globular cluster X-ray sources, but SMBH emission cannot be excluded (Hou & Li 2016; Pandya et al. 2016). Given that 106−7 M SMBHs do seem to be present in UCDs, if these are accreting at the typical Eddington ratios seen for early-type galaxies (Lbol /Ledd ∼ 10−6 Ho 2009), we would expect the UCDs to have detectable X-ray sources of ∼1038 ergs/s. As discussed further in the next section, the radio emission from LMXBs is much lower than that expected for emission from SMBHs, and thus a detection of both X-ray and radio emission from a source would provide strong evidence for SMBH accretion. There are two separate observations of M59 (encompassing both M59-UCD3 and M59cO) and six observations of M60 that cover M60-UCD1; these are summarized in Table 3.5. We downloaded these observations from Chandra data archive and reprocessed them using CIAO 4.9 and CalDB 4.7.6. We used a 1.500 extraction radius around each source and measured the background in a larger nearby source-free area before normalizing the counts to the Figure 3.12: VLA images of a 1000 (800 pc) box around three massive UCDs with dynamical evidence for SMBHs: M59-UCD3, M59cO, and M60-UCD1. The red crosses mark the optical positions of the UCDs. M59cO has evidence for an associated radio source as discussed in the text, while M59-UCD3 and M60-UCD1 are not detected in the VLA images. 2068 8074 2068 8074 785 8182 8507 12976 12975 14328 all M59-UCD3 M59-UCD3 M59cO M59cO M60-UCD1 M60-UCD1 M60-UCD1 M60-UCD1 M60-UCD1 M60-UCD1 M60-UCD1 2000-04-20 2007-01-30 2007-02-01 2011-02-24 2011-08-08 2011-08-12 ... Epoch Date 51654.148669 54130.521317 54132.122880 55616.730009 55781.313337 55785.067169 ... Epoch Date (MJD) Effective Time (ksec) 24.8 5.3 24.8 5.3 38.1 52.4 17.5 101.0 84.9 14.0 307.9 ×10−15 Flux erg s−1 cm−2 ) 3.1+2.7 −1.7 ... < 0.5 < 4.4 1.8+1.3 −0.9 3.7+1.7 −1.3 3.6+3.7 −2.2 7.5+1.6 −1.4 2.8+1.2 −1.0 2.4+3.3 −1.9 3.3+0.8 −0.7 Luminosity (×1037 erg s−1 ) 10.1+8.8 −5.5 ... < 1.7 < 14.3 5.9+4.2 −2.9 12.1+5.5 −4.2 +12.1 11.7−7.2 24.4+5.2 −4.6 9.1+3.9 −3.3 +10.8 7.8−6.2 10.8+2.6 −2.3 All limits and uncertainties are at the 95% level and over the energy range 0.5–10 keV. Luminosities assume a distance of 16.5 Mpc. Obs. ID UCD ID Table 3.5: Chandra X-ray Constraints and Measurements 75 76 source extraction region size. We initially determined all counts in the 0.3–10 keV range for maximum sensitivity, but report results in the 0.5–10 keV range for appropriate comparison to the fundamental plane. For both galaxies, we fix NH = 2 × 1020 cm−2 (taking extinction from Schlafly & Finkbeiner 2011 and conversion from Bahramian et al. 2015). All spectral extractions were performed with CIAO task specextract, and spectral analysis was done using Xspec 12.9.1n (Arnaud 1996). We assumed Wilms et al. (2000) abundances and Verner et al. (1996) absorption cross-sections. M59cO is not detected in the 2001 or 2008 observations. Assuming a power-law with Γ = 1.5, in the 2001 data we find a 95% upper 0.5–10 keV unabsorbed flux limit of < 5.3 × 10−16 erg s−1 cm−2 , equivalent to LX < 1.7 × 1037 erg s−1 . The shorter 2008 data are less constraining and give a limit of LX < 1.4 × 1038 erg s−1 using the same assumptions. M59-UCD3 is detected at > 2σ in the 2001 Chandra data, with a 0.5–10 keV unabsorbed −15 erg s−1 cm−2 , equivalent to 1.0+0.9 × 1038 erg s−1 (uncertainties are flux of 3.1+2.7 −0.6 −1.7 × 10 at the 95% level). Unsurprisingly it is not detected in the factor of ∼ 5 shorter 2008 data. In addition, it is located near a chip gap in the 2008 observations, which makes it difficult to determine a valid upper flux limit. Here, we have assumed Gehrels statistics for all of the upper limits (Gehrels 1986). M60-UCD1 is detected in all six observations. The total merged dataset, representing 308 ksec of Chandra data, is deep enough to allow spectral fitting. After binning to 20 counts per bin, we fit the spectrum to a power-law in XSPEC using cstat, a modified +0.2 , version of the Cash statistic (Cash 1979).6 The best-fitting power law index is Γ = 1.8−0.3 consistent with the Γ = 1.5 value assumed. Hence, for consistency, we assume Γ = 1.5 for all the flux measurements for M60-UCD1. The individual unabsorbed 0.5–10 keV fluxes for M60-UCD1 range over 1.8–7.5 × 10−15 erg s−1 cm−2 (LX = 0.6–2.4×1038 erg s−1 , depending on the epoch. The average flux is −15 erg s−1 cm−2 (L = 1.1+0.3 × 1038 erg s−1 ). The individual and merged (3.3+0.8 X −0.7 × 10 −0.2 fluxes are listed in Table 3.5. There is compelling evidence for X-ray variability of M60-UCD1, but only at a single epoch: five of the six epochs are consistent with the mean flux, while one (Obs ID 12976) is ∼ 6σ higher compared to the mean flux. Due to the shorter exposure times and smaller number of epochs for M59-UCD3 and M59cO, we have no useful constraints on X-ray variability for these other sources. 6 https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html 77 3.4.3 Fundamental Plane of BH Accretion We can combine X-ray and radio detections and nondetections described above with the dynamical black hole mass estimates to see if these observations are consistent with the fundamental plane of BH accretion (Merloni et al. 2003; Falcke et al. 2004). We note that significant variability is seen in M60-UCD1, and long-term variability of low-luminosity AGN seems to be common (e.g., Maoz et al. 2005; Hernández-Garcı́a et al. 2014), and in general UCD X-ray sources appear to be variable (Hou & Li 2016; Pandya et al. 2016). Because our radio and X-ray observations are not contemperaneous, this variability adds to any intrinsic scatter present. We use the fundamental plane of Plotkin et al. (2012): log(LX ) = (1.45 ± 0.04) log(LR ) − (0.88 ± 0.06) log(MBH ) − 6.07 ± 1.10 (3.1) and plot the combinations of detections and upper limits for the three UCDs in Figure 3.13. We find that the radio upper limits in M59-UCD3 and especially M60-UCD1 fall well below the radio luminosities expected for objects lying on the fundamental plane. Similarly, the X-ray nondetection of M59cO is below the expectation based on its radio luminosity. However, these nondetections do not provide strong constraints on whether an accreting BH is present. This is because of the order of magnitude scatter in the radio luminosities relative to the fundamental plane observed in similar low-Eddington systems with known BH masses (Gültekin et al. 2009), as well as the lack of simultaneous radio and X-ray observations. Future, simultaneous detections of X-ray and radio emission in UCD BHs could provide important confirming evidence for SMBHs in these systems. Finally, we note that from our NIFS observations, we can also constrain the presence of hot dust emission in these systems. This hot dust is typically thought to be an AGN accretion signature and is found to be quite luminous in both quiescent and actively accreting BHs (Seth et al. 2010; Seth 2010; Burtscher et al. 2015). A clear correlation is seen between this emission and the X-ray and midinfrared emission, although this correlation seems to depend significantly on AGN type (Burtscher et al. 2015). Whether this correlation continues to lower luminosity AGN like those observed here is not yet clear. 3.5 Discussion and Conclusions In this paper, we have presented the results of three separate dynamical models on the most massive UCD, M59-UCD3, and discussed the radio and X-ray observations of UCDs as a way to infer the presence of SMBHs. Detections of SMBHs in UCDs provide evidence 78 Figure 3.13: The fundamental plane for UCDs with known BHs from (Plotkin et al. 2012). The radio upper limits in M60-UCD1 and M59-UCD3 and the X-ray upper limits in M59cO (arrows) are all below the predictions from the fundamental plane, but these differences are small relative to the scatter seen for similar quiescent black holes with dynamical BH mass determinations shown in gray from Gültekin et al. (2009). that they are the tidally stripped remnants of once more massive galaxies (Seth et al. 2014; Ahn et al. 2017; Afanasiev et al. 2018). Furthermore, the effects of the central massive BH can explain the elevated dynamical-to-stellar mass ratios detected in the most massive UCDs (Mieske et al. 2013). Nondetection of BHs can suggest either that the object is not a stripped nucleus (and instead is a massive star cluster), or that the nucleus stripped lacked detectably massive central BHs (Voggel et al. 2018) . 3.5.1 Summary of Main Results For our analysis, we combined adaptive optics assisted Gemini/NIFS kinematic data with high resolution HST imaging. The Gemini/NIFS data were used to determine the full 79 LOSVD, which includes the radial velocity, velocity dispersion, skewness, and kurtosis. We found the integrated (r < 0.7500 ) barycentric radial velocity to be V = 434.5 ± 0.6 km s−1 and velocity dispersion σ = 65.7 ± 0.6 km s−1 . The HST images were used to construct a mass density, and luminosity profile by fitting a PSF convolved triple component Sérsic profile to the data. These models were used to calculate the total luminosity (within the central 300 ), which we found to be LF 814W = 1.19 × 108 L and an effective radius of 0.3400 (27 pc). The model fits suggest the outer component of the UCD is somewhat bluer than the central components, and we modeled this stellar population variation using SSP models. We combined these mass and luminosity density profiles with the kinematic measurements to test for the presence of a central massive BH using three dynamical modeling techniques. We summarize the results of this modeling in (Section 3.3.4); our final conclusion 6 is that M59-UCD3 appears to host a BH with a mass of MBH = 4.2+2.1 −1.7 ×10 M . We derive a best-fit Γ from the JAM and axisymmetric Schwarzschild models of 0.65 ± 0.04, which corresponds to an M/LF 814W,dyn = 1.61 ± 0.10 and M/LV,dyn = 2.73 ± 0.17. Therefore, the total dynamical mass is Mdyn = 1.9 ± 0.1 × 108 M . 3.5.2 Implications of a Central Massive BH We can compare our best-fit models with a BH to those without to look at how large a change is caused in the stellar M/L. From our JAM models, the best-fit zero mass BH has a Γ of 0.74 (M/LV,dyn = 3.11). This represents our equivalent to the inferred mass from integrated dispersions used to determine the ratio of dynamical to population values (e.g., Mieske et al. 2013). Therefore, inclusion of a BH in the model reduces the M/L by ∼12%. By contrast, the M/Ls drop by > 40% in all of the other UCDs with central massive BH detections (Seth et al. 2014; Ahn et al. 2017). We show the effect on Γ of including the BH in M59-UCD3 in Fig. 3.14, which, by our determination, does not have an elevated global dynamical M/L even if there is no BH present. We note that Liu et al. (2015) estimated M/LV,dyn to be 4.9±0.5 based on an integrated dispersion of 77.8±1.6 km s−1 and ref f = 25 pc, shown as an open star in Figure 3.14; this estimate is significantly higher than our estimates with or without a BH. Voggel et al. (2018) showed that measurements of the dynamical-to-stellar mass ratio are easily overestimated using lower spatial resolution data, which demonstrates the importance of high spatial resolution integral field unit (IFU) data for determining the dynamical M/L. This is likely the reason for the discrapancy between M/Ls in Liu et al. (2015) and our derived value. In M60-UCD1, VUCD3, and M59cO, the central BH constituted 10–20% of the total mass of the system (Seth et al. 2014; Ahn et al. 2017). M59-UCD3 is quite different, with a 80 Figure 3.14: Dynamical-to-stellar mass ratio Γ vs. total dynamical mass. Grey points represent GCs and UCDs, where mass estimates are based on integrated velocity dispersions assuming mass traces light from Mieske et al. (2013) and references therin (with the exceptions of UCD3, UCD 320, UCD 330 and M59-UCD3; Frank et al. 2011; Afanasiev et al. 2018; Voggel et al. 2018). Here, stars represent the seven UCDs and two GCs for which AO-assisted data has or will be analyzed. The colored stars represent updated stellar mass measurements after accounting for the central massive BH. Arrows illustrate the effect of the BH. The open grey star represents the M59-UCD3 estimate from Liu et al. (2015). In the case of UCD 320 and UCD 330, Voggel et al. (2018) found their initial dynamical-to-stellar mass ratios were overestimated and did not detect a central massive BH in either object. 81 central BH of only ∼2% of the total mass. A small, ∼4% mass BH has also been found in Fornax UCD3 (Afanasiev et al. 2018) but due to much lower S/N data (with similar spatial resolution); this BH was only detectable after assuming isotropy and fixing the stellar M/L to expected stellar population values. In the context of the other UCDs with high BH mass fractions, it is not surprising that the BH mass in M59-UCD3 is less well determined. Using the conventional definition of the BH sphere of influence (rinfl = GMBH /σ 2 ), we find rinf l = 0.03 − 0.0700 , assuming σ = 65.7 km s−1 , for M59-UCD3. Here, the BH sphere of influence range is calculated assuming the best-fit MBH for the axisymmetric Schwarzschild models and JAM models, respectively. Therefore, the BH sphere of influence for M59-UCD3 is approximately an order of magnitude less than what was found for VUCD3, M59cO, and M60-UCD1 (Seth et al. 2014; Ahn et al. 2017). Furthermore, since our adaptive optics PSF has a core FWHM of 0.16500 (and the diffraction limit is 0.0700 ), this likely explains why we can’t constrain the lower mass limit of the central BH with all three dynamical models. However, the factor of ∼10 range in BH masses in UCDs with similar stellar mass is not particularly surprising given the comparable range of BH masses seen in lower mass (especially spiral) galaxies at a fixed dispersion or stellar mass (e.g., Kormendy & Ho 2013; Greene et al. 2016; Nguyen et al. 2017). 3.5.3 Progenitor Galaxy Properties The combination of M59-UCD3’s high luminosity and apparent BH strongly suggests it is a tidally stripped remnant of a once more massive galaxy. We can try to estimate the progenitor mass in several ways. Assuming that UCDs follow the same MBH vs. Bulge Luminosity relation as galaxies (e.g., McConnell & Ma 2013; Mieske et al. 2013), we obtain a progenitor bulge mass of ∼ 109 M . This corresponds to the total galaxy mass if we assume an early-type galaxy; in this case, the total galaxy mass would only be ∼3× that of the current UCD, which would make the galaxy remarkably compact, similar to M32’s current size and mass. We can also estimate the mass of the NSC by assuming it is the central two Sérsic components of the UCD (given that a King + Sérsic model provides an equally good fit; Pfeffer & Baumgardt 2013); this component has a mass of MN SC ≡ 28 × 107 M and an effective radius of ref f = 21 pc (0.2600 ). Galaxies with similar NSCs in the Georgiev et al. (2016) sample have stellar masses between 4×109 and 1011 M . A similar total stellar mass range is found for galaxies with BHs between 106 and 107 M in Reines & Volonteri (2015), while the M − σ relation implies a galaxy dispersion of ∼100 km s−1 (e.g., Kormendy & Ho 2013). Overall, it appears likely that the original mass of M59-UCD3’s progenitor was 82 of order 1010 M . M59 has two UCDs whose progenitor galaxies masses have been estimated to be ∼ 109 − 1010 M . Therefore, we would expect the stars of these disrupted galaxies to be deposited in the outer halo of M59. Using the Sérsic fits from Kormendy et al. (2009), we calculated the total luminosity outside the projected radius of M59-UCD3 and M59cO (both ∼10 kpc), which was found to be L(> 12700 ) ∼ 1.2 × 1010 L . Here we have assumed a constant axis ratio and position angle of 0.7 and 164◦ , respectively. We also note that, Liu et al. (2015) reported a plume with similar luminosity to M59-UCD3 itself, which may be a tidal feature associated with it. Unless M59-UCD3 was remarkably compact before stripping, it appears likely that this plume, if associated with M59-UCD3, represents just a small fraction of the total amount of mass stripped off the galaxy. 3.5.4 Further Evidence and Future Prospects for SMBHs in UCDs In this paper, we have also presented deep radio observations of three UCDs around M59 and M60, which showed a detection in only one case (around M59cO). The radio upper limits in M59-UCD3 and M60-UCD1 are lower than expected for AGN emission based on X-ray detections in these systems, although given the scatter in the fundamental plane, these measurements do not necessarily argue against an SMBH origin for the X-ray emission. M59-UCD3 represents the fifth UCD known to host a SMBH. All UCDs with measured BH masses have been shown to have near solar metallicities and α-enhancement in the range of [Mg/Fe]∼ 0.1-0.5 (Chilingarian & Mamon 2008; Francis et al. 2012; Sandoval et al. 2015; Janz et al. 2016). Both UCDs consistent with no central BH have been shown to have subsolar metallicities and solar [α/Fe] ratios. This could indicate that solar–supersolar metallicity UCDs with alpha enhancement may serve as a secondary indicator of a central massive BH. With current telescopes and instrumentation our dynamical modeling effort has, thus far, been limited to either the brightest UCDs at the distance of Virgo/Fornax, or some fainter UCDs at the distance of nearby large galaxies such as Centaurus A. This is due to a combination of the need to resolve the BH sphere of influence, as well as limits on the faintest sources observable with adaptive optics corrections. As the next generation telescopes come online, such as the Extremely Large Telescope (ELT), we will be able to significantly increase the number of UCDs for which we have the capability to run dynamical models. For example, the new High Angular Resolution Monolithic Optical and Near-infrared Integral Field Spectrograph (HARMONI), which is to be mounted on ELT, 83 has an estimated adaptive optics corrected angular resolution of 0.01 − 0.0300 (Cunningham et al. 2008). Assuming two hypothetical BHs with MBH = 107 M and MBH = 106 M and corresponding integrated dispersions of ∼50 km s−1 and ∼30 km s−1 , the BH sphere of influence would be ∼16 pc and ∼5 pc, respectively. If we require the BH sphere of influence to be at least 0.0500 in size for us to resolve it with ELT/HARMONI, then we could theoretically resolve these BHs out to 19 Mpc (the distance of Fornax/Virgo) for the 106 M BH and an astonishing 66 Mpc (the distance of the Perseus cluster and many others) for the 107 M BH. CHAPTER 4 CONCLUSIONS AND FUTURE PROSPECTS 4.1 Conclusions I’ve presented the detection of SMBHs in three UCDs: VUCD3, M59cO, and M59UCD3, which brings the total number of UCDs with dynamical evidence for central massive BHs up to five (Seth et al. 2014; Ahn et al. 2017, 2018; Afanasiev et al. 2018). These detections, combined with their structural information, suggests that all of these UCDs are likely the tidally stripped remnants of 109 −1010 M galaxies. This information has provided us with clues about the role UCDs play in overall galaxy formation and evolution, especially in cluster environments where UCDs are common. However, with only five UCDs that have dynamical BH detections, we are limited by small number statistics. Furthermore, most UCDs that we are capable of modeling dynamically have been examined. This is in large part due to current technology limits, which restrict dynamical modeling of these compact objects to the brightest UCDs at the distance of Virgo/Fornax. 4.2 Future: Galaxy Cluster Assembly With the upcoming new generation of telescopes, such as the 40 m ELT, we will be able to find SMBHs in UCDs out to three times the distance of our current limitations. The sheer number of UCDs within the volume encompassed by the ELT capabilities will provide a statistically significant sample of how many UCDs host SMBHs in local clusters. Knowing that these UCDs with SMBHs are the former nuclei of galaxies, we can use UCDs to trace back how many galaxies have been stripped into galaxy clusters. Therefore, we hope to use UCDs as a new type of direct tracer to understanding galaxy cluster assembly history in the near future. 85 4.3 Future: Scaling Relations and UCD Formation In the past, scaling relations have been used to look for patterns and trends among galaxy properties, which conceal underlying physics. For example, as discussed in Chapter 1, the mass of the central BH correlates with the host galaxy velocity dispersion. This property suggests some sort of feedback mechanism that governs the growth of SMBHs and the growth of host galaxy bulges, or coevolution. With a statistically significant sample of the BH occupation fraction in UCDs, we can begin to understand where these compact objects fit into the overall galaxy formation and evolution scenarios. Our current measurements imply that UCDs most closely resemble nuclear star clusters at the centers of most massive galaxies (Mieske et al. 2013; Georgiev et al. 2016), albeit with some remaining extended component. Since massive UCDs are likely the tidally stripped remnants of 109 − 1010 M galaxies, it is not surprising that the nuclear star cluster could remain intact, which is confirmed by tidal stripping simulations (Pfeffer et al. 2014, 2016). However, as stated above, as the number of UCDs with dynamical BH mass estimates increases, we will further constrain the role UCDs play in galaxy formation and evolution. REFERENCES Abbott, B. P., et al. 2016a, PhRvL, 116, 061102 —. 2016b, ApJL, 833, L1 Afanasiev, A. V., et al. 2018, MNRAS, 477, 4856 Ahn, C. P., et al. 2017, ApJ, 839, 72 —. 2018, ApJ, 858, 102 Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco, CA: ASP), 17 Bahramian, A., Heinke, C. O., Degenaar, N., et al. 2015, MNRAS, 452, 3475 Baumgardt, H. 2017, MNRAS, 464, 2174 Begelman, M. C., Volonteri, M., & Rees, M. J. 2006, MNRAS, 370, 289 Bekki, K., Couch, W. J., & Drinkwater, M. J. 2001, ApJL, 552, L105 Bekki, K., Couch, W. J., Drinkwater, M. J., & Shioya, Y. 2003, MNRAS, 344, 399 Boccas, M., Rigaut, F., Bec, M., et al. 2006, in Proc. SPIE, Vol. 6272, Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, 62723L Briggs, D. S. 1995, PhD thesis, The New Mexico Institute of Mining and Technology Brodie, J. P., Romanowsky, A. J., Strader, J., & Forbes, D. A. 2011, AJ, 142, 199 Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 Burtscher, L., Orban de Xivry, G., Davies, R. I., et al. 2015, A&A, 578, A47 Cappellari, M. 2002, MNRAS, 333, 400 —. 2008, MNRAS, 390, 71 —. 2017, MNRAS, 466, 798 Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345 Cappellari, M., & Emsellem, E. 2004, PASP, 116, 138 Cappellari, M., Neumayer, N., Reunanen, J., et al. 2009, MNRAS, 394, 660 Cappellari, M., Verolme, E. K., van der Marel, R. P., et al. 2002, ApJ, 578, 787 Cappellari, M., Bacon, R., Bureau, M., et al. 2006, MNRAS, 366, 1126 87 Cappellari, M., McDermid, R. M., Bacon, R., et al. 2010, in American Institute of Physics Conference Series, Vol. 1240, Hunting for the Dark: The Hidden Side of Galaxy Formation, ed. V. P. Debattista & C. C. Popescu (Melville, NY: AIP), 211–214 Cash, W. 1979, ApJ, 228, 939 Celotti, A., Miller, J. C., & Sciama, D. W. 1999, Classical and Quantum Gravity, 16, A3 Chilingarian, I. V., & Mamon, G. A. 2008, MNRAS, 385, L83 Colless, M., et al. 2001, MNRAS, 328, 1039 Conroy, C., & van Dokkum, P. G. 2012, ApJ, 760, 71 Côté, P., et al. 2006, ApJS, 165, 57 Cunningham, C., Evans, C., Monnet, G., & Le Louarn, M. 2008, in Proc. SPIE, Vol. 6986, 69860K Da Rocha, C., Mieske, S., Georgiev, I. Y., et al. 2011, A&A, 525, A86 Dabringhausen, J., Fellhauer, M., & Kroupa, P. 2010, MNRAS, 403, 1054 Dabringhausen, J., Hilker, M., & Kroupa, P. 2008, MNRAS, 386, 864 Dabringhausen, J., Kroupa, P., & Baumgardt, H. 2009, MNRAS, 394, 1529 Dabringhausen, J., Kroupa, P., Pflamm-Altenburg, J., & Mieske, S. 2012, ApJ, 747, 72 den Brok, M., Seth, A. C., Barth, A. J., et al. 2015, ApJ, 809, 101 Drehmer, D. A., Storchi-Bergmann, T., Ferrari, F., Cappellari, M., & Riffel, R. A. 2015, MNRAS, 450, 128 Drinkwater, M. J., Gregg, M. D., Hilker, M., et al. 2003, Natur, 423, 519 Drinkwater, M. J., Jones, J. B., Gregg, M. D., & Phillipps, S. 2000, PASA, 17, 227 Eggen, O. J., Lynden-Bell, D., & Sandage, A. R. 1962, ApJ, 136, 748 Emsellem, E., Monnet, G., & Bacon, R. 1994, A&A, 285 Evstigneeva, E. A., Gregg, M. D., Drinkwater, M. J., & Hilker, M. 2007, AJ, 133, 1722 Evstigneeva, E. A., Drinkwater, M. J., Peng, C. Y., et al. 2008, AJ, 136, 461 Falcke, H., Körding, E., & Markoff, S. 2004, A&A, 414, 895 Feldmeier, A., Neumayer, N., Seth, A., et al. 2014, A&A, 570, A2 Feldmeier-Krause, A., Zhu, L., Neumayer, N., et al. 2017, MNRAS, 466, 4040 Fellhauer, M., & Kroupa, P. 2002, MNRAS, 330, 642 —. 2005, MNRAS, 359, 223 Ferrarese, L., & Ford, H. 2005, Space Science Reviews, 116, 523 Ferrarese, L., & Merritt, D. 2000, ApJL, 539, L9 88 Firth, P., Evstigneeva, E. A., & Drinkwater, M. J. 2009, MNRAS, 394, 1801 Forbes, D. A., Norris, M. A., Strader, J., et al. 2014, MNRAS, 444, 2993 Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 Francis, K. J., Drinkwater, M. J., Chilingarian, I. V., Bolt, A. M., & Firth, P. 2012, MNRAS, 425, 325 Frank, M. J., Hilker, M., Mieske, S., et al. 2011, MNRAS, 414, L70 Gebhardt, K., Rich, R. M., & Ho, L. C. 2005, ApJ, 634, 1093 Gebhardt, K., et al. 2000, ApJL, 539, L13 Gebhardt, K., Richstone, D., Tremaine, S., et al. 2003, ApJ, 583, 92 Geha, M., Guhathakurta, P., & van der Marel, R. P. 2003, AJ, 126, 1794 Gehrels, N. 1986, ApJ, 303, 336 Genzel, R., Eckart, A., Ott, T., & Eisenhauer, F. 1997, MNRAS, 291, 219 Georgiev, I. Y., Böker, T., Leigh, N., Lützgendorf, N., & Neumayer, N. 2016, MNRAS, 457, 2122 Ghez, A. M., Klein, B. L., Morris, M., & Becklin, E. E. 1998, ApJ, 509, 678 Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, Vol. 5, No. 1, p. 65-80, 2010, 5, 65 Gould, N. I. M., Orban, D., & Toint, P. L. 2003, ACM Trans. Math. Softw., 29, 353 Graham, A. W. 2016, in Astrophysics and Space Science Library, Vol. 418, Galactic Bulges, ed. E. Laurikainen, R. Peletier, & D. Gadotti, 263 Graham, A. W., & Spitler, L. R. 2009, MNRAS, 397, 2148 Greene, J. E., Seth, A., Kim, M., et al. 2016, ApJL, 826, L32 Gültekin, K., Cackett, E. M., Miller, J. M., et al. 2009, ApJ, 706, 404 Haşegan, M., Jordán, A., Côté, P., et al. 2005, ApJ, 627, 203 Haggard, D., Cool, A. M., Heinke, C. O., et al. 2013, ApJL, 773, L31 Hartmann, M., Debattista, V. P., Seth, A., Cappellari, M., & Quinn, T. R. 2011, MNRAS, 418, 2697 Hernández-Garcı́a, L., González-Martı́n, O., Masegosa, J., & Márquez, I. 2014, A&A, 569, A26 Herriot, G., Morris, S., Anthony, A., et al. 2000, in Proc. SPIE, Vol. 4007, Adaptive Optical Systems Technology, ed. P. L. Wizinowich, 115–125 Hilker, M., Infante, L., Vieira, G., Kissler-Patig, M., & Richtler, T. 1999, A&A Supp., 134, 75 89 Ho, L. C. 2009, ApJ, 699, 626 Hou, M., & Li, Z. 2016, ApJ, 819, 164 Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6 Ibata, R., Bellazzini, M., Chapman, S. C., et al. 2009, ApJL, 699, L169 Janz, J., Forbes, D. A., Norris, M. A., et al. 2015, MNRAS, 449, 1716 Janz, J., Norris, M. A., Forbes, D. A., et al. 2016, MNRAS, 456, 617 Jarrett, T. H., Chester, T., Cutri, R., Schneider, S. E., & Huchra, J. P. 2003, AJ, 125, 525 Jeans, J. H. 1915, MNRAS, 76, 70 Kamann, S., Husser, T.-O., Dreizler, S., et al. 2018, MNRAS, 473, 5591 Kimmig, B., Seth, A., Ivans, I. I., et al. 2015, AJ, 149, 53 Kissler-Patig, M., Jordán, A., & Bastian, N. 2006, A&A, 448, 1031 Kızıltan, B., Baumgardt, H., & Loeb, A. 2017, Natur, 542, 203 Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 82 Koleva, M., Prugniel, P., de Rijcke, S., & Zeilinger, W. W. 2011, MNRAS, 417, 1643 Kormendy, J., Fisher, D. B., Cornell, M. E., & Bender, R. 2009, ApJS, 182, 216 Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511 Krajnović, D., Cappellari, M., de Zeeuw, P. T., & Copin, Y. 2006, MNRAS, 366, 787 Krajnović, D., McDermid, R. M., Cappellari, M., & Davies, R. L. 2009, MNRAS, 399, 1839 Krajnović, D., Cappellari, M., McDermid, R. M., et al. 2018, MNRAS, arXiv:1803.08055 Lawson, C. L., & Hanson, R. J. 1974, Solving least squares problems (Englewood Cliffs, NJ: Prentice-Hall) Leung, G. Y. C., Leaman, doi:10.1093/mnras/sty288 R., van de Ven, G., et al. 2018, MNRAS, Liu, C., Peng, E. W., Toloba, E., et al. 2015, ApJL, 812, L2 Luo, B., Fabbiano, G., Strader, J., et al. 2013, ApJS, 204, 14 Maoz, D., Nagar, N. M., Falcke, H., & Wilson, A. S. 2005, ApJ, 625, 699 Maraston, C. 2005, MNRAS, 362, 799 Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS, 351, 169 Markwardt, C. B. 2009, in Astronomical Society of the Pacific Conference Series, Vol. 411, Astronomical Data Analysis Software and Systems XVIII, ed. D. A. Bohlender, D. Durand, & P. Dowler (San Francisco, CA: ASP), 251 McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184 90 McGregor, P. J., Hart, J., Conroy, P. G., et al. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, ed. M. Iye & A. F. M. Moorwood, 1581–1591 Mei, S., Blakeslee, J. P., Côté, P., et al. 2007, ApJ, 655, 144 Merloni, A., Heinz, S., & di Matteo, T. 2003, MNRAS, 345, 1057 Meylan, G., Sarajedini, A., Jablonka, P., et al. 2001, AJ, 122, 830 Mieske, S., Frank, M. J., Baumgardt, H., et al. 2013, A&A, 558, A14 Mieske, S., Hilker, M., & Infante, L. 2002, A&A, 383, 823 Mieske, S., & Kroupa, P. 2008, ApJ, 677, 276 Mieske, S., Hilker, M., Jordán, A., et al. 2008, A&A, 487, 921 Miller-Jones, J. C. A., Wrobel, J. M., Sivakoff, G. R., et al. 2012, ApJL, 755, L1 Misgeld, I., & Hilker, M. 2011, MNRAS, 414, 3699 Moffat, A. F. J. 1969, A&A, 3, 455 Moore, B. 1994, Natur, 370, 629 Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJL, 524, L19 Murray, N. 2009, ApJ, 691, 946 Nguyen, D. D., Seth, A. C., den Brok, M., et al. 2016, ArXiv e-prints, arXiv:1610.09385 Nguyen, D. D., Seth, A. C., Neumayer, N., et al. 2017, ArXiv e-prints, arXiv:1711.04314 Norris, J. E., Freeman, K. C., Mayor, M., & Seitzer, P. 1997, ApJL, 487, L187 Norris, M. A., Escudero, C. G., Faifer, F. R., et al. 2015, MNRAS, 451, 3615 Norris, M. A., & Kannappan, S. J. 2011, MNRAS, 414, 739 Norris, M. A., Kannappan, S. J., Forbes, D. A., et al. 2014, MNRAS, 443, 1151 Nowak, N., Saglia, R. P., Thomas, J., et al. 2008, MNRAS, 391, 1629 —. 2007, MNRAS, 379, 909 Noyola, E., Gebhardt, K., Kissler-Patig, M., et al. 2010, ApJL, 719, L60 Pandya, V., Mulchaey, J., & Greene, J. E. 2016, ApJ, 819, 162 Paudel, S., Lisker, T., & Kuntschner, H. 2011, MNRAS, 413, 1764 Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266 Pfeffer, J., & Baumgardt, H. 2013, MNRAS, 433, 1997 Pfeffer, J., Griffen, B. F., Baumgardt, H., & Hilker, M. 2014, MNRAS, 444, 3670 Pfeffer, J., Hilker, M., Baumgardt, H., & Griffen, B. F. 2016, MNRAS, 458, 2492 91 Phillipps, S., Young, A. J., Drinkwater, M. J., Gregg, M. D., & Karick, A. 2013, MNRAS, 433, 1444 Plotkin, R. M., Markoff, S., Kelly, B. C., Körding, E., & Anderson, S. F. 2012, MNRAS, 419, 267 Reines, A. E., Greene, J. E., & Geha, M. 2013, ApJ, 775, 116 Reines, A. E., & Volonteri, M. 2015, ApJ, 813, 82 Richstone, D. O., & Tremaine, S. 1988, ApJ, 327, 82 Rix, H.-W., de Zeeuw, P. T., Cretton, N., van der Marel, R. P., & Carollo, C. M. 1997, ApJ, 488, 702 Saglia, R. P., Opitsch, M., Erwin, P., et al. 2016, ApJ, 818, 47 Sandoval, M. A., Vo, R. P., Romanowsky, A. J., et al. 2015, ApJL, 808, L32 Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103 Schödel, R., Merritt, D., & Eckart, A. 2009, A&A, 502, 91 Schwarzschild, M. 1979, ApJ, 232, 236 Scott, N., & Graham, A. W. 2013, ApJ, 763, 76 Searle, L., & Zinn, R. 1978, ApJ, 225, 357 Seth, A. C. 2010, ApJ, 725, 670 Seth, A. C., Cappellari, M., Neumayer, N., et al. 2010, ApJ, 714, 713 Seth, A. C., van den Bosch, R., Mieske, S., et al. 2014, Natur, 513, 398 Shanahan, R. L., & Gieles, M. 2015, MNRAS, 448, L94 Shapiro, K. L., Cappellari, M., de Zeeuw, T., et al. 2006, MNRAS, 370, 559 Soltan, A. 1982, MNRAS, 200, 115 Springel, V., et al. 2005, Natur, 435, 629 Strader, J., Caldwell, N., & Seth, A. C. 2011, AJ, 142, 8 Strader, J., Seth, A. C., Forbes, D. A., et al. 2013, ApJL, 775, L6 Taylor, J. H., Hulse, R. A., Fowler, L. A., Gullahorn, G. E., & Rankin, J. M. 1976, ApJL, 206, L53 Taylor, J. H., & Weisberg, J. M. 1989, ApJ, 345, 434 Taylor, M. A., Puzia, T. H., Harris, G. L., et al. 2010, ApJ, 712, 1191 Thomas, J., Saglia, R. P., Bender, R., et al. 2004, MNRAS, 353, 391 Tollerud, E. J., Bullock, J. S., Graves, G. J., & Wolf, J. 2011, ApJ, 726, 108 Toomre, A., & Toomre, J. 1972, ApJ, 178, 623 92 Tremaine, S., et al. 2002, ApJ, 574, 740 Valluri, M., Ferrarese, L., Merritt, D., & Joseph, C. L. 2005, ApJ, 628, 137 van den Bosch, R., de Zeeuw, T., Gebhardt, K., Noyola, E., & van de Ven, G. 2006, ApJ, 641, 852 van den Bosch, R. C. E. 2016, ApJ, 831, 134 van den Bosch, R. C. E., & de Zeeuw, P. T. 2010, MNRAS, 401, 1770 van den Bosch, R. C. E., van de Ven, G., Verolme, E. K., Cappellari, M., & de Zeeuw, P. T. 2008, MNRAS, 385, 647 van der Marel, R. P., & Anderson, J. 2010, ApJ, 710, 1063 van der Marel, R. P., Cretton, N., de Zeeuw, P. T., & Rix, H.-W. 1998, ApJ, 493, 613 van Dokkum, P. G., & Conroy, C. 2010, Natur, 468, 940 Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487 Verolme, E. K., Cappellari, M., Copin, Y., et al. 2002, MNRAS, 335, 517 Villaume, A., Brodie, J., Conroy, C., Romanowsky, A. J., & van Dokkum, P. 2017, ApJL, 850, L14 Vogelsberger, M., et al. 2014, Natur, 509, 177 Voggel, K. T., Seth, A. C., Neumayer, N., et al. 2018, ArXiv e-prints, arXiv:1803.09750 Volonteri, M. 2010, A&AR, 18, 279 Wallace, L., & Hinkle, K. 1996, ApJS, 107, 312 Whalen, D. J., & Fryer, C. L. 2012, ApJL, 756, L19 Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914