The Structure and Evolution of Interacting Binary Galaxies Citation Borne, Kirk Daniel (1983) The Structure and Evolution of Interacting Binary Galaxies. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/8PH6-6R18. https://resolver.caltech.edu/CaltechETD:etd-09082008-132748 Abstract A numerical code has been constructed for the study of the evolution of interacting binary galaxies. This "multiple three-body" algorithm (MTBA) essentially involves N concurrent three-body integrations; it is the sum of the interactions of the perturber with each of the N particles comprising the primary "galaxy" that determines the binary orbital evolution. MTBA incorporates a violent relaxation phase that allows the particles to redistribute themselves in the gravitational field of the perturber prior to the full binary orbital evolution calculation. This redistribution is important for systems with an initially strong tidal potential; their predicted merger times are 50-100% larger than previously estimated. Tabulated merger times for both circular and elliptical orbits demonstrate that typical close binary galaxies will merge in about twice their initial binary orbital periods, roughly 10 9 years. This timescale depends strongly on the binary separation, weakly on its mass ratio, and less on initial galaxy rotation than has heretofore been suggested. A specific interacting binary simulation is described in detail in the first paper. Surface brightness maps, structural parameters, and rotation properties are described at various times during the evolution, followed by a discussion of the properties of the merger remnant. Many of the results reported here are consistent with those obtained from the larger, more expensive N-body simulations. For the second paper, MTBA is altered so that each "galaxy" is represented by a configuration of test particles. Simulated tidal interactions between these "galaxies" are matched to the measured properties of real binaries in an effort to determine the physical state of each observed pair. The catalog of isolated galaxy pairs prepared by Karachentsev has been culled for its E-E constituents. CCD pictures and long-slit spectroscopic observations have been obtained at KPNO for the most photometrically distorted of these systems. Such pairs very likely represent strongly interacting, physically associated, binary elliptical galaxies. Radial variations of rotation velocity and velocity dispersion are extracted from the spectroscopic data for each of the two galaxies in a given pair. The combined spectroscopic and photometric data for that binary can be reproduced by an appropriate projection of a particular numerical simulation of the system. Some constraints can be placed on the masses of the galaxies, their internal dynamical configurations, and the properties of their relative orbit, including the three-dimensional orientation of the pair. Detailed results are presented for two pairs, K99 and K564. NGC 1587, the brighter component of K99, has the highest rotation rate known for an E2 galaxy, with a value v rot /σ = 0.6. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: (Astronomy) ; Galaxies, Mergers, Galaxy Evolution, Colliding Galaxies, Numerical Simulations, Merging Galaxies Degree Grantor: California Institute of Technology Division: Physics, Mathematics and Astronomy Major Option: Astronomy Thesis Availability: Public (worldwide access) Research Advisor(s): Gunn, James E. Group: Astronomy Department Thesis Committee: Oke, J. Beverly (chair) Gunn, James E. Sargent, Wallace L. W. Goldreich, Peter Martin Defense Date: 16 June 1982 Non-Caltech Author Email: kirk.borne (AT) gmail.com Funders: Funding Agency Grant Number NSF UNSPECIFIED Caltech UNSPECIFIED Record Number: CaltechETD:etd-09082008-132748 Persistent URL: https://resolver.caltech.edu/CaltechETD:etd-09082008-132748 DOI: 10.7907/8PH6-6R18 Related URLs: URL URL Type Description https://doi.org/10.1086/162709 DOI Article adapted for Chapter 1. https://doi.org/10.1086/166453 DOI Article adapted for Chapter 2. ORCID: Author ORCID Borne, Kirk Daniel 0000-0001-5608-255X Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 3390 Collection: CaltechTHESIS Deposited By: Imported from ETD-db Deposited On: 12 Sep 2008 Last Modified: 23 Jul 2025 17:12 Thesis Files Preview PDF - Final Version See Usage Policy. 10MB Repository Staff Only: item control page

THE STRUCTURE ÄND EVOLUTION OF INTERACTING BINARY GALAXIES Thesis by Kirk D. Borne In partial fulfillment of the requirements for the Degree of Doctor of Philosophy California Institute of Technology Pasadena, California 1983 (Submitted June 16, 1982) Ad Majorent Dei Gloriam "Lift your eyes and look. Who made these stars if not he who drills them like an army, calling each one by name? So mighty is his power, so great is his strength, that not one fails to answer." Isaiah 40:26 Ill ACKNOWLEDGEMENTS It gives me great pleasure to thank all those who made this work possible. Thanks are especially extended to those whose contributions are most memorable: — foremost, to Jim Gunn, for his invaluable advice, counsel, and encouragement over the last several years. — to Bev Oke for many hours of enlightening discussion on all manner of binary systems. — to Wal Sargent, Peter Goldreich, Kip Thorne, and Bohdan Paczynski for suggestions concerning this and other research endeavors along the way. — to Roger Griffin for opening my eyes to the understated elegance of coude photographic spectroscopy, for teaching me how to do it, and for countless taffies. -- to the Director of the (then) Hale Observatories for generous allocations of telescope time. — to the Director of the Ritt Peak National Observatory for the telescope time required to bring this undertaking to fruition. — to the fine Observatory staff and mountain crews for their assistance. — to Barbara Zimmerman for her help at many different computing facilities. — to Helen Holloway for caring and for helping me get into this place. — to Lilo Hauck and to Helen Knudsen without whom the place would surely be at a loss. iv — to the National Science Foundation for a graduate fellowship. — to the California Institute of Technology for additional financial and material support. — to Al Hiltner and the University of Michigan for employment during the latter stages of the work described here. -- least, but not last, to Stirling L. Huntley for his inimitable style and wit. — to my fellow graduate students for sharing their knowledge and enthusiasm, especially Richard "Bade", John "Hassler", Doug R o , Russell, the Donz (even if he is a Ot‡S f∏ Cornhusker fan), Wild Bill, PJ, Matt, Bruce the Brit, Abi, and comrade Alex. — and finally, yet most emphatically, to my wife, Colleen, for much more than words could ever express. V ABSTRACT A numerical code has been constructed for the study of the evolution of interacting binary galaxies. This "multiple three-body" algorithm (MTBA) essentially involves N concurrent three-body integrations; it is the sum of the interactions of the perturber with each of the N particles comprising the primary "galaxy" that determines the binary orbital evolution. MTBA incorporates a violent relaxation phase that allows the particles to redistribute themselves in the gravitational field of the perturber prior to the full binary orbital evolution calculation. This redistribution is important for systems with an initially strong tidal potential; their predicted merger times are 50-100% larger than previously estimated. Tabulated merger times for both circular and elliptical orbits demonstrate that typical close binary galaxies will merge in about twice g their initial binary orbital periods, roughly 10 years. This timescale depends strongly on the binary separation, weakly on its mass ratio, and less on initial galaxy rotation than has heretofore been suggested. A specific interacting binary simulation is described in detail in the first paper. Surface brightness maps, structural parameters, and rotation properties are described at various times during the evolution, followed by a discussion of the properties of the merger remnant. Many of the results reported here are consistent with those obtained from the larger, more expensive N-body simulations. VI For the second paper, MTBA is altered so that each "galaxy" is represented by a configuration of test particles. Simulated tidal interactions between these "galaxies" are matched to the measured properties of real binaries in an effort to determine the physical state of each observed pair. The catalog of isolated galaxy pairs prepared by Karachentsev has been culled for its E-E constituents. CCD pictures and long-slit spectroscopic observations have been obtained at KPNO for the most photometrically distorted of these systems. Such pairs very likely represent strongly interacting, physically associated, binary elliptical galaxies. Radial variations of rotation velocity and velocity dispersion are extracted from the spectroscopic data for each of the two galaxies in a given pair. The combined spectroscopic and photometric data for that binary can be reproduced by an appropriate projection of a particular numerical simulation of the system. Some constraints can be placed on the masses of the galaxies, their internal dynamical configurations, and the properties of their relative orbit, including the threedimensional orientation of the pair. presented for two pairs, K99 and K564. Detailed results are NGC 1587, the brighter component of K99, has the highest rotation rate known for an E2 galaxy, with a value vrot∕σ = 0.6. VI 1 CONTENTS Acknowledgements iiï Abstract v Preliminary Remarks 1 Chapter 1. A NUMERICAL MODEL AND PRELIMINARY RESULTS 2 Introduction Details of the Calculation Equations of Motion Energy Conservation Units Coordinates and Symmetry Numerical Intergration The Galaxy Model Initial Conditions Mixing and Relaxation Phases Some Results of the Calculation Merger Times for Circular Orbits Effects of Initial Conditions Effects of Rotation Merger Times for Non-Circular Orbits A Specific Example The Surface Mass Distribution The Structural Parameters Velocity Measurements The Spin Parameter Summary and Discussion Appendix A. The Energy Integral Appendix B. Postscript on NumericalProcedures Tables References Figures 3 12 12 13 16 16 17 18 21 24 28 28 34 36 39 42 42 47 49 54 57 65 67 68 72 76 MATCHING MODELS TO OBSERVATIONS 92 Chapter 2. Introduction Observations and Reductions Spectroscopy Pictures Summary of Observations Results from the Data Analysis R99 (NCC 1587/88) K564 (NGC 7236/37, 3C442, Arp 169) 93 95 95 98 99 101 101 109 VI 1 1 The Physical Model The Multiple Three-BodyAlgorithm (MTBA) The Galaxy Potential Flattening Rotation Projections Free Parameters and the MatchGame Simulations of K99 Simulations of K564 On the Evolutionary Statusof BinaryGalaxies Tables References Figures 114 114 117 118 118 121 125 134 139 143 147 154 155 1 PRELIMINARY REMARKS This thesis is the result of investigations of interacting binary elliptical galaxies. The work is described in two chapters, each a self-contained article to be submitted to The Astrophysical Journal. The first chapter contains a complete description of the numerical model and presents some of the theoretical results derived from its use. The second chapter describes spectroscopic and photometric observations for two select cases of interacting ellipticals and provides an analysis of these data in terms of specific numerical models. The observational work described in the latter chapter was done in collaboration with John G. Hoessel. 2 CHAPTER 1 A NUMERICAL MODEL AND PRELIMINARY RESULTS Submitted to The Astrophysical Journal in March, 1982 Author: Kirk D. Borne 3 I. INTRODUCTION Gravitational interactions between galaxies have played an important role in the evolution of the observed universe. The presently observed distributions and structures of cluster galaxies suggest that their histories have been dominated by interactive gravitational dynamics. White and Rees (1978), for example, have presented a scenario for the formation and clustering of galaxies that proceeds through both gravitational and gas-dynamical processes. Oemler (1974) has stated that the galaxy distribution and content of a given cluster may be a measure of the dynamical age for that system. Though the results of Melnick and Sargent (1977) and of Dressier (1978,1980) argue against simple interpretations of cluster data, apparent evidences for dynamical evolution are again noted. The observations of Dressier (1978) and others (e.g., Rood and Sastry 1972; Austin and Peach 1974) indicate that the flattening of cD galaxies is commonly aligned with that of their local galaxy distribution, which fact probably has its explanation in the dynamical friction hypothesis (Binney 1977). That hypothesis has been employed by Ostriker and his co-workers (Ostriker and Tremaine 1975; Ostriker and Hausman 1977; Hausman and Ostriker 1978) to describe the brightest galaxy in some clusters as an aggregate of merged stellar systems. Observations of brightest cluster galaxies by Oemler (1976), Rood and Leir (1979), and Hoessel (1980) are consistent with this hypothesis. Other observational 4 studies have supported the "structural evolution by means of gravitational interaction" hypothesis: in clusters (e.g., Strom and Strom 1978), in groups (e.g., Hickson, Richstone, and Turner 1977), and in pairs (e.g., Faber 1973; Kormendy 1977). The lack of diffuse background light in some groups studied by Rose (1979) argues against the occurrence of much dynamical evolution there. In binary galaxy systems tidal processes affect the evolution through the infusion of orbital energy into the internal degrees of freedom. In the close encounter of two galaxies a redistribution of matter, energy, and angular momentum will affect the observables for each galaxy. For example, the peculiarities in the colors of some spiral pairs (Sharp and Jones 1980, and references therein) can probably be explained in terms of the internal response of the gas to external gravitational perturbations (Larson and Tinsley 1978). Various kinematic properties of binaries are similarly describable in terms of a simple gravitational interaction hypothesis (e.g., that of Toomre and Toomre 1972, hereafter TT). Comparisons between specific examples of gravitationally interacting galaxies and such models are reported by Lynds and Toomre (1976), Combes (1978), Rots (1978), Byrd (1978), van der Hulst (1979), and Stockton and Bertola (1980). Past galactic encounters can also explain the observed curvature in some extragalactic radio jets (e.g., Blandford and Icke 1978; Vallee, Bridle, and Wilson 1979). 5 If the encounter of two galaxies is close enough (a few radii), then the highly inelastic nature of the tidal interaction will render a high probability for merger in a remarkably short time. Significant merger rates have been found in a variety of imaginative simulations of interacting galaxies (e.g., Holmberg 1941; Alladin 1965; Toomre 1974, 1977; van Albada and van Gorkom 1977; Yabushita 1977; White 1978; Alladin and Parthasarathy 1978). Binary orbital evolution and subsequent merger are particular results predicted by the dynamical friction hypothesis, another possible result of which is the observed correlation of luminosity with separation in a sample of galaxy pairs (Turner and Ostriker 1977; Ostriker and Turner 1979; see, however, White and Valdes 1980). TT suggested further that most ellipticals may be the remnants of past merger events among disk galaxies (see also Toomre 1977; Tinsley and Larson 1979; compare with Marchant and Shapiro 1977). Because massive galactic halos should have significant cross-sections, their existence would imply an even higher predicted frequency of mergers among interacting galaxies (White and Sharp 1977). It is therefore possible that a large number of pairs existed in the past that have since coalesced, with the presently observed binaries soon to follow suit (Toomre 1977; White and Sharp 1977). To this diminishing population of binaries are added new binaries formed through random, inelastic two-galaxy encounters. 6 It appears that one or more of the following conclusions is unavoidable: either (i) dynamical friction has been a dominant factor in the evolution of binary galaxies, so that many more binaries existed in the past than exist today, with the result that a non-negligible fraction of all observed galaxies are actually remnants of past merger events (Toomre 1977); or (ii) massive galactic halos are uncommon, thereby strongly diminishing the predicted frequency of mergers; or (iii) the merger process itself is not nearly so efficient as is now believed. Because statements such as these strongly influence current extragalactic and cosmological research, it is essential that we further our understanding of the dynamical friction mechanism. Reported here are new attempts (a) to calculate the efficiency (i.e., the timescale) of this process within the context of binary galaxy systems and (b) to find physical parameters that describe both the evolutionary state and the ultimate remnant of these interacting pairs. Detailed evolution calculations for individual galaxies in clusters are difficult since each galaxy probably encounters several different galaxies in a Hubble time; idealized theoretical studies of this problem include those of Gallagher and Ostriker (1972), Richstone (1975, 1976), Hausman and Ostriker (1978), Knobloch (1978), and Da Costa and Knobloch (1979). The problem of the evolution of an isolated pair of galaxies is more tractable, and the results found here will be applicable in some measure to the 7 problem of the evolution of galaxies in clusters. Theoretical studies of this type are numerous (e.g., Wright 1972; Biermann and Silk 1976; White 1979b; Roos and Norman 1979; Miller and Smith 1980 and references therein). The Why, then, should yet another study be pursued? two major reasons are: 1. Currently accepted dynamical friction timescales 8 9 are disquieting. In times as short as 10 or 10 years significant changes in the orbit and observed structure of some binaries' are predicted. It is important that these startling numbers be recalculated and corrected if necessary. That is accomplished here. Merger times are presented for a wide range of initial binary orbital parameters. The results in this paper demonstrate that the numbers in current use are probably not too far wrong, although an upward correction of 50-100% seems to be in order for the more strongly interacting systems. 2. Techniques now in use for the study of interacting galaxies fall roughly into two categories: (i) expensive, yet physically realistic simulations; and (ii) inexpensive, yet physically unrealistic calculations. A method is needed which will bridge the gap between these two sets. The algorithm described here provides such a tool. It is labeled the "multiple three-body algorithm".(hereafter MTBA) by Borne (1979). Its high efficiency permits many simulations to be run at low cost. These provide a solid base on which to build conclusions about the fate of 8 interacting binaries. Through the approximation of some hitherto ignored physical effects the model presented here also provides a closer approximation to physical reality than did many of the earlier studies. For example, only the dynamically self-consistent calculations have included the instantaneous back-reaction of the perturbations in the test "galaxy" onto the trajectory of the perturbing "galaxy"; MTBA handles this easily. Another overlooked physical necessity that is met here is the requirement that the initial "star" distribution be in equilibrium with the nearby perturber. The use of spherically (or otherwise) symmetric initial models in strongly-interacting binary simulations is physically unsatisfactory. Efforts to study interacting galaxies are of three types. There are analytic calculations (e.g., Toomre 1977), semi-analytic calculations (e.g., Sastry and Alladin 1977; Da Costa and Knobloch 1979), and numerical simulations (e.g., TT and Miller and Smith 1980, hereafter MS); the latter provide a more faithful physical representation of the dynamics than do studies of the other two types. It is possible to further categorize each technique according to its major physical approximation: impulsive, restricted three-body, or self-gravitating. The present simulations are more realistic than were those that were based on the "impact approximation" (e.g., Biermann and Silk 1976), yet MTBA is just a simple physical modification of the otherwise physically unreasonable restricted three-body code employed 9 by TT and others (e.g., Gutowski and Larson 1976). Yabushita (1977), who also worked with the latter code, nearly discovered the MTBA "trick" in his study of the capture hypothesis for the formation of binary galaxy systems. Because MTBA does not include the self-gravity of the test particles (see details in §11), it cannot compete on the same grounds with codes that do (e.g., MS and Villumsen 1982). For example, no consideration is given to the effects of the coherent contraction of the mass distribution that occurs during close passage (see MS). Nor can the present investigator respond favorably to the additional warnings of MS that the test "galaxy" and the perturber are not always physically distinct systems. It is not clear however that these are disadvantages, since N-body codes do not always adequately suppress two-body relaxation. MTBA specifically avoids the two-body effects, while to some degree the global violent relaxation is provided for. The advantages that MTBA has over the self-gravitating models of MS and others (e.g., White 1978, hereafter W78) are: (i) these latter investigators use spherically-symmetric initial "galaxy" models that are not in equilibrium with the nearby companion, whereas MTBA introduces a "relaxation phase" to provide for this (see Sllh); calculations are (ii) the self-consistent model of necessity large-scale and expensive endeavors, whereas MTBA is inexpensive and of moderate size; (iii) N-body codes either permit very few particles (e.g. , Lauberts 1974, Roos and Norman 1979) or have limited phase- 10 space coverage (e.g., MS), whereas MTBA has no such (iv) MTBA easily reproduces many of the limitations; and results obtained from those simulations, such as merger times, rotation-flattening correlations, and the remnant mass distribution. Its efficiency not only makes MTBA an attractive tool, but the efficiency also renders the model amenable to comprehensive comparisons with observational data. This is possible since dozens of test cases can be run at low cost on a reasonably-sized computer and because there is no major constraint on the number of test particles that can be used; the statistical significance of the numerical results clearly increases with N, the number of test particles. Chapter 2 will present comprehensive surface photometric and spectroscopic observations for two real pairs of elliptical galaxies (from the catalog of Karachentsev 1972) and will compare those data with the numerical simulations. Later papers in this series will present (i) a detailed discussion of a large number of the physical properties that describe the simulated binaries and their remnants, (ii) complete analyses for about six other Karachentsev pairs, and (iii) a compilation of photometric and spectroscopic parameters for nearly all of the pairs of ellipticals listed by Karachentsev. In Sil of this paper the details of the calculation are presented. Included are the MTBA equations of motion, a description of the model used for the test "galaxy", a 11 discussion of energy conservation, units, coordinate conventions, a look at how the initial conditions were prepared for an encounter, and a few remarks on the numerical computations. in Sill. Preliminary model results are given Included are a listing of merger times for a wide variety of initial orbital configurations, a description of the effects of both the random initial conditions and initial galactic rotation, and an in-depth look at a specific example of an interacting binary system. A summary of the main results and a discussion of how these relate to the results of other investigators are presented in SIV. Appendices provide a derivation of the energy integral from the MTBA equations of motion and some specific remarks about the numerical procedures. 12 DETAILS OF THE CALCULATION II. In this section the basic features of the numerical The equations of motion are calculation are to be found. presented and discussed in Sila. The meaning of the corresponding energy integral is investigated in §IIb. Units are defined in Sllc. Coordinates and symmetries are specified in Slid. Facts about the numerical integration are listed in Slle. Details on the derivation and use of the model for the primary "galaxy" are given in Sllf. The Monte Carlo particle initialization procedure and its limitations are described in Sllg. Finally, mixing and relaxation phases are introduced in §IIh; there they are discussed as means by which the transient effects of the initial conditions are curbed. a). Equations of Motion The numerical technique is essentially a simultaneous integration of N three-body equations of motion. In a few test cases N = 2000; otherwise, N = 1000. An initially spherical primary "galaxy" of mass M^ comprises N test particles i, each with mass m^ = M^/N; the secondary "galaxy" (or perturber) has mass M2∙ The equations of motion that define the three-dimensional "multiple three- body" problem are: d⅜i mi dfc2 mi M1 7i ψi1 ^ mi M2 7i φi2 ' (1a) 13 d2?. Μ and where 2 dt2 N = - m∙i Μ m2 if1 72 φ2i = mi m2 ? 7i φi2 ' d⅜1 d⅛2 M1 —2^^ = - m2 —ΪΓ ' dt' dt' Ob) Oc) i-∣, and ?2 arθ fcke position vectors of the test particle i, of the center-of-force in the primary "galaxy", and of the secondary center-of-mass, respectively. physical meaning of equations (1) is as follows. The Changes in the trajectories of test particles are determined by the smooth unperturbed fields of the primary and secondary "galaxies", as derived from the potentials respectively (equation 1a). and φ^2, The force exerted on the secondary is the sum of the individual forces exerted on it by the N test particles (equation 1b). It is in the calculation of this force that the present model contrasts sharply with previous models that were based on either the impulsive or the restricted three-body approximation. Finally, the motion of the center-of-force of the primary is such that the vector quantity + m2^2 rema^∙ns fiχe<3 (global conservation of linear momentum; equation 1c). b). Energy Conservation A simple manipulation of equations (1) results in M 2 N Σ Si] = mi M1 Σi 7i φi1 ' (2) 14 where Φ denotes a velocity. If a self-consistent model for the force field within the perturbed primary were used, then the right-hand side of equation (2) would be equal to zero (Newton's third law) and ‰ - <V = ⅛ pi ∙ (3) As detailed in Appendix A, equation (3) implies the constancy of the energy integral defined by E = i z ιi 1 ‰>' 1 m1v21 ⅛ 1 2 m2v2 mi∑(Mιφiι+ M2φi2).(4) Because the present model is not fully self-consistent, it follows that E =≈ constant if and only if ^-∣ = ⅛ ? ∙ (5) Deviations from a constant value for E are not to be confused with inaccuracies in the numerical integration scheme; the latter are discussed in Appendix B. Nor can these variations in E be identified with physical deviations from total energy conservation since equations (1) do not describe a fully self-consistent system. Test cases show that deviations from a constant value for E are nearly always small (≤2%). Larger deviations occur at times prior to merger during which the secondary is most strongly interacting with the primary. Following the merger event, the system of test particles relaxes to a nearly-symmetric state that satisfies the right-hand side of equation (5). 15 It is at times of intensive tidal interaction prior to merger that the forces exerted by M2 act most to break the approximate symmetry that would guarantee the near-equality in equation (5). It is usually then that several test particles become unbound relative to t which renders an inequality in equation (5) simply because the right-hand side of that equation is no longer meaningful. Consequently, energy is not conserved as it is expressed in equation (4); the value of the parameter E is no longer physically equivalent to the total system energy. These difficulties are ignored here since MTBA cannot adequately represent the dynamics at the moment of the actual merger event so long as the model ignores both the soft, distended stellar distribution within the secondary and the growing perturbations in the force fields of the two "galaxies". Although the MTBA equations of motion resemble those for N independent unrestricted three-body problems, equations (1) differ from such a representation in the following ways: (i) MTBA represents N sets of intermingled three-body equations; (ii) the secondary M2 does not interact directly with the center-of-force in the primary M1, as would be the case in the unrestricted three-body integration; and (iii) M1 does not interact directly with the N test particles. If the last two restrictions were removed by altering equations (1b) and (1c) such that M1 interacts explicitly with both M2 and the test particles, then (i) a new energy integral would exist that would be 16 conserved exactly at all times, and (ii) there would be a two-body relaxation effect among the N test particles. The former is easy to derive from the equations of motion. The latter would exist simply because each particle would then affect the motions of all other particles through their mutual interaction with the center-of-force in . Since N is so much smaller than the number of stars in a real galaxy, any form of two-body relaxation in the model should be and has been avoided, even at the expense of precise model energy conservation. c). Units A dimensionless system of units is used in which G = M1 = r∙, = ,b where R1 is the cutoff radius of the primary "galaxy" (i.e., the radius within which the Monte Carlo algorithm initially places all test particles). The following equations demonstrate the translation from dimensionless to physical units for the case R1 = 20 kpc and M1 = 1011 M@: R 1 time unit = (pη½-)1∕2 = 1.3 x 108 years ; ∙GM (6a) GM ∙, I /ο -1 velocity unit = (—= R—) z = 150 km sec (6b) d). Coordinates and Symmetry The program is built around a non-rotating, stationary Cartesian coordinate system whose origin is fixed 17 by: M1⅞1 + M2?2 = 0∙ τhθ centers of the "galaxies" are initially located on the x-axis and are thereafter constrained to lie in the xy-plane (i.e., z^ = z2 - 0 at all times). This constraint presupposes reflection symmetry through the xy-plane. The solution z1 = z2 = 0 is rigorous only when the mass density around the position (x^, yi, z∣) equals that around (xi, yi, -zp, something which is not true in a small-N calculation. But suppose that for every particle with a given position and velocity there was another particle whose position and velocity were just those of the first particle reflected through the xy-plane. Upon integrating the equations of motion for this configuration, one would find a solution identical to the one found originally by just ignoring the z-components of equations (1b) and (1c) as suggested. e). Numerical Integration A fifth-order predictor-corrector algorithm based on Numerov's formula is used to integrate equations (1) (see Lapidus and Seinfeld 1971). As Appendix B describes, the value chosen for the step size is h » 1/32 (time units). An average binary simulation is integrated for about 20 time units. It is executed in about 50 seconds on the Lawrence Berkeley Laboratory CDC 7600 computer. An equal amount of computing time is expended in the generation of detailed listings of the binary diagnostic parameters, most of which will be discussed in a later paper. Some of the more immediately useful results are discussed in Sill of this 18 paper. The rest of this section is devoted to a description of the initial galaxy model, how it is generated by a Monte Carlo algorithm, and how it is fed into equations (1). f). The Galaxy Model A realistic potential Φ(r) = φ∙1 is desired. The choice of an analytic form for Φ(r) was made in the manner described in this subsection. Assume, first of all, that the phase space density of particles in the primary "galaxy" is given by a spherical, isotropic distribution function of the form f(E) = A exp (-E∕σ2) and f (E) = 0 where and E esc for E < Egsc (7a) for E ≥ Eer, , esc ' (7b) E = ∙∣ v2 + Φ(r) (7c) ⅛ vese<r> + 4lr> * 0 ' (7d) This is the distribution function used by Wooley and Dickens (1961) in their study of model stellar systems. Although King (1966) does not recommend the use of this truncated isothermal distribution function in studies of real stellar systems, the deviations from reality so imposed are definitely no more significant than those, for example, that arise from the assumptions made implicitly through the use of equations (1). An integration of equation (7a) over all velocities determines the space density as 19 4r3 p(r) = 4 π A exp[-- vescir∖ ∫ e×P(-z⅛) dv 2σ (8) The normalization constant A is determined by M1 = 4ιr ∫1 p(r) r2 dr ≡ 1 . 0 For vgsc 2 (9) 2 >> σ , the integral in equation (8) is practically independent of vegc. In that case the density pQ = α2∕2πr2 and the potential Φq = 2p logθ(r) are approximate solutions to the combination of equation (8) and the Poisson equation, 2 V Φ = 4πp. A more general integration of equation (8) yields p(r) = where x 2 exp (x2) [ erf (x) ----- exp (-x2) ] f = -Φ(r)∕σ 2 (10) and erf is the standard error function. With an accuracy of a few percent over the whole range of x, an approximation for the error function is erf (x) ≡ -γyj ∕θ e^t e^χ2 + (1 - e~χ2)2. dt ≈ (11) Roughly, then, the space density is p(r) = ⅛ e* (1 ^ e^χ2)2 . (12) For x2 ≈ -Φ0∕σ2 ∙ -2 loge(r), equation (12) becomes p(r) ≈ -½-≈- (1 - r2}2 2πr for r ≤ 1 , (13) 20 with which the corresponding potential can be derived from the Poisson equation. This is Φ(r) - 2 k [ loge(r) - j2 + ⅛4 ÷ ⅛ ] , 2 where Φ(1) = 0 and k = σ (r=0). (14) The choice of zero point for Φ(r) guarantees that all particles will be contained within the sphere of radius R1 = 1 (i.e., vesc = 0 at r = = 1; see equation 7d). In order to avoid a singularity at r = 0 and to maintain reasonable numerical accuracy, a softening parameter a1 is introduced into the expression for the potential given above. This change yields: φi1 = Φ(r) = k 1 [loge(r2+a 12)- ∣r2+ γ⅞ ] + k2, for r*1'* {15a) φi1 = Φ(r) = 1 - ∙p , r > 1 . (15b) If M(r) denotes the mass contained within the sphere of radius r, then dΦ∕dr = M(r)∕r. Hence, k1 and k2 are determined from the boundary conditions, Φ(1> = 0 and M(1) = M1 = and 1. A value a1 = 0.05 is used, so that k1 ≈ -J ( ------2---- 2 - ~ )~1 ≈ 0.9419 (16a) k2 β k1 C ⅛ ^ lo9e (1 + a12) 3 ≈ θ∙5314 ∙ (16b) Equations (15) give the potential φi1 as it is used by MTBA, where r = force field is The density that corresponds to this 21 r2 + 3a.2 P<r> = 7j [ (r2 ÷ a,2)2 1 p(r) , - 2 + r = 0 ] ( 17a) for r ≤ 1 ; ( 17b) for r > 1 . A softened point-mass potential is used for φ^2' given by Φi2 « - [(?; - + a 2 ]"1/2 ( 18) The force that arises from this potential has a maximum magnitude F2(max) = 2M2∕(32^2a2^) at r^ = 2 ^2a2∙ The maximum magnitude of F1 = -Μ^^φ^ 1 occurs at r^1 ≈ a1, with a value F-,(max) ≈ kιa∙j~1 ~ 18∙84∙ reasonable to assume that, for M1 = ri2 ≈ aτ lt is physically = 1, F2(max) ≈ 20 at Equality of F1(max) and F2(max) for M1 = M2 requires a2 ≈ 0.02. So that a2 would be smaller for less massive "galaxies" it was arbitrarily decided to choose a22 = 0.02M2∙ This value of a2 has the property that F2(max) is independent of M2 and has a value about equal to 19.25. g). Initial Conditions Initial positions and velocities for each of the N test particles are obtained by a Monte Carlo selection procedure. From a rectangular distribution of pseudorandom numbers, the radial coordinate r, cosθ (where θ is the polar angle), and φ (the azimuthal angle) are tentatively chosen The coordinate r is so chosen because -2 the density law (equation 17a) varies approximately as r for a given particle. 22 multiplied by an energy cutoff term, as shown in equations (8) and (13). In order to accommodate the energy cutoff, the first step is to choose particle velocity components v , v , and v from a Gaussian distribution with variance y z 2 σ = k1 (from equation 16a). From equations (7c) and (15a), the corresponding energy E is calculated. If E ≤ 0, then these coordinates (position and velocity) are assigned to the particle; otherwise, they are rejected and new coordinates are selected, subject to the same negative- energy condition. Figure 1 presents a comparison between the particle density distribution generated by the above procedure and the theoretical density curve given in equation (17a). Systematic deviations are to be expected since that equation is only a rough approximation to equations (8) and (10). The latter equations form the basis for the Monte Carlo selection algorithm described in the previous paragraph. Plotted in Figure 1 are: (i) N(R) = 1000 [M(r+δr) - M(r)] (the theoretical curve), specifying the number of particles in a shell of thickness δr = 0.04 with mean radius R = r + fir/2; (ii) open circles, indicating the number of particles placed in each shell by the Monte Carlo algorithm, averaged over five different sets of initial conditions, with error bars indicating the standard deviation of a given value about the mean; and {iii) solid circles, indicating the number in each shell as determined by an average over 85 time steps during the mixing phase of one particular initial 23 configuration (see the next subsection for a description of the mixing and relaxation phases.) The Monte Carlo "galaxy" simulates the theoretical density distribution very well, except for a couple of discrepancies. For one, the number of particles in the innermost shell greatly exceeds the theoretical value. Second, the points at mean radii R > 0.4 are systematically below the curve. That the selection algorithm preferentially places particles in shells at small radii at the expense of those shells at greater radii is a result of the following approximation used in the selection procedure. The Monte Carlo technique defined above actually simulates a "galaxy" whose density law is where x 2 p(r) = —⅛-≈ [ erf(x) exp(-x2) 3 , 2πr2 ir1/2 (19) = -Φ(r)∕k, and Φ(r) is given in equation (15a). A check of the distribution function N(R) that corresponds to equation (19) indeed finds that it agrees with the placement of the points in Figure 1, including the discrepancies. The difference between equations (19) and (17a) is in the way the core of the "galaxy" is treated. The core appears explicitly in equation (17a), whereas it appears only implicitly in equation (19) (in the energy cutoff term). It is therefore not surprising that the simulated galaxy has a higher density of particles at small radii than is expected from equation (17a). Because of the added complexity that would attend the inclusion of the core in the Monte Carlo algorithm, it was decided to ignore the problems with the 24 simulated "galaxy" and to proceed with the integration as is. Since the particles respond in time to a potential that explicitly contains the core (equation 15a), their distribution tends to approach the density that is the source of that potential. The filled circles in Figure 1 demonstrate this. Two useful numbers that can be derived from the data plotted in Figure 1 and from the density law in equation (19) are the "galaxy" half-mass radius R^ and the projected half-mass radius Rβ (i.e., the "effective radius"). Their values are Rh ≈ 0.26 h). and Rβ ≈ 0.20 . (20) Mixing and Relaxation Phases In order to minimize any systematic effects of the initialization procedure the initial model of the "galaxy" is run through a mixing phase where the secondary "galaxy" is totally ignored, the motion of the primary center-of- force is nil, and the origin of the coordinate system is set at the center of M1∙ For a pre-specified time, typically 2-5 time units, the N test particles are therefore allowed to move exclusively in the gravitational field of the primary. The filled circles in Figure 1 represent the average particle density distribution during such a mixing phase. A kinematic description of the mixed, isolated "galaxy" is reserved for a later paper. 25 The full binary simulation is performed with the origin of the coordinate system at the binary center-of-mass (see Slid). Because the initial values of the position and velocity coordinates for the particles refer to an origin at the center of M1, the values of the position and velocity coordinates of M1 are added to those of each particle at the end of the mixing phase. The system of N particles that comprises the primary "galaxy" is thereby set in motion, thus marking the beginning of the interaction phase. Aside from statistical fluctuations, this system is in a non­ rotating configuration. A couple of binaries were run with a spinning primary "galaxy". The latter rotated synchronously with the binary revolution; to the peculiar post-mixing velocity of each particle was added a velocity component equal to that of a rigid rotation of the primary around the binary center-of-mass at the binary angular velocity. Some results from these runs with rotating particle configurations are discussed in SlIIc. One of the difficulties attending earlier numerical binary interaction studies is the problem of transients; such problems were mentioned in SI. It is not physically realistic to place a massive perturber on a trajectory that brings it very close to another "galaxy" whose mass and velocity distributions correspond to an equilibrium configuration at infinite separation. Unphysical transient effects arise here as soon as the force field of the secondary is turned on; how much these have affected the 26 merger rates calculated by W78 and others is not yet known. The present study attempts to answer this last question. For a fixed period of time the MTBA equations of motion are replaced by the equations of the restricted three-body The duration of this relaxation phase is called problem. the relaxation time Trθlaχ ; its value is chosen after the initial binary orbital parameters are selected. T Γ6ι3X n is the elapsed time during which M1 and M2 respond only to the total mass of the other as though all of the mass in each were confined to a point at its center. respond to M1 and The test particles as shown in equation (1a). A discussion of the procedure by which the best value for 7relax can est^matθ^ is given in SlIIa. During the relaxation phase the binary follows a simple Keplerian trajectory. Before the relaxation phase is actually begun, however, the solution of the two-body problem is used to take the secondary from its initially chosen position to the point on its trajectory that precedes that initial location by an elapsed time equal to the relaxation time. This guarantees that following the relaxation phase the full binary interaction phase will begin with the orbital parameters of the secondary precisely equal to those initially prescribed. The relaxation phase begins as soon as the perturber is placed back on its Keplerian trajectory by a time Trelaχ. It is at this point that the restricted three-body calculation is initiated. During that calculation the test particle mass, energy, and 27 angular momentum distributions change in such a way as to make accommodation for the nearness of the perturber. Thereby the primary "galaxy" enters into a reasonable equilibrium with the secondary, and thereafter the full interaction phase is ready to begin. The latter phase corresponds to the integration of equations (1) and continues until either a merger event has occurred or a prescribed time has elapsed. For every test run there is a run in which τreιaχ = 0 and one in which Tre^aχ > 0. Comparisons between the "relaxed" and the "unrelaxed" runs are made in SSlIIa-d. In summary, mixing and relaxation phases are introduced into the binary evolution calculations in order to avoid those non-physical transient phenomena that depend on initial conditions. Before the full simulation of the binary is undertaken, a restricted three-body integration is carried out with a point primary M1, a point perturber M2, and N temporarily massless test particles. This permits the particle mass, energy, and angular momentum distributions within the primary to accommodate the presence of the perturbing secondary in a manner more nearly physical than heretofore attempted. 28 III. SOME RESULTS OF THE CALCULATION In this section are described preliminary results of the integration of equations (la-c) for some 30 different initial orbital configurations. Sections IIIa-d will present a discussion of merger times and SlIIe will describe the "observables" that have been derived from a specific test run. Merger times are tabulated for both circular (Sllla) and non-circular orbits (Sllld). The effects of the random initial conditions and initial galactic rotation are discussed in §IIIb and SlIIc, respectively. For the specific model simulation of SlIIe, there are discussions of surface mass distributions (Sllle.i), radial ellipticity variations (Sllle.ii), rotation velocity and velocity dispersion measurements (Sllle.iii), and angular momentum properties of the merger remnant (Sllle.iv). a). Merger Times for Circular Orbits A merger time is defined to be the time elapsed between the end of the relaxation phase (i.e., the beginning of the interaction phase) and the moment at which the perturber becomes energetically bound to the core of the primary, roughly at a binary separation equal to 0.25. This time is usually well-determined since the merger event is typified by a marked increase in the binding energy of the secondary relative to the primary within a narrow interval of time. Table 1 presents merger times for runs with initially circular orbits. The times are presented as a 29 function of both perturber mass M2 and initial binary separation a, in units as defined in Sllc. Each case that is entered in Table 1 has a relaxation phase duration equal to one initial binary orbital period. Merger times are nearly independent of M2, except in those cases of initially -3 weak tidal interaction. A "weak tide" limit, M2a = 0.20, is established empirically from the location of the points that are plotted in Figure 2. That figure demonstrates the effect of the relaxation phase on numerically-determined merger times. Merger times T2 are those given in Table 1 (that is, for interactions with τreιax equal to one initial orbital period); merger times T1 are for those same configurations (i.e., M2 and a unchanged), but with no relaxation. One sees therefore that the non-physical merger times T1 are usually underestimates of the more nearly physical times T2∙ The relation T2 = 1.5T1 provides a fair estimate of the magnitude of the effect for the strongly interacting binaries (filled circles). For values of the -3 tidal force strength Mja < 0.20, T1 approximates T2; these are the open circles in Figure 2. has a strong dependence on M2. It is also found that T1 For example: if a = 1.5, then T.j = 27.5, 19.5, 15.5, 10, and 7.5 for M2 = 0.25, 0.5, 0.75, 1, and 2, respectively. A comparison of these numbers with the third line in Table 1 shows that the strong dependence of T1 on M2 is an artifact of the numerical integration scheme that ignores the activity of transients (see discussion in Sllh). 30 It is convenient to have a working approximation for the times T2 that are given in Table 1. For the range 0.20 ≤ M2a^3 <1.0, T2 = 7 (M2 a"3)-0·25 a1,5 . (21) Returning to physical units and using equation (20), one finds that this becomes T2 ≈ 2 [ m2⅞ M, -0.25 ) 1 1/2 ( GM. ) (22) The first factor in equation (22) measures the strength of the tidal force in units of the internal force within the primary "galaxy"; T2 depends only weakly on this ratio, contrary to the behavior found for T1∙ The second factor in equation (22) approximates the dynamical timescale for the binary system as a whole. For M2 ≥ 0.10, a better approximation is log1θ T2 ≈ 0.09 P + 0.5 , (23) where P is the initial circular orbital period of the binary and both T2 and P are measured in dimensionless units. Each of these last three equations demonstrates the strong dependence of T2 on the initial binary separation a, which fact is also obvious in Table 1. From those equations one derives a merger timescale of order 10 years for a pair of galaxies whose relative orbit is nearly circular, whose mass ratio is near unity, whose separation equals a few effective 31 radii, and whose component galaxies are of typical galactic dimensions. A physically meaningful value for the "relaxation time" was first determined empirically. A particular orbital configuration was chosen for a detailed study of the dependence of the merger time τmerge °n τreιax∙ found that I" grew as T lt was . „ was increased in value from 0 to P/2 to P, but was nearly unchanged when τre^ax was increased from P to 3P∕2. The value τre2ax = p was therefore adopted for use in all bound-orbit simulations; various initially-unbound orbits were also run, some of which were unrelaxed and some of which had values of Tre^aχ in the range from 4 to 10 time units. These empirically chosen values for "relaxation time" were later justified physically in a study of the time-dependence of the net torque applied to M2; this is the torque applied by the N test particles to M2 relative to the systemic center of mass. In Figures 3 and 4 the torque applied to M2 is plotted against time for two different orbits. Each point actually represents the negative value of a summed torque that is the sum of the values of the tidal torque calculated at each of 16 consecutive time steps. Each point then corresponds to an average over the elapsed time ∆t = 16h = 0.50. Figures 3 and 4 each present the summed torque plotted against time for two different values of τrejax∙ The left-hand side of each figure shows the variation in 32 torque during the interaction phase of a run with zero relaxation; the right-hand side shows the torque variation both during the relaxation phase (the shaded region) and during the subsequent interaction phase for a run in which τrelax ss P* ln F^9ure 3 the orbital parameters are M2=1 and a=1.5; in Figure 4, M2=2 and a=1.5. The "tidal interaction strength" at time t = 0 for -3 the configuration studied in Figure 3 is M2a = 0.30; this value is close to the "weak tide" limit (see Figure 2). A slight rise in the torque is noted between t = 0 and t = 5 on both sides of Figure 3; these are transients. The sum of the torques applied to M2 during that time is therefore strongly negative, corresponding to a large angular momentum loss by M2. This in turn leads to the rapid merger that is noted in the case with zero relaxation (left-hand side). In the case with non-zero relaxation (the right-hand side), the effectiveness with which the torque reduces the angular momentum of M2 is diminished since the transients have all but disappeared by the time the interaction phase is begun. The instantaneous torque oscillates between positive and negative values for some time after the relaxation phase, hence the decreased capture efficiency. In Figure 3 merger times are 16.5 and 10 for runs with and without relaxation, respectively. A more striking demonstration of the transient effects is shown in Figure 4. The "tidal strength" at t = 0 is M2a~3 = 0.59, twice that of the orbit presented in 33 Figure 3. The peak value of the transient torque in Figure 4 is likewise two times that of Figure 3. There is no doubt in this case that the rise in torque at early times is a non-physical transient effect; there is no other peak of similar height until just prior to merger. Merger times are 16.5 and 7.5 for the runs with and without relaxation, respectively. Between the right-hand and left-hand sides of Figure 4 there is a difference in the peak torque values plotted at the moments of the actual merger events (indicated by arrows). This difference is due mainly to the nearly instantaneous capture of M2 in the case with zero relaxation, as compared to the slower capture of M2 in the second case. The shaded regions of Figures 3 and 4 have already been identified with relaxation phases of duration τrelax = P* since t^e transient torque peaks are completely enclosed within the boundaries of those regions, perhaps a shorter relaxation phase would be allowable. There appear however to be small secondary rises in the torque during the latter portion of each relaxation phase; these may or may not be transient phenomena in the same sense as the primary peaks are transients. In order to provide for some margin or error in our estimate of the best value for Trg^aχ, a value equal to one orbital period is usually chosen. The numerical algorithm can thereby suppress the major transients that develop when the massive perturber is placed near to the initially spherically symmetric particle 34 configuration that comprises the primary "galaxy". The empirical choice of Trelaχ = P is therefore physically justified from an investigation of the torque curves (Figures 3 and 4), and values like one-half of the initial orbital period are probably okay also. The conclusion derived from this is that the merger times found by previous investigators must be modified upwards if the initial tidal forces in their simulations were very large. In most cases the corrections to the merger times amount to about 50% (see Figure 2). b). Effects of Initial Conditions It is important to determine what effect initial conditions have on the merger times that have been calculated by MTBA. For this purpose many tests were run with one of the orbital configurations specified in Table 1. The initial orbital parameters in all of the test cases discussed in this subsection are identical: M2 » 0.5 and a = 1.5, on an initially circular orbit. For these values, M2a"3 = 0.15, which places this interaction in the "weak tide" domain. number. Table 2 lists all of the various runs by Each run is specified by its initial set of particle positions and velocities. These initial conditions are determined by (i) the initial "seed" value in the sequence of random numbers employed by the Monte Carlo algorithm, (ii) the mixing time of the isolated, spherical particle configuration that was generated by that "seed", (iii) the number of test particles, (iv) the inclusion or 35 exclusion of synchronous, rigid rotation in the particle configuration at the angular velocity of the binary system, and (v) the duration of the relaxation phase. In Table 2, , and are measured in units of the initial relax merge binary orbital period P = 3π; Tjniχ is measured in time units both T as specified in Sllc. The contents of this table are examined below. A comparison of the values of τmerge for runs 2, 3, and 5 demonstrates their virtual independence of the "seed" value, a conclusion also to be drawn from the results for runs 11 and 12. An apparently statistical fluctuation in the distribution of particles used for run 13 caused a substantial deviation in its merger time from the values for found in nearly all other runs. Runs.8 and 11 merge together show that τmθrge is basically independent of the T mixing3 time Tmix w,. . Runs 3, 4, 8, and 10 demonstrate that T is independent of the number of particles N. This merge last conclusion is important since it implies that the effects of small N are to some extent under control. Here is a distinct advantage over the N-body codes that can only marginally suppress two-body relaxation, the source of most dependence on N. Statistical variations in numerically-determined merger times are probably smallest for large N. Runs 4 and 10 (with N ≈ 2000) therefore suggest that τmθrge is independent of τreιaχr as expected in the "weak tide" domain. This independence is further supported by a 36 comparison of the values of T_ ___ sr ι∏eeΓy e in runs 6, 7, and 8. The results from runs 14, 15, and 16 are consistent however with an increasing merger timescale for values of the relaxation time in excess of 1.5P. c). Effects of Rotation A major result to be found in Table 2 is demonstrated in a comparison of runs 1 and 9 with other runs. In these two simulations the primary "galaxy" rotates synchronously with the initial circular orbital motion of the binary system. Such rigid rotation is added to the initial peculiar velocities of the individual particles at the end of the mixing phase (see Sllh). When Tre]_ax = θ» the merger times calculated in test cases without rotation (runs 2, 3, 4, and 5) are about 50% larger than that calculated in the case with rotation (run 1). When τrelax = p, k°wθvθr' t*lθ merger times are independent of rotation (compare run 9 with runs 8, 10, 11, and 12). There are a few other cases that are not listed in Table 2 that were also run with initial rotation; their results substantiate the above conclusions. When physically reasonable merger times are computed through the use of simulations in which Tre^aχ is non-zero, times Tmerge are thereby found to be independent of the rotational properties of the component galaxies. White (1979a) found the opposite to be true, that τmθrge depends strongly on the relative directions of the galaxy spin and binary orbital rotation vectors. Although interesting, his results should be viewed 37 with caution since his simulations were run with zero relaxation and are therefore subject to the transient effects described in SlIIa. Transient torques are most likely to be of greatest amplitude in the parallel-rotation simulations (see Keenan and Innanen 1975). This statement can be substantiated by a comparison of the time-dependent torques that have been calculated for runs 9 and 11, which have synchronously- rotating and non-rotating configurations, respectively. The peak absolute values of the total torque during the first half of the relaxation phase’s in runs 9 and 11 are 0.39 and 0.18, respectively; summed total torques are defined in SlIIa and have units in the dimensionless system that was defined in Sllc. These two numbers indicate that strong spin-orbit coupling occurs in the case with parallel- rotation (run 9). An additional peak in the torque is noted during the second half of each relaxation phase, similar to the secondary peaks noted in Figures 3 and 4. The peak absolute values for the secondary transient torques are 0.37 for run 9 and 0.33 for run 11. That the relaxation phase has indeed permitted the development of a particle distribution in near-equilibrium with the perturber is indicated by the approxiate equality of the latter two numbers. In addition to this near-equality, both runs had τmerge = 1θ' Subsequent values of the summed torques were also very similar. At the five time intervals (∆t « 0.5) that immediately preceded the merger event, absolute values 38 (in dimensionless units) were TORQUE ≈ 1.29, 1.99, 1.69, 3.31, and 4.67 for run 9 and TORQUE ≈ 1.01, 2.08, 1.59, 3.25, and 4.28 for run 11 at the identical interaction times T « 16, 16.5, 17, 17.5, and 18, respectively; here, the interaction time T ≈ 0 corresponds to the end of the relaxation phase (i.e., the beginning of the interaction phase). Since the relaxed configurations have far more similarities in their internal dynamics than do the initial models, it seems that the relaxation phase is accomplishing its intended goal of leading the particle distribution into a configuration that is nearly in equilibrium with the perturber. The result of this discussion is that the relative sizes of the earliest transient peaks in the tidal torques of runs 9 and 11 do indeed suggest that a rotating "galaxy" is initially more strongly coupled to the secondary "galaxy" than is a non-rotating "galaxy". The evidence provided by the similarities in the total overall evolution of runs 9 and 11 reveals, however, the probable non-physical nature of this preferential strong initial coupling. The implementation of a non-zero relaxation phase in binary interaction studies leads to the conclusion that the strong spin-orbit coupling suggested by White (1979b) and others (see also TT, and Keenan and Innanen 1975) may be more an artifact of the numerical technique than fact. Because of the importance of this result in both galaxy and globular cluster studies, the tidal-coupling problem should be 39 pursued in greater detail with models that are both selfconsistent and relaxed. d). Merger Times For Non-Circular Orbits Evolution calculations for several test binaries with non-circular initial orbits are mentioned here. These cases offer more probable representations of real binaries than do cases with circular orbits. The results of these calculations can also be directly compared with the results of other investigators who preferred to study binaries in non-circular orbits. Table 3 presents values of τmerge for various test runs. Each tabulated orbit is specified by (i) the perturber mass, (ii) the semi-major axis, eccentricity, and perigalacticon distance of the initial (unperturbed) trajectory, (iii) the initial binary separation, and (iv) the initial binary orbital period. At least two different values of τreιax are used for each orbital configuration. In the simulations with initially elliptical orbital parameters, Treιax is set equal to zero in one run and is set equal to about one initial orbital period in another run. In those cases in which the orbit is initially parabolic, the values of 1,relaχ vary. In one of these runs -, is again set equal to zero. In the other run(s), the relax perturber is placed back on its trajectory to a point that T is separated by an elapsed time τre∣ax from the point at which the perturber initiated the run that has zero relaxation time. The merger times that are given for runs 40 with different relaxation times, yet with identical parabolic orbit parameters, are thereby measured from the same binary position-velocity configuration. In practice, every value of Tmerge that is given in Table 3 for any parabolic orbit is measured from the moment when the perturber first crosses the boundary r = 1.5, regardless of the number of times the perturber may eventually cross that boundary. The relaxation phase does not seem to significantly affect the merger times for initially parabolic trajectories. The relaxation phase does, however, have a distinct effect on those elliptic orbits that have large initial values of the tidal parameter M2a . As was found in SlIIa for circular orbits, Tma ,„„ is larger for the merge 3 relaxed elliptical configurations than for the unrelaxed configurations. The selection of the initial parameters for the binaries listed in Table 3 was based on the test cases run by W78. For his initial models and in his units, one internal crossing time is τcr s = 5*7 ∙ (24a> In the units and for the initial model that are presented in this paper, one crossing time is 41 τ- ∙ fr>= 0∙42 ' ,24b) where the numerically-determined squared velocity dispersion of the initial model is σ = 1.53 + 0.03 (one standard is taken from equation (20). deviation) and time units are related by τw^1^∙te = 14 χ τBorne. These two ψ]ιe values W from this relation and merge in Table 3 have been computed ^ from the merger times listed by W78; his values were obtained from N-body simulations with these same initial orbital parameters. Here is the important result that the w merger times τmerge ate reproducible with the present technique. MTBA does this at a fraction of the expense required by the N-body algorithm. In those cases listed in Table 3 that have the strongest initial tidal interaction (e.g., M2 = a = 1.0), both the effect and the importance of a non-zero value for τrelax man^est∙ This raises a caution signal to those researchers who wish to simulate strongly interacting galaxies. Though the effects of self-gravity are also of great importance, their exclusion from this study is dictated both by the increased complexity and expense associated with their inclusion and by the unknown two-body relaxation effects thereby introduced. MTBA reproduces the merger times and some of the other results that have been found with the more expensive codes (see next subsection). That is enough reason to pursue it as an alternative simulation algorithm in future studies. 42 e). A Specific Example One particular case is presented in this subsection to serve as an indication of what can be obtained from these simulations for use in comparisons with observations. The parameters of the test case chosen for this study are M2 = 0.5 and a = 0.75, on an initially circular trajectory. value of the tidal strength parameter, M2a The = 1.19, identifies this as a strongly interacting binary. The values Tm^χ = 2 and τrg^aχ = 1.5 are used, to be compared with the initial binary orbital period P = 3.3. These numbers indicate that the interaction phase actually begins at an internal numerical integration time t = Tmi∙ X + τrelax = 3∙5∙ τile mer9er time is 4.0. This implies a numerical integration time t = 3.5 + T __ = 7.5 at the merge moment of the merger event; the definition of this event is found in SlIIa. The numerical dynamical investigation is continued for about six additional crossing times after that event, until t « 10.0. i). The Surface Mass Distribution Figures 5 and 6 present sample contour plots of the logarithmic projected particle distribution. Plotted at various times during the integration, these "surface brightness" maps provide both top and side views of the interaction (Figures 5 and 6, respectively). Each view is centered on the center-of-mass of the primary. Presented below are descriptions both of the procedure by which these 43 maps are constructed and of the information contained in the figures shown here. A single contour plot is the result of a sum of the projected mass distributions calculated at 17 consecutive time steps t = tθ + mh, where the step size h = 1/32 and m = 0, 1, 2, ..., and 16. The time t1g is the time T that is given in the upper left corner of each map. Because a crossing time is Tcr = 0.42, each plot is basically an average over a little more than one crossing time. It is for improved statistics that each map is created as a sum over all intervening time steps, inclusive of the endpoints of that interval (∆t ■ 0.5). The sum of the projections is made in the non-rotating frame. This may result in some smearing in the maps, but it is not obvious that there would be less smearing if the sum were to be made in the rotating frame. At each time step tjn the program determines that number of particles that project into small two-dimensional cells comprising three independent, orthogonal planes (the coordinate planes). These numbers are added into an array N^(i,j), which specifies the number of particles that lie in the (i,j) cell of the projection into the x∣ς = 0 plane. After 16 time steps have elapsed, this array is transformed into another integer array, which is defined by L(i,j,k) ≈ Int [2.5 log1θ N∣ζ(i,j)]. The line printer simply reproduces this large array of one-digit integers in the form of three orthogonal two-dimensional projections of the particle mass distribution. Contour lines are then drawn that enclose 44 equal-valued digits. Figures 5 and 6 are representative of the results that follow from this procedure. According to the definition of L(i,j,k), the contour lines are separated by equal one-magnitude intervals. The only exception is the outermost contour, which is drawn two magnitudes below the next higher level; the omitted contour between those two levels is excluded because it would add more confusion than information. Contours are also not drawn around those particles that have become separated from the main body of the primary "galaxy"; their number is sufficiently small as to render their specific locations largely dependent on initial conditions. The maps that are portrayed in Figures 5 and 6 are to be discussed in time-sequential order. marks the end of the mixing phase. The time t = 2.0 Hence, the first map in each figure is a time-averaged "snapshot" of the projected mass distribution in the isolated, mixed primary. The outermost contour in the side view (Figure 6; t = 2.0) appears flattened relative to any other contour in either view. This is presumably no more than a small-numbers statistical effect. In Figure 5 (t' = .2.0) a comparison between the analytic formulation and the Monte Carlo simulation is presented. Dashed circular contours are drawn there at radii that are predicted from the projection of the density law given in equation (17). The mean radii of the two innermost contours of the simulated "galaxy" agree with the theoretical values; dashed contours are not drawn at 45 those levels to avoid confusion. Except for small-numbers statistical fluctuations the positions of the other contour levels in the simulated "galaxy" are also in agreement with those expected. At t = 2.0 the relaxation phase is begun by placing the perturber on the x-axis at a distance a = 0.75 from the center of the primary. Because Trelaχ = 1.5, the time t ≈ 3.0 marks the point at which two-thirds of that phase has elapsed. Plots at t = 3.0 thus portray the time- averaged appearance of the primary during the middle third of the relaxation phase. Strong asymmetries in the primary mass distribution are seen in those maps. The asymmetries, which are obviously induced by the perturber, are of quite different appearance when viewed from the two directions. In the top view (Figure 5) the amplitude, ellipticity, and rotation angle of the perturbations vary with radius, whereas in the side view (Figure 6) only the amplitude and the ellipticity vary radially, as expected. Figures such as these will help to determine the present state (e.g., the separation and orientation) of some of the observed interacting binary galaxies. Particularly useful for this is Figure 5, from which the radial variations of isophotal ellipticity and position angle can be extracted. The time t = 6.0 occurs roughly halfway between the beginning of the interaction phase (t = 3.5) and the time of merger (t = 7.5). Even though the perturber is still at a radius of about 2Rh and has not yet rendered obvious damage 46 to the inner regions of the primary, the outer regions of the primary are already greatly extended and obviously asymmetric. The appearance of "fingers" in the outermost contour lines are indicative of individual particles whose trajectories are traced out during the 17-step plot­ averaging period. Of all the particles that are detached from the main body of the primary (i.e., the ones not identified in these plots), the majority are flung by the perturber in the directions of positive x and y coordinate values. The asymmetries are therefore more extensive than is portrayed. Time t = 9.0 occurs after the merger event has taken place (t = 7.5) and the particles have had time to respond to the increased core mass of the merger remnant. time each plot shows extended outer contours. At that Each also shows a lack of central condensation; the surface "brightness" of the highest contour level observed at the earlier times exceeds that of the highest level observed at t = 9.0 by one magnitude. The core radius of the "galaxy" has therefore increased. Rigorous conclusions about the magnitude of the core radius variations following a merger event are difficult to make, however. This difficulty is due first to the lack of information on the extended mass distribution that necessarily comprises the perturber M2. The second reason is that self-gravity (ignored here) necessarily affects the final adjustment of the mass distribution following a merger event. 47 ii). The Structural Parameters Some idea of the structural changes involved in a collision/merger episode can be obtained from plots like those given in Figures 5 and 6. For example, diameters a and b can be measured across the major and minor axes, respectively, of each contour level drawn in these figures. Then a mean isophotal radius <r> = 0.5(ab) 1 z/2 and an ellipticity ε = 1 ^ b/a can be calculated for that contour. The latter two parameters are plotted in Figure 7; their values are derived from the surface "brightness" distributions in Figures 5 and 6 at t = 9.0. Open circles indicate the variation in surface "brightness" (in units of magnitudes per square coordinate cell) as a function of mean isophotal radius; these are obtained from the top view of the merger remnant. The level that has surface "brightness" equal to 7.0 has a mean isophotal radius equal to 1.11, which is off the graph in Figure 7. Filled circles in the lower graph specify projected brightness variations in the side view of the merger remnant. In the upper graph of Figure 7, filled circles show the variations in ε as a function of <r>, also extracted from the side-view of the merger remnant. The dashed line in the lower part of the figure corresponds to a surface brightness variation like -2 <r> , whereas the solid curve demonstrates the radial behavior of the projection of the spatial density law given in equation (17). There is not much difference in these two curves at intermediate radii. The top view of the merger 48 remnant behaves very nearly like the projection of the initial density law. This would support the work of Hausman and Ostriker (1978) who used homologous models in their study of merging galaxies. There are two exceptions to this These are (i) the absence presumably homologous behavior. of a central brightness peak in the remnant and (ii) the distension of its outer regions. In the side view the inverse-square law seems to provide the better fit, although the other cannot entirely be ruled out. Here again the innermost and the outermost points are not in accord with either curve. As indicated by the upper graph in Figure 7, an observer in the plane of the binary orbit would probably classify the remnant as E2-3. The apparent increase in ellipticity with projected radius may be real. This conclusion is supported by the additional observation that the distribution of those particles that are stripped from the main body of the primary is flattened with ellipticity ε ≈ 0.4 - 0.5. Since a real galaxy has many orders of magnitude more stars than included in these simulations, the stripped stars would probably be in sufficient numbers so as to be detectable. They would appear at very low brightness levels, however, something like 6-7 magnitudes per unit area lower than the central surface brightness. In some cases then where galaxies are observed to have a significant increase in ellipticity with radius, these may be identified with products of the merger process. 49 Another measure of flattening in the particle distribution is determined from the inertia tensor I∙j = A flattening parameter f is defined by f ≡ (25) 1 - 2i33∕(l1l+ I22) The values for these three moments of inertia at t = 9.0 are I11=<x2>≈ 0.254, ∑22=<y2>= 0∙23θ' aπd Σ33=<ζ2>= θ∙193. (26) These numbers are results of averages taken only over those particles that are within unit radial distance of the center of the primary. A value f = 0.22 is obtained. This slightly underestimates the observed ellipticity, as expected from the remarks of Aarseth and Binney (1978). iii). Velocity Measurements The only other set of observations that is available to astronomers is the set of projected radial velocities and velocity dispersions taken as functions of position across the face of the galaxy. A collection of such observables is presented in Figure 8. They are discussed below. In this figure, the time coordinate is displaced in such a way that T = 0 corresponds to the beginning of the interaction phase. This coincides with the moment labeled t = 3.5 earlier in this section. T ≈ 4.0 denotes the time of the merger event. T = 6.0 is a time after the particle distribution has begun to respond to the increased core mass of the merger remnant. Figure 8 is divided into three sections, each of which includes a graph of the instantaneous primary rotation 50 curve at one of the times defined above. Also given in each frame is a value for σθ, the central projected velocity dispersion. This number is the square root of an average over three orthogonal directions of the squared central projected velocity dispersion. Only those particles within a projected radius R = 0.04 of the center of the primary are involved in this calculation. The numbers shown in Figure 8 suggest that there is a significant increase (30%) in the value of σθ between T = 0.0 and T = 4.0. Variations in the value of σθ at times following the merger event appear to be statistical fluctuations. The calculation of σθ includes an rms average over the three values σθ(x), σθ(y), and σθ(z), which are the central projected velocity dispersions along the specified coordinate axes. For comparison with the numbers shown in Figure 8, values can also be given for σθ(z) and σθ(2). The latter is defined as the root mean of the squares of σ^(x) and σ^(y). Values are: σ^(2) = 0.90, 1.26, and 1.15 and σθ(z) = 0.76, 1.02, and 0.98 at times T = 0.0. 4.0, and 6.0, respectively. The percentage variations in σθ(2) equal those in σθ, whereas those in σ^(z) are only roughly equal to those in σ^. The vertical velocity dispersion nevertheless behaves very much like the velocity dispersions in the orbital plane even though its values are some 15-20% smaller than those. Rotation curves like those that are plotted in Figure 8 are generated according to the following procedure. For the graph of a given rotation curve, two projected 51 rotation velocity fields at a single time step are averaged One velocity field represents the variation in v^ together. with increasing x; this is evaluated at equal intervals in X, but only for those particles that are located in the range -1 ≤ z ≤ 1. variation in -v The other velocity field represents the with increasing y; this is evaluated at equal intervals in y, but again only for thoe particles with -1 ≤ z ≤ 1. One given rotation curve in Figure 8 is therefore the average of these two velocity fields at a given time. These velocity measurements obviously do not provide rigorous major-axis rotation curves; there are simply not enough particles available to do the latter with much statistical significance. Small-numbers fluctuations are minimized simply by including the many particles that are located within unit distance of the orbital plane. In Figure 8, R is the projected distance from the rotation axis (z-axis). R increases both with increasing x and with increasing y, as explained above. At time T = 0.0 the primary has just completed its relaxation phase. The relaxation process has induced substantial rotation in the outermost regions of the primary. Very little rotation is evident in the inner regions wherein resides most of the mass of particles. At this time 93% of the total mass in particles is located in the interval -0.5 ≤ R ≤ 0.5. For comparison, 83% and 84% of the total mass is located in that same interval at times T ≈ 4.0 and 6.0, respectively. At T = 4.0 the system of 52 particles has substantial rotation. There may be a slight turnover in the rotation curve at R ≈ 0.8, but its significance is probably minimal since less than 2% of the test particles actually contribute to the calculation of the velocity curve at all of the points beyond ∣R∣ = 0.8. The rotation velocity is nearly independent of radius for values R > 0.1 and has a mean value <vrofc> - 0.7. For values R < -0.2, <vrot> ≈ -0.5, where vrofc seems to be slightly increasing with radial distance. A difference in the mean absolute values of vrot for R < 0 and for R > 0 is already apparent at T = 0.0. 6.0. The difference is most striking at T = When T = 6.0, <vrot> ≈ 0.75 for R > 0.2, and ιvrot> = -0.4 for -0.8 < R < -0.2. Positive values of R correspond either to x > 0 or to y > 0, as explained above, depending on which of the two separate velocity fields is responsible for a particular feature of interest in the plotted rotation curve. At times T = 4.0 and T = 6.0 the cylindrical coordinates (r,φ) of the perturber relative to the center of the primary are (0.139, 51°) and (0.076, 64°), respectively. These coordinates place the perturber at both times in the same quadrant, the one with x > 0 and y > 0. This fact makes it apparent why an asymmetry in the rotational velocity curves is seen at those two times. The asymmetry simply reflects the greater ability of the perturber to enhance the rotational velocities of the nearer particles in comparison to its smaller effect on the more distant particles in the opposite 53 quadrant. Mergers described by the present technique are really capture events, wherein the perturber does not lose its identity after coalescence. The observability of the rotation asymmetry in Figure 8 is thereby enhanced by the on-going post-merger interaction between the primary "galaxy" particle distribution and the massive point perturber. A check of all the post-merger rotation curves reveals a periodic behavior in the appearance of the asymmetry on alternating sides of the primary; the period for this is remarkably close to one time unit. The observed rotation curves at times T = 4.0, 5.0, and 6.0 are very similar to one another and are of opposite parity to the set of mutually-similar curves seen at T = 4.5, 5.5, and 6.5. Cylindrical coordinates (r,φ) of the perturber are (0.139, 51°), (0.115, -1340), (0.085, 48°), (0.064, -1280), (0.076, 64°), and (0.078, -95°) at T ≈ 4.0, 4.5, 5.0, 5.5, 6.0, and 6.5, respectively. The period of this orbital motion is about one time unit, similar to that observed for the periodic behavior in the rotation curves. The correlation between the angular coordinate of the perturber and the side of the primary that shows the enhanced absolute rotation velocity is therefore verified. It is possible that the preservation of the asymmetry over more than five crossing times also implies some form of solid-body rotation in the particle distribution. The inner "galaxy", defined by ∣R∣ ≤ 0.25, which contains 67% of the test particles, does 54 in fact have a linear rotation curve. This is indicative of at least some degree of solid-body rotation in the primary. The rotating configuration described above can be compared with models of rotating elliptical galaxies (e.g., Binney 1978). A projected rotation velocity Vp is defined from Figure 8 as the mean of all values of lvrotI f°r R <-0.3 and R > 0.3 at times T = 4.0 and 6.0; the value is V ≈ 0.6. Also obtained from Figure 8 is the value of the Sr central projected velocity dispersion at those times: Oq ≈ 1.1-1.2. Therefore, (vp∕σθ) ≈ 0.5, a value typical of those found in earlier numerical experiments (White 1979a). Figure 7 shows that for all but the least dense regions of the primary an average ellipticity ε ≈ 0.25 applies. This elliptical "galaxy" therefore occupies a place in Figure 1 of Binney (1978) among his models of rotating, isotropic, oblate spheroids; its properties are also consistent with those of the rotating, isotropic, tri-axial ellipsoids. This is surely coincidental however; a low-rotation "galaxy" with equal flattening can be obtained from the simulation of a nearly head-on collision, as will be described in a later paper. iv). The Spin Parameter The collision and subsequent merger that have been described in detail above involve a transfer of angular momentum from the binary orbit to the trajectories of the individual test particles. Not all of the initial orbital angular momentum of the binary ends up in the bound merger 55 remant however. The initial angular momentum of the binary has a value Lιnιt = 0.354. The value for the z-component of angular momentum of the primary particle distribution at the end of the relaxation phase is L z = 0.030. The angular momentum of the bound merger remnant is L = 0.283, a value 26% smaller than that of the initial configuration. Only 4% of the value for L is contributed by the common revolution of the massive perturber and primary about their mutual center-of-mass. The "missing" angular momentum ∆L » Lιnιt + Lz - L = 0.10 is carried away by the stripped, unbound particles; particle escapes number 10, so that ΔM∕M = 0.010. As a result of mass loss, the remnant has a mass M = 1.490 < M1 + M2∙ An internal binding energy ∣E∣ = 3.28 is also obtained for the final primary-particle-perturber configuration. A value is calculated from these numbers for the dimensionless spin parameter λ: λ = I» ( G2 Μ5 ∣E∣'1 )"1/2 = 0.19 = 0.4 (vp∕σ) . (27) The relationship λ ≈ 0.3(νρ/σθ) found empirically by White (1979b) is consistent with that found here. Of all the models run by White (1979b), the orbital and density parameters of his model G are closest to those of the present simulation, except that his run has unit mass ratio and the one here has M2 = 0.5. He finds near-zero values for mass and angular momentum losses, plus these values ; τmerge = 2,4' ivp∕σ(p = θ∙67' ε s θ∙3θ' and λ ≈ θ∙1θ∙ values for two of these parameters, λ and ε, are nearly The 56 identical to those found here, while the values for the merger time and (Vp∕σθ) are at least consistent with their values here (e.g., ⅛ergθ = 4.0). The degree to which MTBA is consistent with expensive N-body codes is sufficient to recommend the use of MTBA in many more studies of this type. Later papers will present some of the more detailed aspects and properties of collisions and merger remnants. 57 IV. SUMMARY AND DISCUSSION The "multiple three-body" algorithm (MTBA) that is discussed in this paper was developed for the study of interacting binary galaxies. In an N-particle simulation MTBA requires 2N force calculations at each time step. For comparison, a restricted three-body integration like that of TT requires 2N + 1 calculations per time step, and traditional N-body codes require N(N-1)∕2 calculations per step. This demonstrates one of the major advantages of MTBA: its efficiency. MTBA also provides a reasonable approximation to physical reality. For example, the tidal distortions in the primary mass distribution are communicated back to the perturber, thereby forcing the secondary to respond to the damage it renders to the primary. The program does not yet also include the response of the individual test particles to their own perturbed distribution. MS have emphasized the importance of self­ gravity in these kinds of studies and it is hoped that some such modifications to MTBA can be made. To this same end, Villumsen (1982) has reportedly used Tesseral Harmonics in an expansion of the three-dimensional gravitational potential for the study of interacting galaxies (see also Clutton-Brock 1972a,b). Several of the results that are described in Sill concur with the results of other investigators. The merger times of W78 are reproducible, as shown in Table 3. In addition, Figures 5 and 6 confirm the conclusions reached by 58 MS when they pointed out that even the core of a "galaxy" is not immune to tidal effects; the original central surface density peak that distinguished the core of the pre­ collision "galaxy" was not visible in the remnant. Richstone (1975) and W78 found the opposite to be true. White reported that in nearly all of his test runs the core became more tightly bound following coalescence, while Richstone suggested that his determination of a rapid evolution of the core radius to smaller values may be no more than an artifact of his fitting procedure. MTBA also has a few problems here; the core is only partially included in the initial model (see Sllg) and the soft stellar distribution within the perturber is ignored. This paper points to additional agreement with White (1979b), in the form of the remnant spatial mass distribution; an inverse- cube radial dependence is found in both studies. Some of the implications of this latter result are discussed later. The work of W78 and others can be criticized because of their use of symmetric initial models (i.e., configurations that would be in equilibrium only at infinite binary separation). Section IIh describes how it is possible to avoid the transient effects that arise through the use of such models; a "relaxation phase" is included in the evolutionary sequence of each model binary. Because the perturber is forced to follow a Keplerian trajectory for a pre-specified time, the particle mass, energy, and angular momentum distributions within the primary "galaxy" can 59 arrange themselves in such a way as to make accommodation for the nearness of the secondary. The duration of the phase, labeled τreιaχ, is typically one dynamical timescale for the binary (i.e., the violent relaxation timescale). Such a choice for Trelaχ is justified by a study of the torque curves (Figures 3 and 4). The importance of such a relaxation phase in binary studies is shown Figure 2: merger times T1 (without relaxation) are like the merger times of W78; the more physical times T2 (with relaxation) are 50-100% larger. T1 and T2 are most discrepant in the more strongly interacting cases. Other investigators have therefore under-estimated merger times for binaries with initially strong interactions simply because their initial "galaxy" models were symmetric. Merger times are given in Tables 1, 2, and 3 for a variety of initial orbital configurations. Table 1 demonstrates the strong dependence of merger times on initial binary separation and, for the stronger interactions, the weak dependence on perturber mass. Table 2 demonstrates the reproducibility of merger times for collisions with different initial conditions. Table 3 demonstrates how the merger times of W78 are reproducible. His calculations did not include a relaxation phase and, appropriately, his times agree most with those for models that were run with 1,rθ^aχ = 0. These are underestimates of the more nearly physical merger times calculated with non­ zero relaxation. 60 White and Sharp (1977) predict that merger events will occur within one initial orbital time for binaries with separations less than or of order 3R^; R^ is the galactic half-mass radius. Their prediction can be corrected using the results of this investigation. From equation (20) it is seen that the simulations specified in the first row of Table 1 all have initial separation a of order 3R∣l∙ In units of their individual orbital periods, these configurations have merger times τιnerge∕p ~ 2.2, 1.5, 2.0, and 2.3 for M2 = 0.1, 0.25, 0.5, and 0.75, respectively. The first two binary configurations of Table 3 also have a ≤ 3R^. The merger times that are listed for their runs with non-zero relaxation are T_„_ merge ∕P = 2.4 and 1.2. The simple mean of all of the above ratios is Tmerge ∕ ' P (28) 1.9 , double the estimate of White and Sharp (1977). This decreased predicted merger rate can partly diminish the discrepancy between theoretical and empirical merger timescales, as discussed by Rood and Leir (1979). Hoessel (1980) has used the magnitude and core radius data for a sample of brightest cluster galaxies to estimate an empirical lifetime (i.e., time before coalescence) for multiple-nuclei systems; he found a time of order 3x10 years, which indeed approximates binary orbital times. 8 The disquietingly short merger timescales that have been determined by many investigators are thus corrected upwards 61 in this paper, and these new results differ little from current empirical estimates. A brief summary of the main attractive features of MTBA includes: (i) efficiency; (ii) the explicit interaction between M2 and the test particles; (iii) the relaxation phase; (iv) internal reproducibility of results; (v) reproducibility of N-body results; and (vi) calculation of realistic merger times. Of interest in the study of galaxy collisions in clusters are (i) the effect of rotation on merger rates and (ii) the amount of rotation induced by merger events. In SlIIc it was shown that "galaxy" pairs that initially rotate synchronously with their common orbital motion are no more likely to merge than initially non-rotating "galaxies". This conclusion depends upon the inclusion of the relaxation phase in the simulations and therefore differs from that of other researchers who did not include such a phase (e.g., White 1979b). More work is needed to settle this issue; for example, a self-gravitating model that includes some mechanism for relaxation should be developed. It is especially important to obtain the results of this rotationcapture problem so as to predict the probable rotational properties of the satellites of giant galaxies in clusters. Do they prefer retrograde rotation? Are satellites with direct rotation captured most efficiently, as suggested by some studies (White 1979b; Keenan and Innanen 1975)? The suggestion here is that the capture rate is independent of 62 the direction of rotation and, subsequently, satellites should not have any preferred rotation properties. Collisions produce some rotation in galaxies. For example, the remnant of the specific collision that was studied in SlIIe follows the rotation-ellipticity relation found by White (1979b). A conclusion drawn from this example is that such merger remnants will occupy those regions of the theoretical diagrams of Binney (1978) that represent the properties of the rotation-supported oblate spheroids. If it is to be believed that some reasonable fraction of elliptical galaxies are merger remnants (see Toomre 1977; compare with Marchant and Shapiro 1977), then the present data (e.g., Illingworth 1977; Schechter and Gunn 1979) require that those collisions that lead to merger events be of low initial orbital angular momentum. It should be noted that MS found very little remnant rotation. In their simulations particles that exited the studied space through a face of the bounding cube were "banned" from the rest of the calculation. These lost particles may ultimately have carried away a large fraction of the initial orbital angular momentum. It is also possible that the difference between the results of MS (little rotation) and those given here (substantial rotation) is in part due to their preference for studying unbound or just-bound pairs of relatively low angular momentum; the preference here is for bound pairs of maximal angular momentum. 63 Finally, of some interest is the form of the spatial density law in the remnant. Since surface density decreases nearly as the inverse square of the projected radius (see Figure 7), a rough approximation would lead to an inverse cube radial dependence in the spatial density. This is the same radial density behavior as found by White (1979b), who used a variety of different initial particle density configurations. Perhaps the lack of dependence of the final density law in a remnant on the initial binary galaxy parameters is a symptom of the same peculiar property of scale-free gravitational dynamics that has appeared elsewhere in astronomy.. For example, in the study of galaxy correlation functions, it is found that the predicted form of the two-point correlation is practically independent of the initial parameters that describe the Universe (Gott, Turner, and Aarseth 1979; see, however, Efstathiou 1979). As another example, Smith et al. (1979) have found that the predicted radial density and velocity-dispersion relations in model clusters of "galaxies" are independent of the spatial distribution of the unseen mass within the cluster. In conclusion, surface brightness distributions within binary galaxies may indicate their current level of interaction, but such data may be useless in a determination of the past degree of interaction unless it is supplemented by detailed velocity data. It may be that the observational problem of identifying the remnants of past interacting 64 binaries will go only partially solved even with both kinds of data. 65 APPENDIX A THE ENERGY INTEGRAL Let $ denote a velocity, so that - ^1) represents the peculiar velocity of particle i relative to the primary. Then the scalar product of both sides of equation (1a) with followed by a sum over all i, yields d$. d$. dφ∙1 m i ∑ *i ’ 3t^ ^ mi ħ * ∑ Ht~ = ^ mi M1 ∑ -3F^ where ' mi M2 I *i ∙ Viφi2 ÷ mi M2 5>1 . [ Viφi2 , ( A1 ) dφ∙1 3φ∙1 √ ■ 3^ ∙ <V ‰>∙ViΦi, - <*i' *,>∙V,Φi, ∙ <*2> The second equality in equation (A2) follows from the use of the unperturbed field of integration. throughout the numerical Similarly, the use of the unperturbed field of M2 throughout the calculation leads to d*. d«, ∑ ar^= dφiιi ■ mi M1 _∑ -⅛ dφ∙2 d$2 mi m2 ∑ ^"3t~ ^ m2 *2 , Ht” d$5 m2 ħ * 3t~ ' mi ∑ si ∙ aτ- - ”i Φ, (A3) 66 where equation (1b) has been used to rewrite the last two terms. Equation (lc) can be used to obtain at t 5 π,i ∑ vi + Ï miv1 + 5 m2v2 * mi ∑ < M1Φi 1÷ M2*i2, ! mi or d?. ∑ cit~ S1 0 , (A4) mi∏ivi+ Mil + w2φi2) + ⅛ m!v1 + ⅛ m2v2 ^ mi ∕ ( d$. ’ I 3χi ) dt constant (A5) After an integration by parts, the integral term becomes d$· J ( *1 , Σ 3t~ ) dt = $1 ∙ Σ d^1 ^ f ( 3t- ‘ Σ Equation (A5) is therefore equivalent to ) dt ∙ <a6) - mi Σ [ ⅛ < *i- ⅛,)2 ♦ M,φil ÷ M2φi2 ] ÷ 4 M2v2 + Mi ∫ [ d^1 1 ∙ ( jj I ) ] dt = constant . (A7) The energy parameter E is defined by E = ∏>i ∑t i<Vħ>2+ M1*i,t M2*i2 1 + ΐ M1v? + ⅛ m2v2 ∙ <A8) Equation (A7) then permits this statement: and only if = <$·>, where = [ ^∕N. E is constant if 67 APPENDIX B POSTSCRIPT ON NUMERICAL PROCEDURES As stated in Slle, a fifth-order predictor-corrector scheme was used for the integration of equations (1a-c). Various modifications of the basic integration formulae were suggested by Lapidus and Seinfeld (1971), but these did not significantly improve the numerical accuracy. The basic sequence PECE (or Predict-Evaluate-Correct-Evaluate) was checked for numerical accuracy with several values for h (the time step size). Since an isolated primary "galaxy" (i.e., M2 = 0) has a total system energy that is conserved exactly according to equation (4), the deviations from a constant value for the energy parameter E in this case determine the accuracy of the numerical model for a given step size. Table 4 shows the resultant percent variations in E for three different values for h and with varying numbers of integration steps. Because of the competition between efficiency (i.e., number of steps per unit time) and numerical accuracy, a step size h = 1/32 was finally chosen as a good compromise value. Predictor-corrector schemes require an initialization procedure by which the values of the numerical solution at the first few time steps are obtained. One of the Runge-Kutta-Nystrom procedures of Fehlbert (see Bettis 1973) was used to calculate these first few values from the set of values given at time t = 0. 68 TABLE 1 Merger Times for Circular Orbits MASS OF PERTURBER SEPARATION 0.10 0.25 0.50 0.75 75 8.5 5.5 6.5 7.0 10 24.0a 8.5a 8.5 9.5 7.5 22.0a 19.0a 16.0 16.5 16 37.5a 25.5a 26.0a 21 00 60.0a 25 50 >ιooa 50 75 a M2a^3<0.20 1.00 2.00 69 TABLE 2 & The Effect of Initial Conditions on Merger Times Run No. Notes^,,c τrelax τ d merge 1 seed=421, synchronous rotation 0 1.3 2 seed=437 0 2.2 3 seed=717 0 1.9 4 seed=717, N=2000 0 2.1 5 seed=813 0 2.2 6 seed=17 0.2 2.1 7 seed=817 0.5 2.2 8 seed=117 1.0 2.1 9 seed=117, synchronous rotation 1.0 1.9 10 seed=117, N=2000 1.0 2.0 11 seed=117, Tmix . =5 1.0 1.9 12 1.0 2.2 1.0 1.3 14 seed=217,’ Tmix . =5 seed=317, T . =5 mix seed=917 1.5 2.6 IS seed=17 2.0 2.6 16 seed=113 3.0 2.8 13 aThe initial binary orbital parameters were identical in all test runs: Μ^Ο.δ , a=1.5 , and an initially circular orbit. ^Unless otherwise specified, T . =2. r mix cThe "seed" is the initial number in the sequence of random numbers that determines the positions and velocities of the test particles at time t=0. d . and T are measured in units of the initial binary relax merge orbital period (P=9.42). T 70 TABLE 3 Merger Times for Non-Circular Orbits 1,W a merge M2 a e r peri rinit P τrelax 0.1 0.50 1/4 0.375 0.500 2.1 0 2.0 4.0 5.0 4.6 0.5 0.75 1/3 0.500 0.875 3.3 0 3.5 3.0 4.0 2.9 0.5 1.50b 1/3 1.000 1.750 9.4 0 9.5 16.5 12.0 0.5 - 1. 0.500 3.000 - 0 24.0C 10.0 22.5C 9.661 T merge 1.0 1.00 1/4 0.750 0.900 4.4 0 4.5 4.0 10.5 5.0 1.0 1.00 1/2 0.500 1.125 4.4 0 4.5 3.0 6.0 3.9 1.0 - 1 0.500 1.500 - 0 20.5 20.5C »4 2.5 6.735 1.0 - 1 0 4.5 - 3.000 0 4.0 2.5c 9.482 8.0 3.0c 1.500 0 aMerger times determined by Simon White with an N-body algorithm. bM2a'3<0.20 . Merger times measured from the moment when the perturber first crosses the boundary r = 1.5 . 71 TABLE 4 Numerical Accuracy versus Step Size NO. OF STEPS TIME ELAPSED ACCURACY (ΔE∕E) 1/96 384 4.0 0.002% 1/32 384 12.0 0.06% 1/32 960 30.0 0.6% 1/24 288 12.0 0.6% STEP SIZE (h) 72 REFERENCES Aarseth, S. J., and Binney, J. 1978, M.N.R.A.S., 185, 227. Aarseth, S. J., and Fall, S. M. 1980, Ap. J., 236, 43. Alladin, S. M. 1965, Ap. J., 141, 768. Alladin, S. M., and Parthasarathy, M. 1978, M.N.R.A.S., 184 871 . Austin, T. B., and Peach, J. V. 1974, M.N.R.A.S., 168, 591. Bettis, D. G. 1973, Cel. Mech., 8, 229. Biermann, P. 1976, Astr. Ap., 48, 295. Biermann, P., and Silk, J. 1976, Astr. Ap., 48, 287. Binney, J. 1977, M.N.R.A.S., 181, 735. ________________ 1978, M.N.R.A.S., 183, 501. Blandford, R. D., and Icke, V. 1978, M.N.R.A.S., 185, 527. Borne, R. D. 1979, B.A.A.S., 11, 673. Byrd, G. G. 1978, Ap. J., 26, 70. Clutton-Brock, M. 1972a, Ap. Space Sei., 16, 101. ___________________________ 1972b, Ap. Space Sei♦, 17, 292. Combes, F. 1978, Astr■ Ap., 65, 47. Da Costa, L. N., and Knobloch, E. 1979, Ap. J., 230, 639. Dressier, A. 1978, Ap. Jx, 226, 55. 1979, Ap. Jr, 231 , 659 1980, Ap. Jλ, 236, 351 Efstathiou, G. 1979, M.N.R.A.S., 187, 117. Faber, S. M. 1973, Ap. Ji, 179, 423. Gallagher, J. S., and Ostriker, J. P. 1972, A. J., 77, 288. Gott, J. R., Turner, E. L., and Aarseth, S. J. 1979, Ap. J. 234, 13. 73 Gutowski, W. J., and Larson, R. B. 1976, Pub. A. S. P., 88, 374. Hausman, M. A., and Ostriker, J. P. 1978, Ap. J., 224, 320. Hickson, P., Richstone, D. , and Turner, E. 1977, Ap. J. , 213, 323. Hoessel, J. G. 1980, Ap. Ji, 241, 493. Holmberg, E. 1941, Ap. J., 94, 385. Illingworth, G. 1977, Ap. J. (Letters), 218, L43. Jones, B. J. T., and Efstathiou, G. 1979, M.N.R.A.S., 189, 27. Karachentsev, I. D. 1972, Comm. Special Astrophys. Obs. U.S.S.R., 7, 3. Keenan, D. W., and Innanen, K. A. 1975, A. J., 80, 290. King, I. R. 1966, Ai Jjl, 71, 64. Knobloch, E. 1978, Ap. J. Suppl♦, 38, 253. Kormendy, J. 1977, Ap. J., 218, 333. Lapidus, L., and Seinfeld, J. H. 1971, Numerical Solution of Ordinary Differential Equations (New York: Academic), pp. 28, 198. Larson, R. B., and Tinsley, B. M. 1978, Ap. J., 219, 46. Lauberts, A. 1974, Astr. Ap., 33, 231. Lynds, R., and Toomre, A. 1976, Ap. J., 209, 382. Marchant, A. B., and Shapiro, S. L. 1977, Ap. J., 215, 1. Melnick, J., and Sargent, W. L. W. 1977, Ap. J., 215, 401. Miller, R. H., and Smith, B. F. 1980, Ap. J., 235, 421. Oemler, A. 1974, Ap. Ji, 194, 1. ________________ 1976, Ap. Jr, 209, 693. 74 Ostriker, J. P., and Hausman, M. A. 1977, Ap. J. (Letters), 217, L125. Ostriker, J. P., and Tremaine, S. D. 1975, Ap. J. (Letters) 202, L113. Ostriker, J. P., and Turner, E. L. 1979, Ap. J., 234, 785. Richstone, D. 0. 1975, Ap. J>, 200, 535. _________________________ 1976, Ap. Ji, 204, 642. Rood, H. J., and Leir, A. A. 1979, Ap. J. (Letters), 231, L3. Rood, H. J., and Sastry, G. N. 1972, A. J., 77, 451. Roos, N., and Norman, C. A. 1979, Astr. Ap., 76, 75. Rose, J. A. 1979, Ap. ÇH, 231, 10. Rots, A. H. 1978, Ai Ji, 83, 219. Schechter, P., and Gunn, J. E. 1979, Ap. J., 229, 472. Sharp, N., and Jones, B. 1980, Nature, 283, 275. Smith, H., Hintzen, P., Holman, G., Oegerle, W., Scott, J., and Sofia, S. 1979, Ap. J. (Letters), 234, L97. Stockton, A., and Bertola, F. 1980, Ap. J., 235, 37. Strom, S., and Strom, K. 1978, Ap. J. (Letters), 225, L93. Tinsley, B. Μ., and Larson, R. B. 1979, M.N.R.A.S., 186, 503. Toomre, A. 1974, in Proceedings of I.A.U. Symposium 58, The Formation and Dynamics of Galaxies, ed. J. R. Shakeshaft (Dordrecht: Reidel), pp. 347-365. Toomre, A. 1977, in The Evolution of Galaxies and Stellar Populations, ed. B. M. Tinsley and R. B. Larson, (New Haven: Yale University Observatory), pp. 401-416. 75 Toomre, A., and Toomre, J. 1972, Ap. J., 178, 623. Turner, E. L., and Ostriker, J. P. 1977, Ap. J., 217, 24. Vallee, J. P., Bridle, A. H., and Wilson, A. S. 1979, B.A.A.S., 11, 631. Van Albada, T. S., and van Gorkom, J. H. 1977, Astr. Ap., 54, 121. Van der Hulst, J. M. 1979, Astr. Ap., 71, 131. Villumsen, J. 1982, M.N.R.A.S., 199, 493. White, S. D. M. 1978, M.N.R.A.S., 184, 185. ________________________ 1979a, Ap. J. (Letters), 229, L9. ________________________ 1979b, M.N.R.A.S., 189, 831. White, S. D. Μ., and Rees, M. J. 1978, M.N.R.A.S., 183, 341 White, S. D. M., and Sharp, N. A. 1977, Nature, 269, 395. White, S. D. M., and Valdes, F. 1980, M.N.R.A.S., 190, 55. Wooley, R., and Dickens, R. 1961, Roy. Obs. Bull., No. 42. Wright, A. E. 1972, M.N.R.A.S., 157, 309. Yabushita, S. 1977, M.N.R.A.S., 178, 289. 76 FIGURE 1 Number of particles in a spherical shell of the pre­ interaction "galaxy" versus radius. A shell has thickness δR = 0.04 and median radius R. "Unmixed" denotes the average particle distribution generated by the Monte Carlo algorithm at t = 0, where the average is taken over five different . initial particle configurations; an error bar on one of these points specifies one standard deviation about the mean value. "Mixed" denotes the mean distribution of particles in the isolated "galaxy", as averaged over about six internal crossing times. The theoretical curve is calculated from the density law given in equation (17a). 77 R 78 FIGURE 2 Merger time Ι,1 for an initially circular binary with an "unrelaxed" particle distribution versus the merger time T2 for the same, but "relaxed" configuration. Each point represents a fixed set of initial binary orbital parameters. The solid line identifies the relation T∙j = T2∙ Note that T1 approximates T2 in the "weak-tide" domain (open circles), but that T1 underestimates T2 in the "strong-tide" domain (filled circles). 79 001 CM H I— 80 FIGURE 3 Time-dependence of the torque applied to M2 by the N test particles relative to the binary center of mass. Each point specifies the negative value of the torque summed over 16 consecutive time steps. The left panel corresponds to a run in which Treιax = θ, the right panel to a run in which τrelax = p, τ^e s^ιa^e^ region of the right panel indicates the relaxation phase in that simulation. Arrows indicate the times of the merger events in the two runs. There are two non-negative peaks in the transient torque (i.e., during the relaxation phase); these cause the shortened merger time in the "unrelaxed" simulation. The binary parameters are M2 = 1, a = 1.5, and an initially circular orbit. TIM E 81 3∩oaoι 82 FIGURE 4 The same as Figure 3, except that M2 = 2. The principal peak in the transient torque is enhanced by about a factor of two over that in Figure 3, hence the enlarged discrepancy between the values of τmerge f°r the "relaxed" and the "unrelaxed" configurations. 83 9 8 7 6 5 4 3 2 I Ο -I ΙΟ TIME 15 20 25 84 FIGURE 5 Top views (projections into the orbital plane) of the logarithmic surface density of particles at four times during the integration·. Contours are plotted at one- magnitude intervals, except that the outer two levels differ by two magnitudes. Dashed circular contours at T = 2.0 indicate the theoretically-predicted locations of each "brightness" level. The numerical integration time is given in the upper left corner of each frame; the mixing phase ends at T » 2.0, the relaxation phase ends at T = 3.5, and the merger event occurs at T = 7.5. The binary orbit is initially circular, with the perturber mass and separation given in the upper right corner of the figure. The perturber is labeled "P" and is circumscribed by a circle of radius equal to its softening radius. Arrows indicate the directions of the coordinate axes. Reseau marks are plotted at coordinate intervals equal to one unit of distance. The center of each frame identifies the primary center of mass, about which the perturber revolves in a counter-clockwise direction. 85 86 FIGURE 6 The same as Figure 5, except that these are side views of the interaction. 87 88 FIGURE 7 Structural parameters of the merger remnant versus mean isophotal radius. Filled circles correspond to measurements made on the side view displayed at T ≈ 9.0 in Figure 6; open circles correspond to measurements made on the top view displayed at T ≈ 9.0 in Figure 5. Circles in the upper panel represent the radial variations in isophotal ellipticity. Circles in the lower panel represent the radial variations in isophotal surface "brightness", measured in magnitudes. The dashed line corresponds to an inverse-square radial surface "brightness" dependence. The solid line is calculated from the projection of the spatial density given in equation (17a). 89 90 FIGURE 8 Line-of-sight projected velocity measurements taken from a side view of the primary "galaxy" (Figure 6) at three times during the evolution. The interaction phase begins at T = 0.0; the merger event occurs at T = 4.0. Each panel includes both a rotation curve and the value of the central projected velocity dispersion (one-dimensional) at the specified time. Positive rotation velocities at positive projected radii correspond to counter-clockwise rotation in Figure 5 (the same as for the perturber). Note both the flatness of the rotation curves and the differences in the mean absolute rotation velocities on opposite sides of the "galaxy". 91 1.5 --------- 1-------- Γ- ι--------1------- J--- 1--------1-------- 1--------1----- ---------,--------,-------- 1--------,--- 1--- ,--- ,--- ,--- ,--- T=6.0 1.0 ...∙∙ 0.5 . •• •• •• • • •• • ••• • ∙ •• ∙ σθ = 1.10 L -1.5 •• -∣.0 --------------1---------------- ! ----------------- -0.5 ,∙ • . Ο ■ ∙ ∙ ∙ • _ 1.5 T = 4.0 rotation 1.0 > • • ∙∙ ■·· ∙ ∙ • ∙∙∙ ∙ ∙ ∙ ββ∙ ∙ ∙ ∙ _! • 0.5 Ο • -0.5 ' ∙∙∙ ∙ ∙ ∙∙∙* .∙∙∙ ∙ ∙ β∙∙∙∙ ∙∙ .∙ -∣.0 crg≈LI9 -1.5 . . . . 1 . . . . 1.5 ■ - T=0.0 • • 1.0 • ∙ ∙ .∙ ∙∙∙ -0.5 -1.0 • • 0.5 Ο 1 > -∙ .∙∙- ∙∙ • ∙∙ ∙. ∙Σ∙ ” ··«··· ..∙ ............... ∙ • ∙∙ cΓq = 0.85 • 1 -1.5 -1.0 1 I 1 1 ■ 1 1 1 —1—1—1—1—!1 Ο 0.5 PROJECTED RADIUS -0.5 1 . 1» ∣.0 92 CHAPTER 2 MATCHING MODELS TO OBSERVATIONS To be submitted to The Astrophysical Journal Authors: Kirk D. Borne and John G. Hoessel 93 I. INTRODUCTION The study of interacting galaxies has reached the stage where such systems can be simulated realistically with reasonably simple numerical models. Comparisons between the results of those models and the data from observed binaries can provide considerable information about the evolution of pairs of galaxies. One particular numerical model and some preliminary results derived from its use are described in It was shown there that the Q uncomfortably short merger times (e.g., 10 years) found by detail in Paper I (Borne 1982). other investigators are probably correct to within a factor two. This timescale has suggested to White and Sharp (1977) and to Toomre (1977) that many pairs in the observed sample of double galaxies will merge in the next aeon and that a considerable number of pairs have merged in the past. In order to further test this hypothesis as well as test the validity of current gravitational interaction models, detailed analyses of individual binaries are needed. Such analyses require the selection of specific values for the various parameters that define each binary within the context of the particular physical model employed. The necessary observations and model results are now being collected for a sample of interacting pairs of galaxies. This paper presents the results for two of the observed binaries. Information pertaining to data acquisition and reduction procedures can be found in SII. Section III 94 presents the observational results (both spectroscopic and photometric) that are to be compared with various numerical simulations. A discussion of the numerical model is given in SIV, with particular attention given to (i) the improvements and modifications of the model that was described in Paper I and to (ii) the philosophy and method of matching simulations to observations. The multiple three-body algorithm (hereafter, MTBA) of Paper I finds its most fruitful application in the search-and-fit procedure described in this paper (SIVf). That is, many simulations can be run cheaply in the effort to find their best match to the observations of a particular binary. Sections V and VI provide complete details on the best-fit models for the observed binaries K99 (NGCs 1587/88) and K564 (NGCs 7236/37), respectively. The paper closes with a summary of what we have learned so far and what that may suggest concerning the evolutionary status of binary galaxies (SVII). 95 II. OBSERVATIONS AND REDUCTIONS a). Spectroscopy Long-slit digital spectra have been obtained for several pairs of elliptical galaxies using the high-gain video spectrometer (HGVS) on the KPNO 4-m telescope. This instrument has been described by Kormendy and Illingworth (1982). The SIT detector is read out in a 128x512 pixel configuration, with the long dimension parallel to the dispersion. The data discussed here were generally obtained by placing the length of the slit along the line connecting the centers of the two galaxies. Individual exposures were typically 50 minutes in duration, preceded and followed by exposures of the He-Ne-Ar comparison lamp. Dark frames of the SIT target (i.e., with the transfer lens between the SIT and the three-stage image tube closed) were obtained every 2-3 hours. Exposures of the dark image tube, with the transfer lens open, were made at beginning and end of each night, as were long exposures of the quartz "flat-field" calibration source. Additionally, once per night, blank sky was observed for several minutes. The latter exposure was used to determine the large-scale slit response function. All galaxies were observed at a single fixed grating o setting that covered the wavelength range 3875-4590 A. Each of the objects was selected on the basis of photometric evidence for interaction. Nearly all elliptical-elliptical pairs in the catalog of Karachentsev (1972) were observed with the SIT area detector on the Palomar 1.5-m telescope 96 (Borne, to be published). Only those that had some sign of interaction (e.g., an asymmetric light distribution) were observed with the HGVS. All of these were sufficiently small in angular extent that both galaxies could be observed simultaneously with one telescope setting, while some blank sky was exposed concurrently at the ends of the slit. latter was used for "sky substraction". The Each pixel perpendicular to the dispersion corresponded approximately to 1.7 arcsec. The slit width was 2", and the FWHM of a line in the comparison arc spectrum was about four pixels, or three times the projected slit width. The standard reduction procedure begins with the subtraction of the bias level and the dark counts from each picture. This is followed by the application of the calibration for both the individual pixel responses and for the large-scale slit response (i.e., the pictures were "flattened"). A comparison arc taken immediately before (sometimes after) an object exposure is then used to map out and hence remove the SIT pincushion distortion for that particular exposure of the object. The wavelength scale is set at this time, followed by the re-binning of the data into units of log λ. The bin width equaled 99.15 km∕s. The worst residuals for the dispersion/curvature fit were of order 30 km∕s. Typical residuals were 10-20 km∕s, scattered randomly about the best-fit polynomial. The S-distortion is then removed by (i) determining the position of the object (or the brightest portion of the object) perpendicular to 97 the dispersion as a function of position along the dispersion, (ii) fitting a cubic polynomial to these points, and (iii) rewriting the entire data frame with the appropriate shifts perpendicular to the dispersion. Finally, for each data frame, the rows of blank sky at the ends of the slit are identified and added together in a 1x512 array, which is then subtracted row by row from the entire 128x512 picture. Any repeat exposures of a particular object are then summed. One late-type giant star (in the range from G8III to K3III) was observed each night. These were used as templates in the Fourier Quotient Program (or, from now on, FQP). That program determines radial velocity, velocity dispersion, and line strength variations along the slit for each observed galaxy (Sargent et al. 1977). The stars were observed through 2.5 magnitudes of neutral density and were trailed rapidly up and down the slit. This was done to prevent heavy charge buildup at any one location on the SIT target, thus avoiding the complications of beam-pulling within the SIT (Schechter and Gunn 1979). Because the target is read out at 1.6-second intervals, beam-pulling is probably not a problem for any of our observations. Comparison arcs were likewise exposed through heavy neutral densities. The reduction of the standard star observations was handled in the same manner as the galaxy observations, except that there were no corrections for S-distortion and 98 no sky subtraction was necessary. Each star spectrum was contracted to a 1x512 format for use by the FQP. The analysis of the prepared galaxy and star frames was handled by the FQP developed at KPNO by Doug McElroy, who followed the procedure develped by Sargent et al. (1977). One major adjustment to each data frame had to be made—the Call H and K region of each spectrum was effectively removed since it has been suggested that those heavily-saturated lines can adversely affect velocity dispersion determinations (Kormendy and Illingworth 1982). A short discussion of this problem can be found in section Ilia, as can a discussion of how the FQP results depend on which particular template star is used in the analysis. For the purpose of improving the S∕N, adjacent rows in the extended spectrum of a galaxy were summed in groups of two or three. b). Pictures CCD pictures have been obtained for the same interacting galaxies for which HGVS observations exist. These data were collected at the KPNO 2.1-m telescope by G. E. Danielson and one of the authors (JGH) using the Wide- Field/Planetary Camera ground-based system, which consisted of an 800x800-pixel Texas Instruments CCD with a special optical couple that yields 0.49 arcsec/pixel. Each exposure o was made in the r(6500A) pass-band of Thuan and Gunn (1976). Typically, there were two exposures per object: a short exposure (1-2 minutes) for the study of the central regions 99 of the galaxies and a longer exposure (8-10 minutes) to bring up the low surface brightness features in the outermost parts of the light distribution. Absolute photometry does not yet exist for any of It will not be possible to extract the CCD observations. such information for all of the observations because of the non-photometric conditions that prevailed on some nights. Nevertheless, our immediate interest is in the relative photometric variations across the images of the galaxies and not so much in the absolute value of the surface brightness at each point. The latter will be forthcoming in a later paper. Each exposure was reduced in the standard fashion: the bias level (i.e., "erase" frame) was subtracted from each of the pictures, each of which was then divided by the averaged "flat field" exposure for the night, followed by the subtraction of the mean sky level on the given frame. The photographs shown in this paper were obtained from the CCD frames with the KPNO IPPS and Dicomed facilities, and the contour plots were produced with a contouring program provided by D. P. Schneider at Caltech. c). Summary of Observations Two Karachentsev pairs are the subject of this first in a series of papers. These two are K99 and K564. names for K99 are NGC 1587/88, II Zw 12, and Ho 76. fainter of the two (NGC 1588) is Mk 616. Other The Other names for K564 are NGC 7236/37, Arp 169, II Zw 172, and 3C442. A 100 third component of Arp 169 is suggested to be the projection of a galaxy at a slightly different distance (Slllb). HGVS observations of K99 and K564 in November, 1981 spanned nearly four hours and three hours, respectively. CCD observations of K99 were 30, 120, and 600 seconds in duration, while for K564 they were 120 and 600 seconds. CCD data were acquired in October, 1981 in 1"-2" seeing. The 101 III. RESULTS FROM THE DATA ANALYSIS a). K99 (NGC 1587/88) The photograph in Figure 1 was made from the 600- second CCD exposure of K99; the photo measures 4'2 on a The larger galaxy is labeled "1" in this paper, with side. galaxy "2" located at separation 59" and position angle 79° (measured east of north). The surface brightness contours for K99 are shown in Figure 2, where the levels are separated by a constant 0.75 mag. Both Figure 1 and Figure 2 have east and west reversed with respect to the appearance on the sky. The contours depicted in Figure 2 resulted from a combination of the two longest CCD exposures for K99. The 120-second exposure was used to locate the five innermost contour levels, while the five outermost levels were derived from the 600-second exposure. From both Figure 1 and Figure 2 one can recognize that K99-2 has a tail trailing to the south. This can be seen on the Palomar Sky Survey print also. What is now easy to see from Figure 2 is the asymmetric light distribution around K99-1. north. This galaxy is noticeably distended to the This pattern and that for K99-2 are suggestive not only of the presumed tidal interaction but also of motion in opposite directions for the two galaxies. It was this set of distortions in the light distribution that led us to select K99 for detailed observation. Of course, all of the pairs that were selected for study with the HGVS display 102 some peculiarities in their photometric appearance. As described in SV, the suggestions of motion and tidal interaction in K99 are supported by the results from a number of simulations of this binary. Table 1 presents the central values of the redshift and velocity dispersion, and the uncertainties in these values, for the two components of K99, as determined by the FQP. These are tabulated for each template star observed during that run. spurious. The redshifts determined with star 3 seem It was decided then to run all of our remaining galaxy spectra through the FQP with a single star; HD 79214 (star 1) was chosen. The numbers in Table 2 were determined with this star as the template. After further experimentation, it appeared to us that there was no significant dependence of the central velocities on the presence or absence of Call H and K in the digital spectra. Nevertheless, the collective experience of previous investigators has deemed it wise to excise that region of the spectral scans prior to running them through the FQP (Schechter and Gunn 1979; Kormendy and Illingworth 1982); this advice was followed. All of the rotation and dispersion profiles described in this paper were computed without the benefit of H and K. The conventional wisdom is that the widths of these lines do not solely measure the velocity dispersion among the stars that contribute to the integrated galaxy spectrum, but the widths also measure the saturation within the lines themselves. It is also true 103 that these lines are located in a region where there is a strong gradient in the continuum intensity, a gradient that is not perfectly matched between the galaxy and template star spectra. The effective wavelengths of H and K therefore do not accurately measure the recession velocity of the galaxy. Given these concerns, the dominance of the lines in the given spectral range, and our interest in the specific variations of the rotation and dispersion velocities across each galaxy, it was clearly wiser for us to remove these lines than to contaminate the velocity measurements by their presence. The uncertainties quoted in Tables 1 and 2 are the errors estimated by the FQP in its determination of the width and phase shift of the velocity broadening function. The tabulated redshifts and dispersions are actually weighted means of their values found for the few slit positions nearest the galaxy center. The values computed at each radius by the FQP were weighted by the inverse square of the error in the redshift determination at that radius. Such weights are at some level proportional to the total number of photons in the spectral scan at that position. With this weighting scheme it was possible to assign both a realistic bulk velocity and a photometric center to each galaxy on the HGVS scans. Only the numbers in Tables 1 and 2 were derived with this particular weighting scheme; the integrity of the positional velocity information shown 104 in the graphs described below has not been compromised in any way. Figures 3 and 4 show the radial variations of mean line-of-sight velocity and velocity dispersion, respectively, as extracted by the FQP from the reduced HGVS data for K99. The velocity units are km/sec; R increases to the west (to the left in Figures 1 and 2). The center of K99-1 is estimated to be located at R ≈ 32 arcsec, with the center of K99-2 at R ≈ -27 arcsec. As alluded to in the previous paragraph, these estimates are determined from a weighted mean of the locations of the points with the smallest errors in their velocity determination. The relative velocity difference between the centers of the two galaxies is 198 ± 13 km/sec. The velocities in Figure 3 are measured relative to czθ ≈ 3600 km/sec (see Table 2). K99-1 has a velocity range near 250 km/sec, which is uncharacteriscally large for an elliptical galaxy (Schechter and Gunn 1979). Similarly, its value of vrofc∕σθ « 0.6 is beyond the expected range for an E2 galaxy (Illingworth 1981). Unless K99-2 is unusually massive for its luminosity, with (M∕L)2 - 10(M∕L)1, and has a peculiarly low velocity dispersion for its mass, it is unlikely to have had much influence on the rotation properties of R99-1, especially in the inner 15-20 arcsec where the rotation profile is steep. The simulations presented in SV for R99 assume an initial rotation velocity in "galaxy 1" consistent with these observations since there was no success in the 105 attempt to generate this degree of rotation through the application of tidal torques from "galaxy 2". It is of interest to recall from the specific example of an interacting pair presented in Paper I that a rotation curve similar to that seen in K99-1 is not unusual for the remnant galaxy of a previous merger event. K99-2 has no rotation near its center, but there is the suggestion in Figure 3 that its outer regions are moving away from the observer relative to the motion of its center. This apparent recession is seen on both sides of the galaxy, which is impossible for an isolated, undisturbed galaxy. Yet such a U-shaped velocity profile is not unusual for interacting galaxies. Of the nine Karachentsev pairs that we observed with the HGVS, the following have at least one galaxy with a U-or V-shaped rotation curve: K340, and K564 (see Figure 7). K99, K107, The following are less certain candidates for having at least one galaxy with such a curve: K143, K175, K194, and K204. The latter pairs are uncertain because of the diminished S/N in their spectra at the radii of interest. Only K594 is certain not to have such a rotation curve. From the small rotation velocities and the small relative velocity of its component galaxies, no more than 100 km∕sec, we are apparently looking at K594 from an angle close to 90° from its orbital and rotation planes. We therefore do not expect any peculiarities in the shape of its line-of-sight velocity field. K594 and the 106 other six Karachentsev pairs mentioned above will be discussed in more detail in later papers. Given the weight of evidence, it is nearly certain that a Ü- or V-shaped rotation velocity curve is a signature of the tidal interaction in close pairs of galaxies. such a velocity field possible? How is What prevents the stars from wrapping around the galaxy so that their individual velocities add destructively on both sides of the galaxy, thereby destroying this remarkable rotation profile? Why do we see so many interacting pairs with at least one galaxy at this point in its dynamical evolution? The phenomenon clearly must be long-lived relative to the crossing time in a typical galaxy, which is of order 10 years. Otherwise, of the observed pairs, only a few would be expected to be found in this dynamical state. In order to best answer the questions posed in the previous paragraph, it is important to draw attention to the usual effect of tides on opposite sides of a disturbed object. As the center-of-mass responds to the bulk gravitational field of the disturbing object, the outer regions are affected relative to that center-of-mass through the differential (i.e., tidal) gravitational field. In the present case, the outermost stars of one galaxy are being heated up as a result of the gravitational interaction with the other galaxy and are being stripped, if not "boiled", out of their parent galaxy. Those stripped stars that are moving in the same general direction as the parent center- 107 of-mass will follow it and remain bound to it, whereas those that are moving in the opposite direction will be "left behind". It is the stars that have perturbed motions opposite the bulk motion of their parent galaxy that will be least bound to that galaxy and will thus be most largely separated from the galaxy center. The result is a rotation curve in which the outermost points on both sides of the galaxy are displaced in the same sense relative to the velocity at the center of the galaxy. Because the strength of the tidal field is of the same order as the internal gravitational field in those outer regions, the dynamical timescale for the tidal interaction is of the same order as the stellar crossing time within the galaxy. What occurs then is a violent, non-adiabatic process that does not allow the motions of the stars in the outermost parts of each galaxy to settle into a regular pattern where each orbit, whether prograde or retrograde, is populated uniformly around the galaxy. It is therefore perfectly reasonable to expect that pairs that are selected for observation because of distortions in their light distribution will also demonstrate unusual, probably U-shaped, distortions in their rotation curves. To further support the interaction interpretation for these unusual curves (as opposed to suggesting that there might be a flaw in the reduction programs), many of the numerical simulations for both K99 and K564 were found to possess such rotation profiles (see SSV-VI). That the 108 physical model described here (SIV) and in Paper I can reproduce these observed curves assures us that the underlying physical phenomenon is of an ordinary gravitational nature. Figure 4 shows the variations in velocity dispersion across K99. As one would probably expect, given its larger size, K99-1 has higher values for the dispersion than does K99-2. The dispersion profile for K99-1 is relatively flat, with perhaps a slight increase toward galaxy "2". Dispersion values for K99-2 certainly increase from left to right on the plot, "heating up" towards K99-1. The latter pattern is perhaps another signature of the interaction, one that is less dependent than the rotation curve on the particular spatial orientations of the binary orbital and internal rotation vectors. Unfortunately, a positive dispersion gradient is at best a weak indicator of interaction if our sample of strongly distorted pairs are at all representative of the whole population of interacting binaries. Out of the nine observed Karachentsev pairs, only K99, K340, and K564 (see Figure 8) present a non-zero gradient in the dispersion profile of at least one of their component galaxies. Though it is not plotted here, the FQP also determines a spectral line strength index as a function of position along the slit. For K99, the line strength is largest at the centers of the two galaxies and decreases radially to a uniform value equal to about 70% of the 109 central value. It is not clear what significance this holds since the errors in the determination of the index are at least of order 10-20% of the computed value at each point. b). K564 (NGC 7236/37, 3C442, Arp 169) The photograph in Figure 5 was made from the 600second CCD picture of K564; the photo measures 412 on a side. This picture is east-west reversed, as were Figures 1 and 2. The large galaxy to the southeast is labeled "1" in this paper. This galaxy is separated from galaxy "2" (to the northwest) by 35" and is at position angle 126° relative to ”2". The smaller galaxy to the southeast of the pair is labeled K564-3. East and west are again reversed in Figure 6, which is the contour diagram generated from the CCD picture depicted in Figure 5. constant 0.75 mag. The contour level separation is a Nine levels encompass the center of K564-2, while there is one less level at the centers of both K564-1 and -3. A suggestive pattern appears in the photometric distortions of K564. As shown in Figure 6, each of the two main galaxies are distended on one side relative to the appearance of the contours on their opposite sides. distentions are opposed for the two galaxies. These Distended isophotes of this type were also seen in Figure 2 for K99. They suggest motion in the opposite directions, toward which the isophotes are seen to be more closely spaced. simulations described in SVI indeed reproduce these The 110 isophotal patterns whenever the motions of the component galaxies are away from their distended sides. Finally, Figure 6 reveals a high degree of twisting in the isophotes for K564-2. Table 2 presents the central values and uncertainties of the redshift and velocity dispersion for the three components of K564. The center of K564-3 was placed on the slit for a brief exposure following the extensive observations of K564 with the slit along the line connecting the two brightest components. The recession velocity for K564-3 differs from the values for the two main components by 700 km∕sec. Because of this difference and the nearly symmetric pattern to the distortions in the light distribution around the two main galaxies, the third component of K564 will be ignored throughout the remainder of this paper. The relative velocity difference between the two brightest components of K564 (hereafter, just "the two components of K564") is 26 + 16 km∕sec, which is just enough different from zero so as to be useful in matching models to the observations for this binary (see SVI). As seen in Table 2, there appears to be no significant difference in the central velocity dispersions for the two components. In fact, if the uncertainties were to be ignored, then K564-1, the larger galaxy and the one that has the lower central surface brightness, has the smaller of the two tabulated velocity dispersions. In SVI, our simulations identify this 111 galaxy with the larger and hence more massive of the two simulated "galaxies". This was done simply because of its larger observed size. The difference between the mass in the model for "galaxy 1" and that for "galaxy 2" amounts to only 15% (see SVI). Figures 7 and 8 depict the positional dependence of the mean line-of-sight velocity and velocity dispersion, respectively, as determined by the FQP from the reduced HGVS data for K564. are km/sec. As in Figures 3 and 4, the velocity units Unlike Figures 3 and 4, R increases to the east for Figures 7 and 8. The right-most data point in both Figures 7 and 8 is our measurement nearest to the center of The centers of K564-1 and -2 are estimated to be K564-3. located at R = 17 arcsec and R ≈ -18 arcsec, respectively. The total range of velocities in Figure 7 for K564 is less than was seen in Figure 3 for K99. Velocities in Figure 7 are measured relative to czθ ≈ 8000 km/sec (see Table 2). A positive rotation-velocity gradient is evident However, the trend in near the centers of both galaxies. the points to the left of the center of galaxy "2" corresponds to a negative gradient. Given the number of data points involved in the trend, this turnaround in the rotation curve is probably real. The corresponding V-shaped rotation profile that this suggests for K564-2 was alluded to earlier in the discussion of the U-shaped profile determined for K99-2. This pattern is again considered to be a signature of the interaction. The apparent lack of 1 12 such a signature in the velocity profile for K564-1 may be due to its different dynamical orientation relative to the That the flattening of K564-1 is along the slit, observer. whereas the flattening of K564-2 is orthogonal to the slit, could be of some significance. But our preference is to suggest that there really is such a turnaround in the rotation curve for K564-1 and it will only be seen when higher S/N spectra of the outermost regions of this galaxy are run through the FQP. In fact, the right-most point for K564-1 (not the point for K564-3) is consistent with this suggestion. Most of the simulations described in §VI demonstrate the physical reality of these turnarounds in the rotation curves. Figure 8 presents the velocity dispersion variations across K564. There may be a slight tendency in each galaxy for the dispersion values to increase toward the companion galaxy. This was alluded to in the discussion for K99-2, whose dispersion profile displays a significant gradient. Again, this may be due to increased tidal "heating" of the stars in the region between the two galaxies where the gravitational interaction is strongest. For the most part, Figure 8 suggests that the line-of-sight velocity dispersion values have very little radial dependence. As is the case for K99, radial variations in the line strength index for K564 are not shown here even though they do constitute part of the output stream from the FQP. The values for the index are somewhat noisier for K564 than 1 13 the 10-20% uncertainties that were referred to in the discussion for K99. What can be gleaned from the numbers is the somewhat similar radial dependences of the index for the galaxies in both K99 and K564, though the radial variations seem to be less peaked at the centers of the K564 galaxies than at the centers of the K99 galaxies (see the last paragraph of SlIIa). 114 IV. THE PHYSICAL MODEL In this section are described the main attributes of the numerical technique that was employed in our binary "galaxy" simulations. Section IVa consists of both a brief explanation of the algorithm that was used for Paper I and a description of its most significant modifications for use in this paper. Section IVb shows how the perturbing "galaxy" is just a scale model of the primary "galaxy", while the methods by which the model "galaxies" can be flattened and by which they can receive varying degrees of internal rotation are detailed in SlVc and SlVd, respectively. Section IVe provides a description of what can be "seen" whenever the three-dimensional particle configuration that comprises the binary system is projected onto the "sky" of a randomly-oriented observer. Examples of the various observables (e.g., contour diagrams, rotation curves, and dispersion profiles) are included with the discussion. The last section (SlVf) is reserved both for describing the free parameters of the problem and for explaining how one begins to search that parameter space for a simulation that matches all of the observations of a given pair of galaxies. a). The Multiple Three-Body Algorithm (MTBA) The basic equations of motion for MTBA comprise equations (1) of Paper I. A significant modification to those equations now allows for a distribution of test particles within the perturbing "galaxy". For the 1 15 simulations described in this paper, the equations of MTBA are : d*i m i at" z M 2 at- = ■ mi M1 7i φi, - mi M2 vi φ.2 , (1a) = d$ μ at ( 1b) dħ M 1 at" = - m2 sκ ” μ 7r Φ , and d*2 aτ^ ∙ ( 1c ) The first and last equations are the same as those of Paper I. Their meanings are likewise identical to those of that paper: each particle responds only to the global, unperturbed fields of M1 and M2 (equation 1a) and global linear momentum conservation is maintained by requiring that M1i1 + M2^2 = 0 (eQuation 1<=). Equation (1b) differs from that of Paper I in the way that the force on M2 is calculated. This change is necessitated by the addition N2 test particles to the "galaxy" M2, where N2 » (M2∕M1)N1∙ most, N1 = 1000. At In equation (1b), μ is the reduced mass of the binary system, $ = ^2~ vector, and r ≈ |?2~ is the relative velocity } is the separation. In essence, equation (1b) is that of an equivalent one-body problem with Φ = ι∏i ∑ φ∣1 + m∣ ∑ φ∣2 - Φseif-gravity * The following approximation is used for Φ Φ = mi ( .Σ <t>i1 iεM, + .Σ Φi2 ) ιεM. (3) 116 Equation (1b) can then be written as M2 M1+ M2 M1 F M.+ M 2 12 f21 ‡ til m∙ M1 21 (4a) where [ 7 ∙φ∙ι iεM2 and + (il (4c ) 12 F2i (4b) is the force acting on M2 as determined by the sum of the external forces (from M1) acting on the test particles comprising M2∙ F12 is the force acting on M1 as determined by the sum of the external forces (from M2) acting on the particles comprising M1. The negative of F^2 is thus the force of reaction on M2 (Newton’s third). Equation (4a) thus represents the total force acting on M2 as the weighted mean of F21 and -F12, where the weights go as the fraction of the total mass of test particles involved in the calculation of each force term. This weighting scheme simply allows each set of particles equal significance in its influence on the dynamical evolution of the binary. If, out of ignorance, the weights were allowed to equal unity, then the net force between M1 and M2 would have been counted twice, though in two different ways (i.e., with two different sets of particles). Equation (4a) is therefore used as the best possible expression of equation (1b) within the context of MTBA. Again, it is derived simply from the 117 equivalent one-body problem for the reduced mass ∙μ = M-jM2∕ (M1+ M2). A simple manipulation of equations (1) through (4) determines the energy integral. ∙∣ μ v2 + ∑ mi ( The result is: v2 + M∙jφ^1 + M2Φi2 ) where v denotes a speed and v = |$2~ b). = constant . (5) Î∙ The Galaxy Potential The analytic form of the potential Φ∣∙1, derived in Paper I, is essentially that of an isothermal sphere with an energy cutoff. This provides a physically realistic representation for model "galaxies". "galaxy" mass, If denotes a its central velocity dispersion, R^ its cutoff radius, and aQ the softening parameter in the P expression for its gravitational potential, then the following expressions provide the recipe for building "galaxies" of different scale: G = M1 = R1 = 1 , (6a) M2 R22 = M1 R^2 , (6b) «2 a'2 . M1 a;2 , (6c) -4 *4 M2 σ2 ≈ M.J σ .J , φi2 (r) = (σ2∕σ1) 2 (6d) and φi1(a1r∕a2) , ( 6e ) or φi2 (r) = (M2∕M1) (R1∕R2) φ∙1(R1r∕R2) . (6f ) 1 18 The expression for φ∣1 can be found in equation (15) of Paper I. As in that paper, a value a1 = 0.05 is used in the calculation of the potential. Therefore, a value for M2 is sufficient to define all of the properties of that "galaxy" with respect to the primary "galaxy" M-j. c). Flattening Unlike for the model "galaxies" of Paper I, MTBA will now permit the flattening of either M-∣ or ∏2 (or both) by an arbitrary, though uniform, factor along an arbitrary minor-axis direction. If ε is the usual flattening parameter (i.e., ellipticity) for a "galaxy" and if ΰ is the unit vector along the minor axis, then the position and velocity vectors ⅞i of each particle in the initially spherical, isotropic configuration for the given "galaxy" are transformed according to: ⅞.new β §.old _ ε ( 12ε . ⅞.old } ' (7) For example, an E7 "galaxy" is obtained from the initially spherical particle distribution by means of the coordinate transformation in equation (7) with a value ε = 0.7. Each No "galaxy" is transformed independent of the other. corresponding change in the analytic form of the potential is made for any of the flattened configurations. d). Rotation Paper I allowed a rigid rotation to be added to the random motions of the test particles comprising . The 119 angular rate and direction of the rotation were restricted to that of the binary orbit. It is now possible to choose any "rigid rotation" rate and any direction for the rotation. In addition, two new modes of rotation can now be selected, if any rotation is desired for the "galaxy" at all. The last rotation mode to be described below seems to work the best in simulating the galaxies observed for this paper. The first of the new rotation modes is the "circular rotation" mode. The random motion of each particle can be replaced by circular motion at the initial radius for that particle. This motion is around a randomly selected rotation axis that is unique to that particle. The rotation axes for all of the particles are constrained to point into the same half-space. That is, if denotes a unit vector along the net rotation axis for the entire "galaxy" and if β. ≈ if X $£ defines the "rotation axis" for a particular particle, then a0 ∙ a1 > 0 (8) is required of that particle, and so on for each additional particle in that "galaxy". The second of the new rotation modes is generated with a "reflection algorithm". The random motion of each particle is left intact except that the sign of its initial velocity vector is determined by constraining equation (8) to be satisfied. This algorithm in no way changes the 120 initial particle distribution function f(E); the only problem is that the model is no longer an isotropic one, as was originally assumed in building the "galaxy". The form of f(E) does in fact change whenever the "rigid rotation" or "circular rotation" mode is selected. The complications imposed on the physical model by rotating the "galaxies" have been ignored. Since the "reflection algorithm" is the one algorithm that leaves the original physical model most nearly intact, it is both expected and found to be true that rotating "galaxies" built with that approach will most resemble real galaxies. Each "galaxy" can rotate independently of the other. The "no rotation" option can be selected for either or both of the "galaxies". the minor axis. The rotation axis need not coincide with Obviously, both retrograde and prograde rotations are permissible. It is even possible to specify what fraction of the total number of particles in a given "galaxy" will participate in the rotation. With so many rotation options to choose from, a dynamical configuration for the particles can usually be found that matches the observed inner rotation curve for a given galaxy. The observed rotation velocities in the disturbed outer regions of the binary galaxies are generally matched best following a period of tidal interaction between the two model "galaxies". 121 e). Projections In order to facilitate the matching of models to data, MTBA has been set up to run interactively on the University of Michigan Amdahl machine. An interactive graphics package on the system is used to plot projected mass distributions and velocity curves (rotation and dispersion) from the simulations. The "observer" may "look" at the simulated interacting binary from any direction in space and draw a smoothed contour plot, calculated from the projected positions of less than 2000 test particles, or may look at the projected positions of the particles themselves, or may draw a rotation velocity (or dispersion) profile by choosing a slit width, a slit position angle, and a position Any number or combination of for the mid-point of the slit. such diagrams can be plotted at predetermined intervals during the numerical simulation. The "pixel" size chosen both for the contour plots and for the velocity curves is 0.05 in units where = 1. This size equals the value of the softening parameter a1 for "galaxy" M^. Contour levels are found after the application of a 7x7 "triangular" smoothing function to the raw projected mass distribution. A value for the "surface brightness" is calculated at each pixel. A dot is drawn to the plot if the integer part of the "surface brightness" is an even number, and nothing is done if the integer part is odd. In this way each drawn contour level is one magnitude wide and is separated in value from the next drawn contour 122 by two magnitudes per square pixel. The un-drawn contour level between these two is also one magnitude wide. As an example, a specific initial particle configuration is shown in Figure 9, with the corresponding contours depicted in Figure 10. In these figures, N1= 1000 and N2= 200. The scale of Figure 9 is 50% larger than that of Figure 10 in order to show the individual test particles more distinctly. Figure 9 is 4.0 distance units on a side and Figure 10 is 6.0 units on a side. Velocity curves are calculated with a 5-pixel triangular smoothing function applied to the raw projected velocity profiles. Because of the small numbers of particles involved in these calculations, slit widths in the range from 0.5 to 1.0 are typical. Each particle within the slit boundaries is assigned to a particular bin (i.e., a column of pixels) along the long axis of the slit independent of how far that particle is above or below the central line defining the slit position. In this way, with a wide slit, small-numbers fluctuations in the simulated velocity curves are kept to a minimum. A small contour plot is generated and presented simultaneously with each velocity curve. This plot not only shows how the "galaxies" appear from the particular direction of the observer, but the plot also includes a drawing of the slit boundaries relative to the contours. These particular plotting capabilities and other features of the program render it remarkably flexible and user-friendly. 123 A rotation curve for the model "galaxies" depicted in Figures 9 and 10 is shown in Figure 11a for an observer at an oblique angle to the system. The corresponding velocity dispersion profile is shown in Figure 11b. The range of the projected-radius coordinate (i.e., the abscissa) in these figures and in all of the other simulated velocity profiles shown in this paper is -2 to 2 (units defined by equation 6a). The velocity range (i.e., the range of the ordinate) in Figure 11a and in all of the other simulated rotation curves depicted in this paper is -2 to 2 velocity units, whereas the velocity range in Figure 11b and in all of the other simulated dispersion profiles in this paper is 0 to 2. The velocity scales for the rotation curves and the dispersion profiles therefore differ by a factor 2.0 for all of the simulations presented in the paper. The total length of the slit in all cases is 4.0 distance units. If Nf is the total number of particles contributing to one point on the smoothed velocity profile and if v and σr denote the mean line-of-sight velocity and the line-of-sight velocity dispersion, respectively, then the upper and lower curves drawn in Figure 11a correspond to v + αr and vf- α , respectively, where αr = σr∕Nr,. Similarly, the upper and lower curves drawn in Figure 11b 2 correspond to σf+ βr and σf- βj,, respectively, where βr = σr∕2Nr. These curves are disconnected across those pixels where Nf < 6, for N equal to the weighted sum of particles 124 included within the smoothing triangle (i.e., Nr= Nj_2+ 2Nj.1+ 3Nj+ 2Nj+1+ Nj+2 at pixel j). The complete details of the particle configuration that led to Figures 9 through 11 will now be given. This is done simply because that configuration is precisely the initial model used for one of the simulations of K99 to be described in §V (i.e., Model 3 of that section). The initial diameter of M1 is, of course, 2R1 = 2.0 units, the mass ratio is M2∕M1 = 0.20, the separation is R = 1.265, neither "galaxy" is flattened, the secondary M2 has no net rotation, and the primary M1 is rotating prograde about an axis whose spherical coordinate direction angles are (θ,φ) = (30o,180o). The latter rotation includes all N particles and the rotation was generated according to the "reflection algorithm" described in the previous sub-section. The rotation axis for the binary system coincides with the coordinate z-axis at (0o,0o)j the initial binary separation vector is parallel to the positive x-axis at (90°,0°); the observer of Figures 9 and 10 is at (0°,0°); and the observer of Figures 11a and 11b is at (50°,90°), along a direction roughly half-way between the y- and z-axes. The latter observer does not see the full rotation of M1. Nor does that observer measure the full orbital velocities for M1 and M2. As in Paper I, the binary orbital motion is constrained to lie in the xy-plane, albeit somewhat inconsistent with the use of "galaxy" models that are no 125 longer symmetric across that plane. The binary configuration of Figures 9 through 11 is that of a system whose orbital eccentricity is 0.15, whose semi-major axis is 1.1 units, whose current separation is that corresponding to apogalacticon (i.e., R ≈ 1.265), and whose relative velocity vector is of magnitude 0.898 units, parallel to the positive y-axis. The projected velocity difference for the observer of Figures 11a and 11b is 0.688 and the projected separation equals the true separation 1.265. f). Free Parameters and the Match Game Before the best matching simulations for K99 and K564 are presented, it is useful to describe the philosophy and method of fitting these models to the given observations. The philosophy behind our model-fitting procedure is that the abundance of information carried in the data ought to place strong constraints (or restrictions) on the range of acceptable values for a fairly large number of model parameters. Given M1 ≈ R1 ■ ^l ∣ given an analytic form for the gravitational potential, and given a specific particle energy distribution function, the results for a particular simulation will depend on three sets of parameters. These are (i) the binary orbital parameters, (ii) the internal rotation and flattening parameters, and (iii) the angles that define the spatial orientations of all these with respect to the observer. The match between theory and observation ought not and does not depend on 126 which particular statistical realization of the chosen distribution function is selected (Paper I). Of the three sets of parameters, the most successful attempts to constrain the values for one of them occurred for the set of binary orbital parameters. These orbital parameters are R, M2, and the relative velocity vector, where R is the binary separation and M2 is in effect the. mass ratio since M1 = 1. Another way of identifying these parameters would list a, e, M2, and R, where e is the orbital eccentricity and a is the binary semi-major axis. Each of these parameters, except perhaps the eccentricity, is tightly constrained by the data because of the strong mass and distance dependences of the tidal interaction terms and because the velocity scale is set by the internal velocity dispersion measurements. More will be said below on the manner in which these constraints work to bring about a reasonable match to the data. The internal parameters for each galaxy are often only weakly constrained by the observations. If the galaxy has non-zero ellipticity and shows rotation along the slit, then something more can be said about the dynamical configuration for that galaxy than for one with very little rotation or with near-zero ellipticity near its center. Typically, the internal parameters for each galaxy are constrained not only by the appearance of the galaxy, but also by the ability or inability of a given model configuration to develop, over the course of a simulation, 127 bulges and tails that are similar to those that are seen in the picture data or to develop rotational properties that are similar to those seen in the HGVS data. Only certain models will develop into a configuration that matches the data. The internal parameters that can be selected for each model "galaxy" are the flattening, the angles defining the direction of the minor axis, the mode of rotation (see SlIId), the direction of the rotation axis, and the fraction of particles participating in the rotation. Something will be said about each of these below. If a galaxy is intrinsically flat, but the direction of the observer is nearly along the minor axis, then it is difficult to derive the degree of flattening from the simulations described in this paper. In fact, the intrinsic flattening is hard to determine from any angle of observation. A lower limit on its value can be set from the observed ellipticity, which is neither an original suggestion nor a particularly enlightening one. The direction of the minor axis for a "galaxy" with moderate flattening is somewhat better constrained through the behavior and development of the isophotal twists and tails over the course of the numerical simulation. The mode of rotation chosen for most of the "galaxy" models in this paper is that generated by the "reflection algorithm" (Sllld). Of the three modes of rotation described in Sllld, only that one preserves the original physically reasonable energy distribution function. The 128 dynamical development of a "galaxy" with the "circular rotation" mode is incompatible with what is now seen in K99 and K564. When the "rigid rotation" mode is employed, the initial "galaxy" model mixes so efficiently that the rotational velocity data cannot be nearly so well matched late in a simulation as they were early in the simulation. Both the direction of the rotation axis and the fraction of particles actually participating in the rotation for a particular simulation will determine the shape, gradient, and velocity range in the "observed" rotation curve. It is therefore possible to restrict these parameters fairly well in the real binaries, provided that the three-dimensional velocity distribution in the model actually mimics that for the real galaxy to some degree. The spatial variations of the line-of-sight stellar velocity dispersion (e.g., Figures 4 and 8) are not used to any great extent in the attempt to match models to the data, except that their central values are used to calibrate the velocity scale for the binary. Once calibrated, the projected relative velocity of the binary constrains the possible projection angles of the total relative velocity vector. The shapes of the velocity dispersion profiles depend on the velocity distribution function within the galaxy, about which very little is known. The distribution function used in our models is a physically reasonable one, but it is not necessarily the correct one. 129 The location of the observer with respect to the many spatial parameters of these binaries is usually well- determined. To be sure, a simulation is not pursued to any length if a good viewing angle for the observer cannot be isolated. A "good viewing angle" is one that reveals the proper photometric distortions, the correct recession velocities relative to the scale of the central velocity dispersion values, the correct shapes for the rotation curves, and a projected appearance of the "galaxies" compatible with observation. Because the observer orientation and the binary orbital parameters are well determined, it is a simple matter to calculate total masses and true separations for the galaxies in the observed binary sample. From these numbers, estimates of binary orbital periods can be made, as can estimates of the expected merger timescale for the systems. The calculation of the latter assumes that the relevant timescales determined from interacting galaxy simulations are physically sound (see Paper I). In order to isolate a reasonable matching model for a binary, a few physical insights are helpful. If one attempts to simulate the data with a model whose binary separation is too large, for example, then the resulting tidal field in the simulation will be too weak to bring about the degree of distortion that is seen in the light distribution. Yet if one suggests a separation that is too small, then the tidal distortions in the simulation may 130 greatly exceed those that are observed. In the latter case there will also be a crowding of the two "galaxies" and of their isophotes to the point that the picture data just cannot be matched. A reasonable guess can therefore be made for the value of the true separation, with the exception of those cases where the orbital eccentricity is large (i.e., where the time derivative of the tidal field is not small over the course of the integration). latter are discussed below. The effects of the To begin the long process of isolating a good simulation for a particular binary, one can initially choose to assume that the projection angle of the separation vector relative to the line-of-sight vector from the "observer" will be 90°. One can then start a simulation, look at it from such a direction, and then examine the isophotal distortions over the course of the integration. If these distortions are too large relative to the those observed in the data, then the true separation must exceed the projected separation in order that the tides not be so strong. If, on the other hand, the model is "observed" at too small a projection angle, there will again be the problem that the isophotes of the two "galaxies" will crowd one another more so than permitted by the data. According to this line of reasoning, it is possible to place tight restrictions on the range of values that are acceptable for the projection angles of the separation vector. 131 Even though the ellipticity of the orbit cannot be determined very well quantitatively, a reasonable qualitative estimate can be made. If the relative trajectory for the simulated pair is highly eccentric, then the relative velocity of the "galaxies" will be very large at perigalacticon. In this case, it may not be possible both to match the observed velocity difference and to have a good match to the other observables. Furthermore, by the time the projected separation has increased back to that observed there may be far too much tidal damage in the system. The scenario described in the preceding statements does not apply to those cases where the orbital eccentricity is moderate to small. One then must judge from the appearance of that particular simulation if the distortions, the magnitude of the projected relative velocity, and the line-of-sight velocity field of the real galaxies are all consistent with the relative orbit of the pair lying on such a trajectory. Clearly, these represent a non-trivial and somewhat subjective set of comparisons that must be made. Complicating matters is the fact that most of the distortions of interest are being mapped by the very few particles in the outer reaches of the simulated "galaxies". This is not really as severe a problem as it appears since there is in reality no particular upper limit to the number of test particles that can be used with MTBA. In principle, a good, working model can be rerun with many more test 132 particles so as to better check the finer details of the simulation. Models with circular binary orbits are usually among the first of the simulations to be investigated. If the "photometric" distortions in such a simulation are too small, then the effect of a smaller circular orbit is investigated next, and vice versa. If the actual shapes and twisting of the simulated isophotes differ from those seen in the data, then this is usually an indication that some degree of eccentricity is needed in the binary trajectory. One of the most rigid constraints that can be placed on any of the parameters for a simulation is that for the relative velocity projection factor. For example, suppose one chooses an orbital inclination for a particular simulation and the projected relative velocity is calculated. The velocity data for the real binary can then be scaled so that its projected relative velocity equals that of the simulation. Dispersion profiles for the simulation and for the real binary are then plotted on the same scale. If the velocity scale is set too high (i.e., our line of sight to the simulated pair is at too small an angle to the total relative velocity), then the scaled velocity dispersion data for the real pair will exceed the scaled dispersion values for the simulation. If the velocity scale is set too low (i.e., our line of sight intersects the relative velocity vector of the simulation at too large an angle), then the scaled dispersion data will 133 fall below the values for the simulation. It is therefore an easy matter to determine the velocity projection factor to within 10% of its actual value. Because the projections of the relative velocity vector and of the separation vector are both determinable to better than 15%, total masses for observed interacting binaries can be estimated to about the same precision. Large errors in the measured relative velocity for a binary (e.g., that for K564) will not affect the mass calculated for the system as much as it will affect the determination of the three-dimensional orientation of the pair relative to the observer. In the following two sections, a variety of models for K99 and K564 are presented. The differences between the various simulations presented for a given pair measure the uncertainties in the physical results finally adopted for that pair. 134 V. SIMULATIONS OF K99 Four simulations are presented in this section that represent some of the better matches to the data for K99. The characteristics that these have in common are described first. This is followed by a discussion of many of the various model parameters. All of the parameters that define the simulations are listed in Table 3. Specific contour plots and velocity profiles from these simulations are provided in Figures 12 through 17. Values for the parameters that describe the matching models and values that represent the adopted physical solution of K99 are presented in Tables 4 and 7, respectively. The ratio of the velocity dispersion for K99-2 to If their mass-to-light ratios are 4 nearly equal and if the L « σ relationship is assumed valid that for K99-1 is 0.67. here (Faber and Jackson 1976; Sargent et al. 1977), then the mass ratio for the K99 binary is M2∕M1 ≈ 0.20. If equation (6b) holds true, we expect Rj/R-i ≈ θ∙45, which is found to be an accurate representation of at least the intermediate isophotes for the two galaxies in Figure 2. The 1:5 mass ratio was therefore adopted for all of the simulations of K99 listed in this paper. The numbers of test particles were N1 » 1000 and N2 ss 200 for "galaxy 1" and "galaxy 2", respectively. Neither "galaxy" was ever flattened and only M1 was given some rotation. In all cases, the rotation of "galaxy 1" was generated with 100% particle participation in the "reflection algorithm" (Sllld). Only in Model 3 was the 135 rotation axis tilted from the vertical (i.e., with respect to the binary revolution axis). This was done in order to allow the "observer" to detect a larger fraction of the internal rotation in that "galaxy", which seems required by the large observed rotation rate in K99-1. Finally, all simulations began at apogalacticon and all started on trajectories with zero to moderate eccentricity. Table 3 provides a list of the initial values for the parameters defining the four models of K99. Given there are the mass ratio (M2∕M1), the semi-major axis of the orbit (a), the initial orbital period, the orbital eccentricity (e), the initial separation (R), the azimuthal component of the relative velocity (v*b, its radial component (vr), and the various internal parameters that have been discussed at length in the previous section. The angles that locate the rotation axes are those of a spherical polar coordinate system whose polar axis is identical to the axis of revolution for the binary (i.e.,the z-coordinate axis). Table 4 lists the model parameters and the projection parameters that provide the best match of each model to the K99 data: the model time (or elapsed time) at which this match occurs, the binary separation at that time, the position angle of M2 relative to the initial binary separation vector, the location of the observer (in spherical polar coordinates, where the binary revolution axis coincides with the z-axis), the orbital inclination and phase relative to the observer, the numbers for the figures 136 in this paper that depict these configurations, the projected relative velocity (∆Vproj), the projected separation (R 0j), the mass projection factor k (defined in the footnote to the table), and the estimated time remaining for the binary model to appear as a distict pair of "galaxies" (see Paper I). Figure 12 contains four contour plots, one for each configuration defined in Table 4. units where R1 = 1. Each box measures 4x4 in The enclosed rectangle represents the slit for the spectroscopic "observations" shown in Figures 13, 14, 15, and 16 for Models 1, 2, 3, and 4, respectively. The projection angles and binary orbital parameters that give rise to these plots are those listed in Table 4. The velocity data for K99 shown in these figures have been scaled so that the observed (i.e., projected) relative velocity of K99 equals that of the particular model configuration depicted. Figure 17 shows a large-scale top view of Model 3 at a time just past that shown in Figure 12c. The binary motion is counter-clockwise in Figures 12 and 17. There is a remarkably well-defined stream (or "tail") of stars leaving "galaxy 2" in the latter figure. This was seen in every simulation, though not always this well. It is clear that the distortions observed in K99 are explainable by tidal gravitational phenomena alone. The expected motion for the components of K99 in the directions opposing their distended sides is confirmed in every simulated model presented here. Though not obvious in 137 every instance depicted in Figure 12, the trailing edge of K99-1 often would be more distended than the leading side. It can be said, with little room for doubt, that such light distributions (i.e., opposing distentions for the two galaxies) will be a sure indicator of the sense of binary motion in pairs that are so disturbed. The much-discussed U-shaped rotation curve of K99-2 (see SlIIa) is not matched perfectly in any of these plots. The rise in the rotation velocity of the secondary "galaxy" toward the primary "galaxy" is seen near the center of every rotation curve depicted here, just as it is seen in the data for K99-2. The rise expected on the left-hand side of these curves is seen less frequently. This is mostly a result of the small numbers of test particles in that region of space . and is not due to the inability of the tidal interaction to produce such a rotation profile. As in Figure 13a, every simulation showed this particular type of rotation profile at one time or another; they just do not have enough particles occupying that region of phase space to allow a O- shaped curve to be seen at each time step. By comparison to Figure 13b, the velocity dispersion profiles for the simulations are seen to have risen in the region between the two "galaxies"; such a "heating" effect was noted in the data for K99-2. In Table 7 a complete array of physical properties for K99 is given. The lack of absolute photometry precludes listing a value for its mass-to-light ratio. All other 138 parameters listed there are derivable from information provided in this paper, except for the virial estimate of the mass of K99-1. For that estimate, a value of the effective radius equal to 0'.4 was assumed. Most noteworthy in this table is the precision to which the total mass for the system can be determined by the model-matching procedure. Any systematic errors in the velocity data will render such an estimate inaccurate. But, within the context of the present physical model, the errors in the estimated orientation angles for K99 are small and it seems certain that our view of the system is nearly 90° ahead of the separation vector and about 50° from the binary revolution axis. Since the galaxies are not strongly disrupted, it also appears certain that their relative trajectory is one of moderate eccentricity. 139 VI. SIMULATIONS OF K564 Three simulations are presented in this section as good matches to the data for K564. Two are very similar (Models 2 and 3); each of these does equally as well as the The other model (1) other in fitting the observations. reproduces the velocity data and the outermost contours remarkably well, yet fails to match the inner isophotes. Each of the model "galaxies" in the latter simulation developed a prolate shape, even though the initial state for each "galaxy" was a spherical one with a large rotation rate. Model 1 and its prolate components are included here to demonstrate that it is not unlikely for galaxies of this shape to develop in a natural way (i.e., in gravitational encounters) within a very short time relative to the age of the universe. The most unnatural aspect of Model 1 is the mode of rotation chosen for its component "galaxies": the "circular rotation" mode described in SlVd. All of the parameters that define the three models of K564 are given in Tables 5 and 6. Specific contour plots and velocity profiles from these simulations are depicted in Figures 18 through 22. The types of information presented in these tables and figures are identical to those described in the previous section for K99. The final table (Table 7) delineates all of the physical parameters for K564 (as well as for K99) that can be constrained to one degree or another by the information made available in this paper. 140 The ratio of the velocity dispersions for the two components of K564 is consistent with unity (see Table 2). In order to force the larger observed galaxy, K564-1, to be represented by the larger model "galaxy", a mass ratio slightly less than unity has been used in the simulations. The results of the interaction depend very little on the value for the mass ratio, unless the corresponding dispersion ratio (see equation 6d) becomes inconsistent with its allowable range from the data (see Table 2). In these simulations, the total number of test particles was always around 1300, apportioned according to the binary mass ratio. The various types of rotation and the particular rotation axes that were selected for the model "galaxies" are shown in Table 5. It- seems clear that if the non-prolate models (2 and 3) provide the proper physical representation of the dynamical state for the galaxies in K564, then the rotation rate in K564-1 is much less than is required for it to be rotationally flattened. The percentage of its particles that are contributing to the net rotation is of order 20-25%. This is in sharp contrast to the unusually high rotation rate detected in K99-1. Finally, the three simulations presented for K564 each began on a circular orbit. Figure 18 contains three contour plots, one for each configuration defined in Table 6. The enclosed rectangles represent the slits for the spectroscopic "observations" shown in Figures 19, 20, and 21 for Models 1, 2, and 3, 141 respectively. Figure 22 shows a large-scale top view of Model 3 at a time just later than that of Figure 18c. The binary motion is counter-clockwise in Figures 18 and 22. The trailing distended mass of particles seen emanating from both "galaxies" in the all of these figures is just what the data in Figure 6 reveal. As is true for K99, the tails are indeed an indicator of the trailing side of the moving galaxies. Noteworthy in Figures 20a and 21a are the sinuous rotation profiles for the models, which match the behavior of the data quite well. The rise in rotation velocity to the left in these plots is consistent with that seen in Even the rotation curve in Figure 19a matches the K564-2. data well, almost too well. As for the dispersion profiles, these too match the data over a large fraction of each "galaxy". The initial shapes of these profiles were basically the same as that seen in Figure 11b. Tidal heating of the particles between the two "galaxies" is seen in the matching models for K564 just as it was seen in the simulations of K99. From Table 6, it appears that the orientation angles of K564 are fairly well constrained. This would be a safe assessment except that these angles depend on the value of the observable velocity difference between the two galaxies. As stated in SlIIb, the error in the determination of this number is only barely less than value of the number itself. It does seem likely that the line-of-sight to K564 cannot be 142 far from the axis of binary revolution, particularly because the projected velocity difference is so small. Yet it is not clear why the virial mass estimate in Table 7 differs so much from the mass estimate derived from our simulations since the latter mass estimate (also given in Table 7) should not depend on the binary orientation angles. That is, the mass scale (i.e., the relative velocity scale) is set by the velocity dispersion measurements alone. 143 VII. ON THE EVOLUTIONARY STATUS OF BINARY GALAXIES It has been demonstrated quite convincingly that the distortions in the light distributions for at least two pairs of elliptical galaxies are reproducible by a simple gravitational model. Because gravity is scale-free, unlike the gas-dynamical processes active in spiral galaxies, the physical model described here will apply to any other pairs of ellipticals that happen to chance by one another as these have done. It has also been shown that it is possible to constrain the physical parameters that describe the binary. This is accomplished through the use of positional spectroscopic data, such as given by the HGVS at KPNO, and through the use of surface photometric observations. The latter are of greatest value when the faintest parts of the galaxies can be detected; those outer parts are most sensitive to tidal gravitational fields and are thus a good probe of such fields. Because the photometric data are derived from wide-band observations, these data are potentially of more value than the elusive narrow-band spectroscopic data. If only central redshifts, central dispersions, and good surface photometry were available, then it would be possible to apply some of the arguments of this paper (SlVf) to constrain a more limited, though still useful set of physical parameters for a much larger sample of pairs. A data base of this type, consisting of both surface photometry and central velocities, is now being compiled. Asymmetric and other variations in the light 144 distributions of nearly all of the Karachentsev E-E pairs will be presented in a later paper, as will central redshifts and dispersions for these. Also forthcoming is an analysis of the mass-to-light ratios for the nine pairs that have been or will be analyzed by the full procedure described in this paper. From the collections of well- determined binary separations and of good M/L estimates one can then at least say something more definitive than has been said on the probable extent of dark halos surrounding elliptical galaxies. Of further observational and theoretical interest is the usefulness of these simulations in distinguishing prolate from oblate galaxies. As was pointed out in the discussion of Model 1 for K564, the development of a prolate galaxy from an initially axisymmetric state is one of the possible natural outcomes of a tidal encounter. Though the following approach was not taken, it would be possible to begin a simulation with a prolate "galaxy" model, let the interaction proceed, and then see if the appearance of the resulting distorted "galaxy" is consistent with any observed system. If a good case can be found, then detailed spectroscopic observations would be in order. The latter would further constrain the internal dynamical configuration for the galaxy. The binaries described in this paper are clearly not yet merged. Yet the rotation rate of the primary galaxy in K99 is considerably more than that of the typical E2 galaxy. 145 It was shown in Paper I that such high rotation rates can be produced in a merger episode. Perhaps K99-1 is a merger remnant about to "consume" yet another companion. Whatever the origin of its rotation, further spectroscopic observations of this galaxy are clearly necessary. The use of various other slit position angles will lead to a detailed map of the velocity field in this galaxy. The tidal interaction and merger hypothesis can then be tested with these more detailed spectroscopic data. A few ideas on the evolutionary state of known pairs of galaxies can be offered even with the results for only two binaries at our immediate disposal. For example, it is impossible for pairs of galaxies on highly eccentric trajectories to avoid a great splash at their closest approach. Since very few pairs (less than 10% of the Karachentsev E-E's) are as disrupted as that, it is likely that highly elliptical orbits are under-populated at the present epoch. Of course, because of the transience of the configuration, very few pairs will be seen at perigalacticon regardless of the distribution function. But, if one of these pairs has had more than one close approach, there is a good chance that the stars on the most tidally-sensitive trajectories within their parent galaxy will have been stripped from that galaxy in the past. On the next approach, there may be less of a splash, but there will probably be a sea of background stars in which the pair of galaxies is imbedded. Such a light distribution ought to be 146 detectable and distinguishable from what is seen for a pair whose relative trajectory has been a slowly decaying spiral for many orbital periods. It is clear from these arguments and from the appearances of K99 and K564 that these two are most probably on trajectories of low eccentricity. That is, these two systems reveal well-defined bulges and tails; they do not possess large, smoothed-out background light distributions. The study of many more pairs will aid in the determination of the global distribution function for orbital eccentricities in binary galaxies. The simulations described in this paper cannot confirm the validity of the dynamical friction timescales that have been calculated to be very short (i.e., about one or two orbital periods). Because the dynamical friction hypothesis depends on simple gravitational physics, as do the simulations presented here, then the validity of one ought to lend support to the validity of the other. If the present technique is a viable one for studying interacting binary galaxies, then it is probably true that, when two galaxies approach one another on a non-hyperbolic trajectory, the outcome of their encounter is coalescence. 147 TABLE 1 FOURIER-QUOTIENT VELOCITIES FOR THE K99 PAIR c^-Zq) c Star Name Type c(z^-zθ) 1 HD 79214 KO III 105 + 10 206 ± 11 -93 ± 8 138 + 10 2 HD 83820 Kl III 115 ÷ 11 216 ± 13 -102 ± 8 146 ± 11 3 HD 214202 G8 III 51 + 11 226 ± 12 -165 + 8 149 ± 10 aczθ=3600 km/sec for K99. subscript "1" denotes the east component of K99. 0 Subscript ''2" denotes the west component. 148 TABLE 2 ADOPTED VELOCITIES AND DISPERSIONS Galaxy c(z-zθ) cz σ K99-l (East) 105 ± 10 3705 206 + 11 K99-2 (West) -93 8 3507 138 ± 10 K564-l (SE) -134 11 7866 177 ± 14 K564-2 (NW) -160 12 7840 183 ± 15 K564-3 (faint; SE) -843 32 7157 225 42 149 TABLE 3 INITIAL PARAMETERS FOR SIMULATIONS OF K99 Parameter Model 1 Model 2 Model 3 Model 4 m2∕m1 0.200 0.200 0.200 0.200 a 0.884 1.081 1.100 1.250 4.8 6.4 6.6 8.0 e 0.414 0.156 0.150 0.0 R 1.250 1.250 1.265 1.250 VΦ 0.750 0.900 0.898 0.980 0.0 0.0 0.0 0.0 M^ Rotation Mode Refl. Refl. Refl. Refl. M^ Rot. Axis (θ,φ) (0o,0o) (0o50o) (30°,180°) (0o,0o) M^ Ellipticity 0.0 0.0 0.0 0.0 M2 Rotation Mode None None None None M2 Ellipticity 0.0 0.0 0.0 0.0 p orbit V r 150 TABLE 4 PARAMETERS FOR K99 MATCHING MODELS Parameter Model 1 Model 2 Model 3 Model 4 Elapsed Time 3.5 4.5 6.5 9.0 Separation 0.991 1.055 1.264 1.061 Position Angle -51° -91θ -5° 64° Observer (θ,φ) (50°,39°) (50o,-lo) Inclination 50° 50° 55° 50° Orbital Phase -90° -90° -110° -120° Contour Plot Fig. 12a Fig. 12b Fig. 12ca Fig. 12d Velocity Plots Fig. 13 Fig. 14 Fig. 15 Fig. 16 ∆v Pr°J R proj 0.725 0.807 0.688 0.716 0.991 1.055 1.214 0.981 „ . b Projection Factor 2.30 1.75 2.09 2.39 . c Time until Merger 6 8 7 7 3. fo (55°,105°) (50o,184o' Figure 17 depicts Model 3 from (0o,0o) at an elapsed time 7.0 . Projection factor k defined by: cTime = 2∙P ,∙4.-T1 , orbit elapsed k = ( M, + M„) ∆v . 1 2 proj (see Paper I). r , 2 R . proj 151 TABLE 5 INITIAL PARAMETERS FOR SIMULATIONS OF KS64 Parameter Model 1 Model 2 Model 3 m2∕m1 0.833 0.850 0.850 a 1.50 1.50 1.35 ^orbit 8.5 8.5 7.2 e 0.0 0.0 0.0 R 1.50 1.50 1.35 VΦ 1.106 1.111 1.170 0.0 0.0 0.0 Circ. Refl. Refl. 100 20 25 M^ Rotation Axis (θ,φ) (10°,45°) (60°,90°) (75o,90o) M^ Ellipticity 0.0 0.7 0.7 (60°,90°) (75°,90°) V r Rotation Mode Participation (%) 3. M^ Minor Axis (θ,φ) M2 Rotation Mode Circ. Refl. Refl. M2 Rotation Axis (θ,φ) (10°,45°) (15°,90°) (30o,90o) M2 Ellipticity 0.0 0.0 0.0 a Percentage of particles participating in the rotation-generating algorithm 152 TABLE 6 PARAMETERS FOR K564 MATCHING MODELS Parameter Model 1 Model 2 Model 3 Elapsed Time 4.0 12.5 10.5 Separation 1.500 1.337 1.107 Position Angle 169° 178° 169° Observer (θ,Φ) (3?6,259°) (4?0,268°) (3?7,259°) Inclination 3?6 4?0 3? 7 Orbital Phase -90° -90° -90° Contour Plot Fig. 18a Fig. 18b Fig. 18ca Velocity Plots Fig. 19 Fig. 20 Fig. 21 ∆vproj 0.069 0.082 0.082 R 1.500 1.337 1.107 Projection Factor 257 206 249 Time until Merger 13 4.5 4 pro] aFigure 22 depicts Model 3 from (0o,0o) at an elapsed time 11.0 . 153 TABLE 7 ADOPTED PHYSICAL PARAMETERS K99 K564 Mean Redshift (∆λ∕λ) for Binary 0.01203 ± .00004 0.02619 ± .00005 Distance (Mpc h 1) 36.1 78.5 Projected Relative Velocity (km sec 1) 198 ± 13 26 ± 16 Projected Separation (kpc h 1) 10.3 ± 0.2a 13.3 ± 0.4a Mass Ratio (M^/M^) 0.2 iυ 1 Orbital Eccentricity <0.5 Ά, 0 Orbit Inclination 50° ± 5° 3?8 ± 0?2 -90° ± ? -100° ± 20° Orbital Phase Mass Projection Factor 2.1 ± 0.3b Total Mass ( 1011 M 1.9 ± 0.3 4.9 ± 0.6b 30 ? 10 ? 1.6 ± 0.3b 2.5 ± 0.4b Mτ∕Lτ ■L h^1) ( Mo h ∕L.) M1 ( 1011 M h^1) Virial Estimate for M. ( 10^1 M h-1) 1 ® Rotation Curve Estimate ( 10bb M h-1) Virial Estimate for M„ 2 Time until Merger ( 10 ( 1011 M g years ) ® h-1) 1.0 1.3 ? o.ιc 0.3 ± 0.06b M2 ( 1011 M0 h^1) 240 ± 30b 2.5 ± 0.4b 0.3 ? 1.1 4 4 Uncertainties in locating galaxy centers on CCD picture °Error quoted is that of the model-matching procedure. cNo correction was made for the tilt of the plane of rotation. 154 REFERENCES Borne, K. D. 1982, preprint (Paper I). Faber, S. M., and Jackson, R. E. 1976, Ap. J., 204, 668. Illingworth, G. 1981, in Structure and Evolution of Normal Galaxies, ed. S. M. Fall and D. Lynden-Bell (Cambridge University Press), p. 27. Karachentsev, I. D. 1972, Comm. Special Astrophys. Obs. U.S.S.R., 7, 3. Kormendy, J., and Illingworth, G. 1982, Ap, J., 256, 460. Sargent, W. L. W., Schechter, P. L., Boksenberg, A., and Shortridge, K. 1977, Ap. J., 212, 326. Schechter, P. L., and Gunn, J. E. 1979, Ap. J., 229, 472. Thuan, T. X., and Gunn, J. E. 1976, Pub. A. S. P., 88, 543. Toomre, A. 1977, in The Evolution of Galaxies and Stellar Populations, ed. B. M. Tinsley and R. B. Larson, (New Haven: Yale University Observatory), pp. 401-416. White, S. D. Μ., and Sharp, N. A. 1977, Nature, 269, 395. 155 FIGURE 1 Ο Ten-minute red (6500A) CCD photograph of K99 from the KPNO 2.1-m, taken with the Wide Field/Planetary Camera ground-based system. North is up. (reversed from the sky). East is to the right Photo measures 412 on a side. interest is the "tail" on the smaller galaxy. Of 156 157 FIGURE 2 Contour plot for K99 derived from a combination of the CCD picture depicted in Figure 1 and a two-minute CCD exposure. North and east are as marked. measures 30". The horizontal bar Contours represent red surface brightness and are separated by 0.75 mag. The five innermost levels are derived from the short exposure, the five outermost levels from the long exposure. Of interest are the opposing distentions in the light distributions for the two galaxies. 158 LU © <9 * Z <ô & & 159 FIGURE 3 K99 rotation curve, along the line connecting the two galaxies. Velocity units are km/sec, measured relative to cZq = 3600 km/sec. East is to the left. The center of the east component is located at R ≈ -27" and that of the west component at R ≈ 32". The west component reveals an unusually large rotation rate for a galaxy of its type (E2) Of interest also is the peculiar U-shaped rotation profile for the east component (see text for a full discussion). 160 K99 161 FIGURE 4 Spatial variations in the line-of-sight velocity dispersion for K99. Radial positions are the same as for the points shown in Figure 3. East is to the left. Velocity units are km/sec. Of interest is the positive gradient in the dispersion values for the east component. 162 K99 163 FIGURE 5 Ten-minute red CCD photograph of K564 from the KPNO 2.1-m. North is up and east is to the right (reversed from sky). Photo measures 4'.2 on a side. The small component to the southeast is probably not closely associated with the larger interacting pair (see text for discussion). 164 165 FIGURE 6 Contour plot for K564 derived from the CCD picture depicted in Figure 5. North and east are as marked. horizontal bar measures 30". The Contours represent red surface brightness and are separated by 0.75 mag. Of interest are the sweeping distentions on opposing sides of the two main galaxies. 166 167 FIGURE 7 K564 rotation curve, along the line connecting the two galaxies. Velocity units are km/sec, measured relative to czθ = 8000 km/sec. East is to the right. The center of the west component is located at R ≈ -18" and that of the east component at R ≈ 17". Of interest is the unusual V- shaped rotation profile for the west component. 168 K564 300 200 >- 100 -300 - _l___________ I__________ I__________ I___________ L -50 -25 0 25 50 R (fiRCSEC) 169 FIGURE 8 Spatial variations in the line-of-sight velocity dispersion for K564. Radial positions are the same as for the points shown in Figure 7. East is to the right. Velocity units are km/sec. 170 171 FIGURE 9 Projection of particle distribution for an initial binary configuration that has 1000 particles in the primary "galaxy" and 200 in the smaller component. measures 4 distance units on a side. is 1.265 units. The plot The binary separation 172 173 FIGURE 10 Contour plot for the same particle distribution as depicted in Figure 9. Plot measures 6 units on a side. Drawn contours are separated by 2-magnitude intervals. 174 175 FIGURE 11 Velocity profiles for an observer of the "galaxies" in Figures 9 and 10 from a view 90° ahead of the binary and 50° from the vertical. Velocities for the secondary component are on the left. Double solid lines are plus and minus one standard deviation from the mean line-of-sight velocities. from -2 to 2. The abscissa is projected radius and varies The ordinate for the rotation curve varies from -2 to 2, while that for the dispersion curve varies from 0 to 2. 176 177 FIGURE 12 Contour plots for Models 1 through 4 of K99. Details for these configurations are given in Table 4. The lowest contour level for Model 3 is 1.5 mag below that for the others. A few test particles appear isolated from the main bodies of the "galaxies". The long rectangle within the boundaries of each figure represents the slit used for the velocity measurements. Each box measures 4 units on a side and is centered on the binary center of mass. 178 a. Model 1 b. Model 2 c. Model 3 d. Model 4 ♦ 179 FIGURE 13 Velocity profiles for Model 1 of K99. Projection angles and other parameters for these plots appear in Table 4. 180 181 FIGURE 14 Velocity profiles for Model 2 of K99. Projection angles and other parameters for these plots appear in Table 4. 182 a. ROTATION CURVE 183 FIGURE 15 Velocity profiles for Model 3 of K99. Projection angles and other parameters for these plots appear in Table 4. 184 185 FIGURE 16 Velocity profiles for Model 4 of K99. Projection angles and other parameters for these plots appear in Table 4. 186 187 FIGURE 17 Top view for Model 3 of K99 at an elapsed time equal to 8.0. Detailed parameters for this simulation are to be found in Tables 3 and 4. 188 189 FIGURE 18 Same as Figure 12 except for the three models of K564. 6. Details for these configurations are given in Table The lowest contour in Model 3 is 1.5 mag above that for the others. 190 b. Model 2 V « c. Model 3 ♦ ♦♦ * 191 FIGURE 19 Velocity profiles for Model 1 of K564. Projection angles and other parameters for these plots appear in Table 6. 192 193 FIGURE 20 Velocity profiles for Model 2 of K564. Projection angles and other parameters for these plots appear in Table 6. 194 195 FIGURE 21 Velocity profiles for Model 3 of K564. Projection angles and other parameters for these plots appear in Table 6. 196 197 FIGURE 22 Top view for Model 3 of K564 at an elapsed time equal to 12.0. Detailed parameters for this simulation are to be found in Tables 5 and 6. 198