TURUN YLIOPISTON JULKAISUJA, ANNALES UNIVERSITATIS TURKUENSIS SARJA - SER. A I OSA - TOM. 409 ASTRONOMICA - CHEMICA - PHYSICA - MATHEMATICA SMALL AND LARGE SCALE DYNAMICS IN THE MILKY WAY by Esko Gardner TURUN YLIOPISTO UNIVERSITY OF TURKU Turku 2010 From Department of Physics and Astronomy University of Turku FI-20014 Turku Finland Supervised by Chris Flynn Docent Finnish Centre for Astronomy with ESO (FINCA) University of Turku Finland Reviewed by Karri Muinonen Alice Quillen Professor Associate Professor of Astronomy University of Helsinki Dept. of Physics and Astronomy Department of Physics University of Rochester Finland United States Opponent Hans Rickman Professor Space Research Center Polish Academy of Sciences Poland ISBN 978-951-29-4335-2 (PRINT) ISBN 978-951-29-4336-4 (PDF) ISSN 0082-7002 Multiprint - Turku, Finland 2010 Acknowledgements I would first like to thank my supervisor Chris Flynn for his encouraging style and support in the process of this thesis. I would also like to thank my other collaborators, mainly Pasi Nurmi and Seppo Mikkola for guidance and discussion from the very beginning and end of the process. I deeply thank the gracious hosting of Prof. Kimmo Innanen during my visit to York University in Toronto and Prof. Geraint Lewis during my stay at Sydney University. Programming was assisted by PathScale and Intel FORTRAN-compilers. Software was run on the University of Turku grid Topaasi. Additional anal- ysis courtesy of Gnumeric. I acknowledge the financial support of the Finnish Graduate School in Astronomy and Space Physics, the Finnish Cultural Foundation, and the Magnus Ehrnrooth Foundation. Musical entertainment provided by (in alphabetical order): Armin van Buuren, DJ Cotts, Matt Darey, and Sean Tyas, as well as the radios of DI.fm/trance and happyhardcore.com. I also wish to thank the peer-support groups of coffee table@Tuorla, deltan toimisto@ircnet, l5r@undernet, and hardcore@byteirc. Finally I would like to thank everyone who has in some way been a part of this process, you know who you are. Oh, and thanks mom. 3 In memoriam of Leena Ta¨htinen Contents Acknowledgements 3 Abstract 7 List of papers 9 1 Introduction 11 2 Methods 13 2.1 The Milky Way . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Basics of potentials . . . . . . . . . . . . . . . . . . . . . . . 14 2.3 A Milky Way potential . . . . . . . . . . . . . . . . . . . . . 14 2.4 Backwards restricted integration method . . . . . . . . . . . 19 2.5 Computational methods . . . . . . . . . . . . . . . . . . . . 20 3 Motion of the Sun 21 3.1 The Solar motion . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 The effect of the Solar motion on the outer Solar System . . 22 3.3 The Galactic tide’s effect on comets . . . . . . . . . . . . . 26 3.3.1 Cometary influx into the Solar system . . . . . . . . 27 3.3.2 Observing the cometary influx . . . . . . . . . . . . 29 4 Motion of the Sun in comparison to other nearby stars 31 5 The effect of a bar on the central region of the Galaxy 35 6 Bar-induced local effects 39 6.1 The bars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 6.2 The Hercules stream as a proxy for bar parameters . . . . . 42 6.2.1 The Galactic bar . . . . . . . . . . . . . . . . . . . . 43 6.2.2 The long bar . . . . . . . . . . . . . . . . . . . . . . 45 5 6.2.3 Both bars . . . . . . . . . . . . . . . . . . . . . . . . 46 6.3 The Hercules stream as derived from simulations . . . . . . 46 7 Discussion and future work 49 Bibliography 53 6 Abstract This thesis contains dynamical analysis on four different scales: the Solar system, the Sun itself, the Solar neighbourhood, and the central region of the Milky Way galaxy. All of these topics have been handled through methods of potential theory and statistics. The central topic of the thesis is the orbits of stars in the Milky Way. An introduction into the general structure of the Milky Way is presented, with an emphasis on the evolution of the observed value for the scale-length of the Milky Way disc and the observations of two separate bars in the Milky Way. The basics of potential theory are also presented, as well as a developed potential model for the Milky Way. An implementation of the backwards restricted integration method is shown, rounding off the basic principles used in the dynamical studies of this thesis. The thesis looks at the orbit of the Sun, and its impact on the Oort cloud comets (Paper IV), showing that there is a clear link between these two dynamical systems. The statistical atypicalness of the orbit of the Sun is questioned (Paper I), concluding that there is some statistical typicalness to the orbit of the Sun, although it is not very significant. This does depend slightly on whether one includes a bar, or not, as a bar has a clear effect on the dynamical features seen in the Solar neighbourhood (Paper III). This method can be used to find the possible properties of a bar. Finally, we look at the effect of a bar on a statistical system in the Milky Way, seeing that there are not only interesting effects depending on the mass and size of the bar, but also how bars can capture disc stars (Paper II). 7 8 List of papers This thesis consists of a review of the subject and the following original research articles: I J. A. Robles, C. H Lineweaver, D. Grether, C. Flynn, C. A. Egan, M. B. Pracy, J. Holmberg, and E. Gardner, A Comprehensive Comparison of the Sun to Other Stars: Searching for Self-Selection Effects, ApJ (2008), Erratum: ApJ (2008) arXiv:0805.2962. II E. Gardner, K. A. Innanen, and C. Flynn, Dynamical effects of the long bar in the Milky Way, in The Galaxy Disk in Cosmological Context, Proceedings of the International Astronomical Union, IAU Symposium, Volume 254. Poster Session. Cambridge: Cambridge University Press (2009) arXiv:0808.0498. III E. Gardner, and C. Flynn, Probing the Galaxy’s bars via the Hercules stream, MNRAS (2010), Erratum: Submitted arXiv:1002.0551. IV E. Gardner, P. Nurmi, C. Flynn, and S. Mikkola, Tidal influence of a variable galactic potential on the orbits of Oort cloud comets, JOURNAL (YEAR), arXiv code. 9 10 Chapter 1 Introduction Motion is an interesting subject: from early on in the physical development of human beings they have the need to move. In physics we call the study of movement ‘dynamics’. Since mankind realised that the phases of the Moon must be related to the movement of the Moon around the Earth, we have sought to quantify, identify, and learn what causes different kinds of motion. In astronomy, especially, we study gravitational motion. The study of dynamics can be of something that is very small, like the orbit of a Bohr model electron around the proton. Or something very large, with galaxies moving around each other in a cosmic dance. The dynamics featured here are more close to home, in our own Galaxy, known as the Milky Way (from Latin: via lactea). Staring at the sky, we can identify many different celestial bodies, from our own Sun, the other planets in the Solar System, comets and meteors. We can start to use telescopes and scientific analysis to find even larger systems, streams of stars moving as one, globular clusters, open clusters being born in large nebulæ, and giant molecular clouds. Finally, we can group whole systems of stars and gas together to form structure, such as discs, bars, bulges, spiral arms and halos. First, the basic methods of the research are approached: the Milky Way itself, the techniques, such as potential theory and numerical computation; as well as a model for the Galactic potential, and a statistical method to analyse the simulated motion of stars, known as the backwards restricted integration method. After the basics I will concentrate on the four different types of dynamics this thesis contains. From the smallest scale to the largest: Oort cloud comets, the Sun, the Solar vicinity, and lastly the Galactic disc. Each of these dynamical systems has been analysed under a different context. The first is studied via the influx of comets into the 11 12 Introduction Solar System, as generated by the movement of the Solar System around the Galaxy. The second is the motion of the Sun in comparison to other stars around us, studying statistically what our own movement is like in comparison to our nearby companions. The third is analysing how the bar can create dynamically associated streams in the stars of our local neighbourhood, as well as a study of the two different observed bars in the Milky Way. The fourth concentrates on the effect a bar would have on the Galactic disc, looking at how a distribution of stars would be changed by this effect. Lastly, a summary of the research as a whole will be presented with a discussion of the matter studied in this thesis, encompassing the whole range of dynamics and techniques studied here, as well as some future considerations. Chapter 2 Methods 2.1 The Milky Way The Milky Way is what we call our own home Galaxy, and it has been the subject of study in astronomy ever since the human race gazed at the night sky and began to wonder. During what we call the modern era, we have tried to study our own Galaxy in order to try to find what it looks like, just as we aim our telescopes at other ‘island universes’ in order to see what they look like. The general consensus of the structure of our Galaxy includes a disc, a bulge, a stellar halo, a dark matter halo, a bar, and spiral arms. There is a clear problem with trying to look at the structure of our own Galaxy. We cannot view our Galaxy from ‘outside’, so we cannot directly see the light and gas distribution for the whole Galaxy at once. Being quite close to the mid-plane of the Galactic disc also affects our ability to get a reliable image of large-scale structure. Work from the COBE/Diffuse Infrared Background Experiment (DIRBE) identified that our Galaxy has a bar in it (Binney et al., 1997; Bissantz & Gerhard, 2002), this feature was dubbed the ‘Galactic bar’. Later work by Benjamin et al. (2005) and Lo´pez-Corredoira et al. (2007) found yet an- other bar in the middle of our Galaxy. It did not match the morphology of the Galactic bar, and it was named the ‘long bar’, as it was calculated to be longer than the Galactic bar. It is difficult to say what exactly is going on in the centre, since, again, we cannot view the Milky Way from any other vantage point and because distance uncertainties and extinction make it difficult to construct an accurate map. The disambiguation is handled, although very lightly, in Paper III. 13 14 Methods 2.2 Basics of potentials The basic equations used in potential theory (see e.g. Binney & Tremaine 2008 for a wealth of information on the topic) are: F = −∇Φ (2.1) ∇2Φ = 4piGρ , (2.2) where F is force, Φ is the potential, R is the radius (e.g. distance from the Galactic centre), G is the gravitational constant and ρ is density. As one can understand from the equations, the potential represents some density distribution, that will generate, in our case, a gravitational force. Equation 2.2 is better known as Poisson’s equation. The circular velocity of the system, Vc is the rotational velocity a ‘tracer’ must have in order to maintain a circular orbit: V 2c = R|∇Φ| . (2.3) This equation is only approximate, as it is only valid for axisymmetric systems, problems arise from considering non-axisymmetric systems. Since equation 2.3 fails to accurately represent circular velocity, since there is no constant force at a certain radius (R). This requires more thought on what kind of orbit is circular, as true circular orbits do not really exist in a system where the density distribution is not symmetric. This was handled in Paper III. 2.3 A Milky Way potential The full equations of the axisymmetric potential are: Φ = ΦH +ΦC +ΦD +Φbar , (2.4) ΦH = 1 2 V 2h ln(r 2 + r20) , (2.5) Methods 15 ΦC = − GMC1√ r2 + r2 C1 − GMC2√ r2 + r2 C2 , and (2.6) ΦD = 3∑ i=1 −GMdi√ (R2 + (adi + √ (z2 + b2))2) , (2.7) where r = x2 + y2 + z2, and R = x2 + y2. ΦH ,ΦC ,ΦD, and Φbar are the respective potentials of the dark halo, the bulge, the disc, and the bar. The parameters for this system are in Table 2.1. The negative mass disk (ΦD2) is required to constrain the scale-length to lower values. Table 2.1: Parameters of the axisymmetric model. Parameter Value Unit Vh 220 km s −1 r0 8.5 kpc MC1 3 10 9 M rC1 2.7 kpc MC2 16 10 9 M rC2 0.42 kpc Md1 77.04 10 9 M ad1 5.81 kpc Md2 −68.48 10 9 M ad2 17.43 kpc Md3 26.75 10 9 M ad3 34.86 kpc b 0.3 kpc The first stage of this work was to study the new information we have on the generalised structure of the Galactic disc. It was quickly noticed that times have changed, and that what is called the scale-length (hR), i.e. 16 Methods the trend in the surface density (Σd) of the Galactic disc, has dramatically reduced from some initial works. For example, van der Kruit (1986) and Lewis & Freeman (1989) quote respective values of hR = 5.5 and 4.4 kilo- parsecs (one kiloparsec being about 3216 light years, or 3.1 × 1016 km). Most recently, the value of hR has settled around 3.3 kiloparsecs (kpc) from Lo´pez-Corredoira et al. (2002); Ruphy et al. (1996); Ojha (2001). But the Sloan Digital Sky Survey (SDSS) found an even shorter scale-length of 2.6 kpc (Juric´ et al., 2008). A more complete review of this information was gathered and presented in Paper II. A few assumptions had to be made in order to match the model to observations. This is because the model of the Milky Way disc should represent the actual Milky Way disc. The first one was the fixing of the solar position at (R0,z0) = (8,0) kpc, which is not consistent with the current estimates from e.g. Juric´ et al. (2008), although the values of R0 tend to vary greatly, with estimates ranging from 7 to 9 kpc (see e.g. Bica et al. 2006; Groenewegen & Blommaert 2005). The second was the fixing of goal values for local density (ρ0), and surface density (Σd). These were taken from Holmberg & Flynn (2000), where they state values of ρ0 = 0.102 ± 0.010 M /pc 3 and Σd = 48 M /pc 2. As the method of modifying the disc was mostly trial-and-error, I managed to achieve these values within a few per cent. Overall disc mass was the third assumption, where the current estimated value is around 4 ±1× 109 M (Flynn et al., 2006). This new information sparked modifications to the potential of the Milky Way in Flynn et al. (1996). Due to this being the starting point of this work, I attempted to reduce the disc scale-length from the original 4.4 kpc to something more in-line with modern estimates from e.g. Juric´ et al. (2008); Gould et al. (1997). The Miyamoto-Nagai potentials, which represent the disc (see Binney & Tremaine 2008 and equation 2.7), were modified in order to achieve this goal. The new model for the Galatic disc has a scale-length of around 3 kpc, as seen in Fig. 2.1. The rotation curve of the system is shown in Fig. 2.2. Some additional modifications to the potential were made for the work in Paper IV. This was because the vertical mass-profile of the potential did not match the observed vertical mass-profile from Holmberg & Flynn (2004) with sufficient accuracy for the particular problem adressed in Paper IV. This was achieved by adding an additional, thinner, disc-component to the potential. This represents the gas-layer in the disc of the Milky Way. Methods 17 The original potential was sufficiently accurate for the work in the other papers. Figure 2.1: The surface density of the disc component as a function of Galactocentric radius. The dashed line corresponds to an exponential den- sity falloff of 3 kpc, which is a good fit to the model over a wide range of radii. 3 kpc is also a fairly good representation of our current knowledge of the Galactic disc’s scale-length. Note that the density ends abruptly at 18 kpc. The case of the bar was more complex, but the work of Chandrasekhar (1969) describes the use of triaxial ellipsoid as useful representation for bars. The groundwork for an ellipsoidal potential was laid out by Ferrers (1877) and Dyson (1891) over 100 years ago. What is now known as a Ferrers’ potential is an ellipsoid with the following density law: ρ = ρ0(1−m 2)n m < 1, = 0 m ≥ 1, (2.8) 18 Methods 0 50 100 150 200 250 300 0 10 20 30 40 50 Vc (k m/ s) R (kpc) Dark Halo Spheroid Disk Total Figure 2.2: Circular velocities for the model, and the different contributions of the model components. m2 = x2 a2 + y2 b2 + z2 c2 a > b > c ≥ 0, (2.9) where a, b, and c are the semi-major axes of the ellipsoid. To get from this density based form to the force it generates, or from equation 2.2 to equation 2.1, is straightforward, but tedious. Analysis can be made very simple using constant density ellipses (n = 0), the potential equations for these can be found in both Chandrasekhar (1969) and Binney & Tremaine (2008). A constant density ellipse was used to emulate a bar in Paper II. It was later deemed that the more complex system presented in Pfenniger (1984), which is an n = 2 Ferrers’ potential, would more accurately repre- sent the bar. The main issue is that a constant density Ferrers’ potential has, at its edges, a rigid cutoff in density, which might create artifacts in the stellar motions. There were a few issues with the equations in the appendix of Pfenniger (1984), including a few typographical errors. Verification for Methods 19 the corrections, as well as answers to some other questions were obtained from Pfenniger (private communication). The equations of the potential are lengthy, and are given in Pfenniger (1984). The adopted parameters for the bars are presented later, in Table 6.1. 2.4 Backwards restricted integration method The backwards restricted integration method is the one used by Dehnen (2000) to investigate the local effects of the Galactic bar. The work by Dehnen (2000) is a major contributor to the inspiration for the work in Paper III. The basic idea of the method is to create an orbit library, using a wide range of initial conditions, and to apply a certain weight to each item in the library, calculating a total weight for each orbit. The combined orbits, weighted, lead to the density distribution and velocity distribution at various places in the Galaxy. An orbital library can be generated by integrating an orbit in a potential, recording the position and velocity of what is called a ‘tracer’. The tracer is a substitute for a star in the real Galactic potential. For each recorded point in position-velocity-space you may assign a weight: the product of the weights for the individual position- velocity-pairs, for example, is the weight of the orbit. Mathematically the backwards restricted integration method can be expressed as the following equation: Ω = N∑ i=1 ω(xi, vi), (2.10) where ω is the weight at a certain point along the orbit, Ω is the total weight of the orbit in the library, N is the amount of individual points in that orbit, and xi and vi are the respective position and velocity of the point. In Paper III I simulated ‘local velocity space’ (see Fig. 6.1 for the observed local velocity space). I used the potential described in the previous section to create an orbital library. I sampled the position and velocity of each orbit passing through ‘local velocity space’ for a billion years and chose weighting functions based on observed position and velocity distributions, i.e. on the distribution of the stars in phase space from observations of 20 Methods the Milky Way. This is done so that one can reconstruct the distribution of stars in velocity space at a certain point (the Solar neighbourhood), by essentially asking how many stars from different parts of the Galaxy have the right velocities to reach the solar neighbourhood. 2.5 Computational methods For the computation of the ellipsoid potential representing the bar I had to use a third-degree equation solver, which proved to be quite complex. The algorithm was translated from the Gnu Scientific Library (GSL) (Mark Galassi et al., 2006). The bar also requires the calculation of elliptical inte- grals. The technique by Carlson (1994) provides an accuracy of 10−4, which was deemed accurate enough for our calculations. Upon further analysis I found that the GSL uses the exact same numerical method as described by Carlson (1994) in its own algorithms. The orbital integration was per- formed using a fourth order Runge-Kutta algorithm. All of the simulation and statistical analysis software was written by the author in FORTRAN 90/95. Additional statistical analysis was performed in GNUMERIC. Chapter 3 Motion of the Sun The motion of the Sun around our own Galaxy has been an ongoing rid- dle. It was assumed that our Sun could be on a circular or near-circular orbit, e.g. Mayor (1974) measures a velocity (U, V,W ) for the Sun of (10.3, 6.3, 5.9) ± (1, 0.9, 0.4) km s−1, where U is the velocity towards the Galactic centre, V is the velocity along rotational direction, and W is the velocity perpendicular to the Galactic plane. The data obtained from the Hipparcos satellite (Dehnen & Binney, 1998), confirmed this assump- tion, measuring a velocity for the Sun of (10.00, 5.23, 7.17) ± (0.36, 0.62, 0.38) km s−1 . The motion described in Mayor (1974) and Dehnen & Bin- ney (1998) would put the Sun on a near-circular orbit. The later work of Hogg et al. (2005), working from the same Hipparcos data, confirmed this, measuring a velocity of the Sun of (10.1, 4.0, 6.7) ± (0.5, 0.8, 0.2) km s−1. Lately, though, it seems as the analysis method used in Dehnen & Bin- ney (1998) and Hogg et al. (2005) ’ might be flawed. The work of Scho¨nrich et al. (2010) clearly states that the original work, and method, in Dehnen & Binney (1998) should be revised. This is due to the assumption in Dehnen & Binney (1998) that the Galactic potential is axisymmetric, which, in practice, it is not (McMillan & Binney, 2010). Note that the authors of Scho¨nrich et al. (2010) include the authors of Dehnen & Binney (1998). So why would this matter, what difference does it make whether the Sun is moving on a circular orbit (eccentricity (e) ∼ 0), or an eccentric orbit (e.g. e > 0.1)? The interest comes not from a purely astronomical state, but rather a more recent emerging field called astrobiology. We recognise that our own existence, and life on earth itself is in some ways special. From the dynamical point of view, then, one could argue that the orbit of our Sun, for at least the past few hundred million years, must have kept us out of the more inhospitable regions of our Galaxy. The general statistical 21 22 Motion of the Sun study of our Sun as a host to a planet with life, was studied in Paper I. I was responsible for the dynamical calculations in Paper I. Paper I states that almost all (93 ± 1%) of the nearby stars have an eccentricity higher than the Sun’s, although the work in Paper I does not include the more recent information about the velocity of the Sun in the Galaxy (from e.g. Scho¨nrich et al. 2010; Binney 2010). The revised velocity had not been published yet. The statistics of orbits and the impact of the change in the Solar motion will be returned to in Chapter 4. 3.1 The Solar motion Before we move on to the effects of the Solar motion, we should look at the Solar motion itself. Figure 3.1 shows the orbit of the Sun in the purely axisymmetric (no bar included) potential of the Galaxy described in Section 2.3. The simulation ran for one billion years. The velocity used here, from Scho¨nrich et al. (2010), is (U, V,W )=(11.1, 12.24, 7.25) km s−1. The V -velocity is relative to what is called the Local Standard of Rest (LSR). This is the velocity that a star, at the radius from the Galactic centre, would require in order to stay on a circular orbit. The significant change with respect to Dehnen & Binney (1998); Hogg et al. (2005) is an increase in the Solar V -velocity of approximately 7 km s−1. This implies that the Sun is on a slightly more eccentric orbit, with simulations giving e = 0.059± 0.003 , compared to e = 0.036 ± 0.003 by us- ing the velocity from Dehnen & Binney (1998), and as used in Paper I. The new maximum and minimum Galactocentric radius are 8.9 ±0.5 and 7.9 ±0.5 kpc, respectively, compared to 8.4 ±0.2 and 7.9 ±0.2 kpc. The max- imum distance from the mid-plane is 0.101 ±0.04 kpc, compared to 0.103 ±0.03 kpc. One rotation period, around the Galaxy, is 223 million years, the radial and vertical periods are 153 and 94 million years, respectively. 3.2 The effect of the Solar motion on the outer Solar System What would happen if we were passing through more inhospitable regions of the Galaxy? What makes them inhospitable? How are we able to know that they could be detrimental to life itself on this planet? Motion of the Sun 23 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 7.8 8 8.2 8.4 8.6 8.8 9 z (kp c) R (kpc) Figure 3.1: Motion of the Sun in the axisymmetric model potential, in Galactocentric radius (R) and vertical height (Z). The initial velocity is from Scho¨nrich et al. (2010). The initial position of the Sun is (R,Z)0 = (8, 0) kpc. One of the more spectacular reoccurring events in the Solar System is comets. The Oort cloud, the parent population of long period comets, was proposed by Oort (1950). Work has been constantly made, e.g. by Fouchard et al. (2006); Rickman et al. (2008) to explain how exactly do the comets in the Oort cloud, which are orbiting at very large distances, evolve dynamically, through perturbations, into comets that we can observe in the Solar System. In Paper IV we looked at how the Solar orbit affects the evolution of long period Oort cloud comets. ‘Long period’ refers to comets that have large semi-major axes (a), of the order of tens of thousands of Astronomical Units (AU), one AU being the distance from the Earth to the Sun. The reason to study these comets is that the Galactic tide, the force of the surrounding 24 Motion of the Sun Galaxy on the volume around of the Sun (including the Oort cloud), is thought to be the one of the main factors for the dynamical evolution of these bodies. The force acting upon the comets is (from Heisler & Tremaine 1986): F = − GM r3 r −G1xxˆ−G2yyˆ −G3zzˆ (3.1) G1 = −(A−B)(3A+B) (3.2) G2 = (A−B) 2 (3.3) G3 = 4piGρ0 − 2(B 2 −A2), (3.4) where A and B are the Oort constants that depend on the underlying rotation curve of the Galaxy (see Binney & Tremaine 2008 for more details) and ρ0 is the local matter density. G is the gravitational constant. In the work presented in Paper IV, I used the improved Milky Way model and the knowledge of the Sun’s orbit and proceeded to see what the effect of the varying Galactic tide would have on the Oort cloud. The work requires a few assumptions, the first one is a potential model of the Milky Way, the second is the velocity of the Sun around the Galaxy. The first was presented in Section 2.3. The velocity I have taken from Scho¨nrich et al. (2010), as this seems to best reflect the current knowledge of the Solar motion. The model does not produce an identically flat rotation curve, as shown in Fig. 2.2. This is customary to adopt, so that G1 = G2 = 0. The Oort constants, A and B, will also vary over time due to the nature of the model’s rotation curve. As stated by Fouchard et al. (2006), the ρ0 term, and thus the G3 term will dominate the tidal force, as it is around ten times larger than the G1 and G2 parameters. Figure 3.2 shows the variation of the G-parameters over the calculated Solar orbit. The velocity used here, from Scho¨nrich et al. (2010), is (U, V,W )=(11.1, 12.24, 7.25) km s−1. The largest variation of the G3-parameter is caused by the radial motion, R (in a cylindrical (R, θ, Z)-system), where the minor variations are due to the vertical, Z, motion. Since the model in Paper IV does not include the bar, it is completely axisymmetric, and as such there is no dependence on azimuthal angle. Motion of the Sun 25 -7.8 x 10-16 -7.6 x 10-16 -7.4 x 10-16 -7.2 x 10-16 -7.0 x 10-16 -6.8 x 10-16 -6.6 x 10-16 -6.4 x 10-16 0 200 400 600 800 1000 G 1 (1/ yr2 ) Time (Myr) 6.4 x 10-16 6.6 x 10-16 6.8 x 10-16 7.0 x 10-16 7.2 x 10-16 7.4 x 10-16 7.6 x 10-16 7.8 x 10-16 8.0 x 10-16 8.2 x 10-16 0 20 40 60 80 100 G 2 (1/ yr2 ) Time (Myr) 3.0 x 10-15 3.5 x 10-15 4.0 x 10-15 4.5 x 10-15 5.0 x 10-15 5.5 x 10-15 6.0 x 10-15 6.5 x 10-15 0 20 40 60 80 100 G 3 (1/ yr2 ) Time (Myr) Figure 3.2: Variation of the G-parameters along the Solar orbit, using the Solar motion determined by Scho¨nrich et al. (2010). 26 Motion of the Sun 3.3 The Galactic tide’s effect on comets Once I determined how the Galactic tide changes over time, or just by choosing some constant value for the Galactic tide, I could input this force into a model of how the comets’ orbits evolve over a period of time. In Paper IV, the method described by Mikkola & Nurmi (2006) was used to integrate the motion of the comets. This method prefers to directly modify the orbital parameters of the comet, such as inclination and eccentricity, instead of classically integrating the orbit in position and velocity. The reasoning for this is simple. To gather the orbital parameters from position and velocity (or vice versa) one must solve Kepler’s equation, which re- quires one to iterate numerically. Although it is not a very time-consuming process to solve once, in simulations like this with 106 comets, with one simulation point per 105 years, it still equates to solving Kepler’s equa- tion 109 times. Thus, not solving Kepler’s equation, for every point in the simulation, makes the simulations run quicker. The study takes comets with certain semi-major axes (a), with the rest of the orbital parameters, such as inclination (i), longitude of the ascending node (Ω), argument of the periapsis (ω), and true anomaly (ν) are com- pletely randomly distributed. The values for a were chosen to represent the approximate location of the Oort cloud, 10000, 20000, 30000, 40000, 50000, and 60000 AU, corresponding to a maximum distance of 120000 AU, when e ∼ 1. At this point, it would be good to notice that whenever one talks about a ‘comet’ or a ‘star’ in the context of a simulation, it means that it is a certain point in position-velocity space that is affected by our simulation system, such as the potential, or, in this case, the force of the Galactic tide. The integration of each orbit was performed for a billion years, tracing the orbit of both the Sun and the comets over this time period. In another case, the effect of the Solar motion was ignored (i.e. the Sun was on a pure circular orbit, without vertical motion), and a mean value of the G-parameters was chosen to represent the Galactic tide. The rate of comets in the inner Solar System can be predicted by inte- grating the orbits of the long period comet population. The orbits of these comets can evolve so that they come close enough so that the planets start to affect the orbits of the comets. At this point the semi-major axis of the comet’s orbit might change so that it may be captured and become a ‘short-period’ comet, like Haley’s comet, or it might be transferred into a Motion of the Sun 27 hyperbolic orbit, never to be seen again. Our chosen limits were the per- ihelion distance (q), less than 30 AU (just inside Neptune’s orbit), and a distance from the Sun to be less than 1000 AU. Additionally, the Galactic tide is negligible inside 1000 AU, so the orbit does not evolve due to the tide, while the comet is within 1000 AU of the Sun. At a distance of 30 AU, the effect of the planets in the Solar System becomes non-negligible (Heisler & Tremaine, 1986; Fouchard et al., 2006). The last choice is so that q does not evolve to 30 AU, or less, and then evolves to a value that is higher than 30 AU, without ever coming near to the Sun. The latter limit is high due to the size of the integration steps, as following an exact motion, in time-scales of months, in simulations that try to handle millions of years is not feasible. The inner Solar System, including planets, is not included in the model, this is another reason why the q limit is set to 30 AU. 3.3.1 Cometary influx into the Solar system What can be first seen from the simulations, shown in Fig. 3.3, is that the comets take a while, approximately 200 million years, to achieve a state where they are more susceptible to drop into the inner Solar System. This is what can be described as a relaxation time. After the relaxation time it is possible to see the effects from the Galactic tide. As expected, when one keeps the values for the G-parameters constant, as seen in the bottom panel of Fig. 3.3, the in-fall of comets is regular, and does not have large-scale variations. The dynamic case, where the G-parameters vary with the motion of the Sun, on the other hand, follows fairly well the variation inG3. What this tells us is that the influx of comets, as a steady state event, does depend closely on the position of the Sun in the Galaxy. Other events, such as passing through high density regions, such as giant molecular clouds (Wickramasinghe & Napier, 2008) or spiral arms (Gies & Helsel, 2005) affect the influx as well, but they were not possible to test here, as our model does not explicitly include them. Other events include the close passage of another star, as described in Rickman et al. (2008). The only difference between the two simulation cases happens in the time-domain, showing that the basic case is valid. More comparisons between the two cases are shown in Paper IV. 28 Motion of the Sun 0 200 400 600 800 1000 Time (Myr) 0 100 200 300 400 500 600 700 800 F lu x 0 200 400 600 800 1000 Time (Myr) 0 100 200 300 400 500 600 700 F lu x Figure 3.3: The time evolution of the cometary flux into the inner So- lar System, the dynamic case (top panel), and the constant case (bottom panel). The solid line represents the flux of comets in 1 million year bins, the white line a 10 million year moving average. The dashed line repre- sents the evolution of the G3-parameter, one can see how the flux evolution corresponds to the evolution in G3. Motion of the Sun 29 3.3.2 Observing the cometary influx What is the general purpose of studying the cometary flux? And what does this information tell us now? To address the last point first, we can see two distinct periods in the cometary influx, one having a larger effect than the other. The radial period, of 149 ± 1 million years is the primary driver for long-term large-scale change in the basic influx-signal, while it is modulated by the shorter-term small changes from the 85 ± 4 million year vertical period. Since the disc is symmetric, it means that the actual signal is half of this period, i.e. 43 ± 2 million years. Now, could such a modulated signal in the cometary flux be found on Earth? For this we could take the cratering of the Earth or the fossil record Chang & Moon (2005); Napier (2006). The main problem with cratering is the low number statistics. Chang & Moon (2005) found a ∼26 million year period with 90 craters, while Napier (2006) found either a ∼24, ∼35, or ∼42 million year period with 40 craters. Trying to find such a good periodic record is not easy, if at all possible. Fossil records would be an indicator of mass-extinction, where a cometary impact, or similar high-energy event has destroyed large number of living species (Napier, 2006). Since the period of mass-extinctions, 30 ± 1 million years, approximately corresponding to the period in the cratering records, 33 ± 3 million years (Rampino & Stothers, 1984). As became apparent from Paper IV, it is possible to drive a long-term, orbital signal into the cometary flux. Therefore, if one wanted to arbitrarily choose a circular orbit, with only vertical movement, the signal will purely reflect that motion. The main issue seems to be that the proxy chosen to represent the change in the Galactic tide does not match the simulated Solar period. There must be some other process in question, ones that are not included in the model, such as passage through spiral arms, which might generate a suitable signal for the aforementioned cratering studies. Bailer- Jones (2009), though, in a recent review of this whole topic concludes that there is little evidence for any intrinsic periodicity in biodiversity, impact cratering, or climate on timescales of tens to hundreds of millions of years. I find that the Solar motion gives vertical periods in the range of 81−89 million years; the motion results in quasi-periodicity. Detecting such a periodicity in, e.g. the cratering record, would be difficult, and we agree with Bailer-Jones (2009) that the evidence of the vertical motion of the Sun in the cratering record is weak. 30 Motion of the Sun Chapter 4 Motion of the Sun in comparison to other nearby stars Paper I contains an analysis on how typical the Sun is. The tested char- acteristics included the mass, age, metallicity, two elemental abundance ratios, the rotational velocity, Galactocentric radius, eccentricity, and max- imum height above the Galactic plane of the Sun, as well as the mass of the host galaxy and the host galaxy group. All of these are compared to several thousand nearby stars. The overall conclusion is that the Sun is not particularly special, and that none of the tested characteristics is nec- essary for a star to host a planet with life. For the work in Paper I, I was responsible for the determination of the eccentricity (e) and the maximum height above the Galactic plane (Zmax) for the Sun and other stars in the Solar vicinity. There are a few improvements that can be made to the dynamical anal- ysis in Paper I. Work on the potential mentioned in section 2.3 had not been completed, or published, at the time when Paper I was being written. The ‘original’ potential from Flynn et al. (1996) was used then. Addition- ally, there was no bar in the model, relying on a pure axisymmetric system. The Solar velocity in Paper I was from Dehnen & Binney (1998), as well as from Hogg et al. (2005). Both of the velocities mentioned (Dehnen & Binney, 1998; Hogg et al., 2005), are not that different from each other, especially compared to Scho¨nrich et al. (2010). The statistics of the orbit of the Sun, and other stars, were reanalysed with the new potential for the disc. Furthermore, the simulation was run 31 32 Motion of the Sun in comparison to other nearby stars with and without the bar, where the bar case utilises the long bar, and not the Galactic bar. Additionally, the velocities, distances, and positions were used from the revised catalogue of Holmberg et al. (2009). The selection criteria are that a star is within 40 parsecs and it has a U , V , and W -velocity, this amounts to 1924 stars. The simulation was run for 2 billion years, sampling the position of the star every 10000 years. Eccentricity (e) is defined as (Rmax −Rmin)/(Rmax +Rmin), where Rmax is the maximum Galactocentric radius and Rmin is the minimum Galactocentric radius. Table 4.1: Statistical analysis of the Solar orbit in comparison to other nearby stars using both the original model from Flynn et al. (1996) and the new model mentioned in 2.3. Two different solar velocities, Dehnen & Binney (1998) and Scho¨nrich et al. (2010) were used. The percentile represents the amount of stars which have a higher value than the Solar value. Parameter Solar value Sample range Percentile Note e 0.061± 0.003 0.0008− 0.971 72.9± 2 % a,c,d e 0.058± 0.003 0.0021− 0.978 78.5± 2 % a,d e 0.036± 0.003 0.0014− 0.996 92.5± 1 % b,e e 0.058± 0.003 0.0020− 0.988 78.5± 2 % b,d Zmax 111.1 +6.1 −5.9 pc 3.7− 12778 pc 58.0± 2 % a,c,d Zmax 110.7 +6.2 −6.0 pc 3.7− 10431 pc 58.1± 2 % a,d Zmax 103.8 +5.9 −5.9 pc 2.5− 11285 pc 58.7± 3 % b,e Zmax 107.8 +6.0 −5.7 pc 3.6− 10292 pc 58.0± 2 % b,d a. Axisymmetric potential from section 2.3 b. Axisymmetric potential from Flynn et al. (1996) c. Barred potential added d. Solar velocity from Scho¨nrich et al. (2010) e. Solar velocity from Dehnen & Binney (1998) Looking at Table 4.1 it can noticed that there is only a slight difference between the barred and un-barred case. The conclusions of Paper I are Motion of the Sun in comparison to other nearby stars 33 unchanged whether or not one includes a bar in the Galaxy. These differ- ences are due to the perturbing effect of the bar. This effect will be fully covered in the following two Chapters. It is also interesting to note that the bar causes one of the stars to adopt an extremely circular orbit, with an eccentricity of only 0.0008. The greatest difference, despite the effect of the bar, in comparison to Paper I, is in the percentage of stars with e and Zmax higher than the Sun. The relative amount of stars on more eccentric orbits than the Sun has been reduced to 78.5 ± 2 %, where it was previously 92.5 ± 1% (first and second row, Table 4.1. Furthermore, the relative amount of stars on higher- reaching orbits than the Sun has decreased from 58.7 ± 3 % to 58.1 ± 2 % (fifth and sixth row, Table 4.1), this change is completely insignificant in comparison to the change in eccentricity. There are a few reasons for these changes. The main key is in the improvements made to the Geneva- Copenhagen survey from the first version (Nordstro¨m et al., 2004), to the second version (Holmberg et al., 2007), and finally to the third version used here (Holmberg et al., 2009). At the time of writing Paper I, only the first version was available. Holmberg et al. (2007) states that the velocities have been revised based on the recalibration of the metallicities and tempera- tures. Holmberg et al. (2009) additionally states that the velocities have been yet again re-calibrated to reflect the improved distance measurements. This also accounts for the difference in the number of simulated stars, from 1987 in Paper I to 1924 here. For these reasons, a reanalysis of the original simulations from Paper I is important. We do this by using the same initial conditions as mentioned above, but using the velocities from both Dehnen & Binney (1998) and Scho¨nrich et al. (2010) in the original disc model from Flynn et al. (1996). The simulations were reanalysed, with errors, in an unbarred case, and are also presented in Table 4.1 (the third, fourth, seventh, and eighth rows). Running a case where all initial Z-positions were offset by 20 pc had no effect on the statistical analysis. This tested what effect the adoption of a 20 pc Solar height has on the statistical outcome. There are minor differences between the statistics derived from the adopted models and Solar velocity. Changing the scale length of the disc has a very small effect on the solar value of Zmax, although the percentile stays within the same range. The greatest change comes from the differ- ence in the adopted Solar velocity. Eccentricity increases by over 60% (the first and second rows of Table 4.1), although it is still on a fairly circular 34 Motion of the Sun in comparison to other nearby stars orbit. The statistical change in eccentricity, from adopting a different So- lar velocity, is fairly significant, with the percentile value decreasing by 14 %. Overall, the effect of using the Solar motion determined by Scho¨nrich et al. (2010) makes the Sun more typical, as the eccentricity of the Solar orbit becomes less unusual than if one used the Solar motion determined by Dehnen & Binney (1998). Chapter 5 The effect of a bar on the central region of the Galaxy Benjamin et al. (2005) and Lo´pez-Corredoira et al. (2007) recently proposed that there is a second, ‘long bar’ in the Galactic centre. Benjamin et al. (2005) discovered the long bar in the Galactic Legacy Mid-Plane Survey Extraordinaire (GLIMPSE) survey, while Lo´pez-Corredoira et al. (2007) used observations from the Two Micron All Sky Survey (2MASS) to find the long bar. The long bar is longer, flatter, and much more dense than the traditional Galactic bar. The long bar lies at a different position angle, but weighs about the same as the Galactic bar. Because the long bar has only been recently detected, its dynamical effect on the disc of the Galaxy has not yet been studied. Here I discuss the effect of the long bar on disc stars passing through it. This also allowed me to develop a program for treating a rotating bar and stellar orbits near and far away from the bar. There is a small undocumented feature pertaining to bar models. Dehnen (2000); Minchev et al. (2007) both use a ‘quadrupole’ to approximate the bar. This is fine, as long as the stars do not get too close to the bar. When analysing both a quadrupole and a Ferrers’ potential, one can notice that the rotating Ferrers’ potential acts closely the same as the quadrupole, at large (R > 5 kpc) distances. Inside the bar, and closer in towards the tip of the bar, the effects change. Even a constant density Ferrers’ potential, used here, is a more accurate description of a bar-like density distribution than a quadrupole. Using the constant density (n = 0) Ferrers’ potential to emulate the bar, I simulated its impact on a distribution of stars. Each star was integrated for two billion years. Two cases were tested, one is a ‘thin’ case, equal 35 36 The effect of a bar on the central region of the Galaxy to the long bar, with axis sizes of (3.9:0.6:0.1) kpc. The other has been made ‘thicker’, with axis sizes of (3.9:2:1) kpc. Due to a constant density in the model, we know that the respective densities in the two cases are 6 M /pc 3 and 0.2 M /pc 3. The initial velocities represent a kinematically cold system, where the U , V , andW -velocities were selected from a normal distribution with a standard deviation of 10 km s−1. An 3 kpc exponential distribution in R represents the initial Galactic disc. A total of 10000 stars were simulated. Fig. 5.1 shows how the distributions of the radius (R) and the height (Z) evolve under the effect of a bar. The initial distribution stays almost the same in the model without a bar, with only a slight shift inward in Z. The radial profiles show more evolution. The thin case exhibits a clear shift inwards in the distribution from logR = 3.0 − 3.6. The thick case has a smaller shift at logR = 3.0. These correspond to radii of 1 kpc to 4 kpc. What this means is that the bars are capturing disc stars. A mere 10000 stars will not tell how they organise themselves in respect to the bar, but what this does tell us is that the general distributions depend on the morphology of a bar. A thinner bar will move stars inwards and upwards, while a thicker bar will only have a slight inward change and a systematic movement upwards in comparison to the case without a bar. In this case, there is far less evolution in the Z and R distributions, from their initial configuration. Still, the dynamical effects for the ‘thin’ long bar-emulating case are quite interesting, as the morphology of the inner disc could be changed by the long bar and that the long bar can affect stars very close to the Solar circle. Additionally, a bar with a density of 6 M /pc 3 is not so dense or so extreme, that it would have drastic effects on the centre of the Galaxy. The effect of a bar on the central region of the Galaxy 37 2 2.4 2.8 3.2 3.6 4 4.4 4.8 Log(R) 0 1000 2000 3000 C o u n t Initial configuration No Bar Thin Bar Thick Bar <0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Log(z) 0 500 1000 1500 2000 2500 C o u n t Initial configuration No Bar Thin Bar Thick Bar Figure 5.1: Initial and final distributions of radius (top panel) and height (bottom panel). The solid line represents the initial distribution, the dot- ted line represents the distribution without a bar, the dashed-dotted line represents the distribution with long bar dimensions, and the dashed line represents a thicker bar. The vertical line in the upper panel marks the position of the Sun. 38 The effect of a bar on the central region of the Galaxy Chapter 6 Bar-induced local effects Using my analytical model, I now explore the local effects of bar structure on the Galaxy. The statistical evidence in Chapter 4 shows that there is not much impact on the eccentricity or maximum height on 1924 stars within 40 pc of the Sun. However, such a small sample does not tell the full story on how stars are moving near the Sun. If I take a much larger sample, all the stars, with measured velocities, from the Geneva-Copenhagen sur- vey (Holmberg et al., 2009) and supplement that set with the stars from Schuster et al. (2006), one obtains a total of 14474 stars. Fig. 6.1 shows the resulting distribution of those stars if one arranges them in a two di- mensional grid, with a U range from −100 to 100 km s−1 and a V range from −150 to 50 km s−1 in bins of 2 km s−1. Note that this range only includes 13579 stars, as there are stars with very high velocities which do not fit into the figure. The contours represent percentages of the bin with the most stars in it, in steps of 5%. Fig. 6.1 shows over-densities at (−35,−5), (0,−10), (20,20) km s−1, and something extending out from the main area of contours centred at (−25,−30) km s−1. These are identified as the Hyades (Eggen, 1958), Pleiades (Eggen, 1972), Sirius (Eggen, 1958), and Hercules (Eggen, 1971b) streams. The Sun is marked in Fig. 6.1 as a white star. A good question to ask is where the streams come from, what causes them, and why? Paper III concentrates on the Hercules stream, a stream that has been suggested to be caused by a bar (Dehnen, 2000; Fux, 2001). If the connection to the bar is correct, one could use it to measure the mass, rotational velocity, and angle of the bar in the Milky Way. This would especially constrain the dynamical properties of such a bar, in particular the pattern speed and the bar’s rotational period compared to the Sun’s orbital period, as shown in Paper III. 39 40 Bar-induced local effects ★ Figure 6.1: Local velocity space, from the Geneva-Copenhagen survey Holmberg et al. (2009) and Schuster et al. (2006). The Hercules stream is the prominent feature area around (−25,−30) km s−1. The velocities have been calculated using the Solar velocity from Scho¨nrich et al. (2010). The Sun is marked with a white star. 6.1 The bars There is a clear problem with identifying the Hercules stream as being caused by the bar. As touched upon earlier, there are now thought to be two observed bars, the Galactic bar (Bissantz & Gerhard, 2002) and the long bar (Benjamin et al., 2005; Lo´pez-Corredoira et al., 2007). The long bar differs from the Galactic bar by being longer, highly flattened, and more dense. How two such massive bars could function dynamically, in the central region of the Galaxy, is uncertain. But first, we should consider looking at how the bars have been constrained observationally. Compiled in Table 6.1 are the observational constraints of the bars, the angle, dimensions, and mass, along with some of the other simulation parameters, which will be explained later. The dimensions of the Galactic Bar-induced local effects 41 bar are from Bissantz & Gerhard (2002), the angle is derived from Dehnen (2000), and the mass we adopt for the Galactic bar is constistent with Zhao (1996); Weiner & Sellwood (1999). The observational constraints of the long bar are exclusively from Lo´pez-Corredoira et al. (2007). There is a review on the topic of bar parameters in Vanhollebeke et al. (2009), especially table 1 in Vanhollebeke et al. (2009) provides various measured parameters for a bar at the Galactic centre. Table 6.1: Default simulation parameters for the bar models. Parameter Value Unit Galactic bar Position angle 25 ◦ Dimensions 3.5:1.4:1.0 kpc Mass 10 109 M Pattern speed 55.9 km s−1 kpc−1 Local circular velocity (Vc) 239 km s −1 OLR1 1.87 Long bar Position angle 43 ◦ Dimensions 3.9:0.6:0.1 kpc Mass 6 109 M Pattern speed 54.9 km s−1 kpc−1 Local circular velocity (Vc) 235 km s −1 OLR1 1.87 1 Outer Lindblad Resonance The aim here is to see if a Galactic bar, a long bar, or both bars can generate a Hercules feature. The use of a Ferrers’ potential as a bar was outlined in the previous Chapter and Paper II, here we modify the bar potential and adopt a more realistic (n = 2) model. 42 Bar-induced local effects Placing a Ferrers’ potential into the axisymmetric system is easy, but calculating the impact of it is not. The first change is the circular velocity: equation 2.3 states that the force will determine the circular velocity of the system. This circular velocity, in our part of the Galaxy, is called the Local Standard of Rest (LSR). If one puts a rotating non-axisymmetric component in a potential, then the force, at a certain radius, does not stay constant over time, and circular orbits do not exist. To deal with this, the local cirular velocity is defined by the V -velocity that produces the least eccentric orbit at a given radius. The local circular velocity (Vc) shown in Table 6.1 is exactly this, the velocity required to be on a near-circular orbit. The next parameter, the Outer Lindblad Radius (OLR), represents the interpretation of the rotation of a bar. The OLR itself is the ratio of the local rotational frequency to the rotational frequency of the bar. At an OLR value, for the bar, of 2, the star rotates around the Galaxy once, while the bar rotates twice around its axis in the plane of the Galaxy. The Hercules stream has been used before to measure the OLR of the Galactic bar, with Dehnen (2000) deriving a value of 1.85 ± 0.15. Minchev et al. (2007) used an independent method, the Oort C-constant, to derive a value of 1.87 ± 0.04. The pattern speed of the bar is the angular velocity of the bar. It can be easily calculated by using the position of the Sun, the LSR, the OLR, and the length of longest axis of the bar. 6.2 The Hercules stream as a proxy for bar pa- rameters Here we vary the parameters of the bar to see if we can model the features of the velocity distribution in the Solar neighbourhood. Using the backwards restricted method as described in Section 2.4, a similar image to Fig. 6.1 is generated. To create a model velocity distribution we must specify a weighting function, i.e. ω. The functions for ω are chosen based on the phase space density of stars in the Galaxy at any given position, R. The phase space density (ω) is the Bar-induced local effects 43 radial probability distribution of stars: ω = ωRωUωV ωR = e −(R−R0) hR , (6.1) ωU = e −U2 R 2σ2 U , (6.2) ωV = e −(VR−V (R)ass) 2 R 2σ2 V , (6.3) where R0 = 8 kpc, hR = 3 kpc, UR is the local U velocity, σU = 105 e −R 4.37kpc km s−1, VR is the local V velocity, σV = 74 e −R 3.36kpc km s−1, and V (R)ass is the assymetric drift at a given radius. The distributions are defined by gaussian fits to observed velocity dispersions and the asymmetric drift made by Lewis & Freeman (1989) over a wide range of radii. The simulations are in 2D, so Z is not taken into account. The asymmetric drift measured by Lewis & Freeman (1989) is a function of Galactic radius, and accounts for the fact that stars, on average, are moving slower than the circular velocity at that radius. A 2D velocity grid, with no vertical motion, with the same boundaries and bins as in Fig. 6.1, was chosen as the initial velocity condition for the simulations. All stars start at an initial Galactocentric radius of 8 kpc. There are no other assumptions made in initial position and velocity, as all of the assumptions are in the dynamical model, and the weights to interpret the orbits. It was assumed that the bar in the Galaxy has been stable for a billion years, and chose to study the orbits over the same billion year timescale. It is interesting to note, that the effect of the bar, over 500 million years, is much less pronounced (see e.g. Fux 2001). Increasing the timescale to two billion years makes the effects, as expected, more pronounced. These two cases can be seen in Paper III. 6.2.1 The Galactic bar The Galactic bar is more massive than the long bar as well as being more voluminous. The other difference is the initial angle. Paper III tries to fully see what the effects are, in all the simulated scenarios, so only a few basic cases are shown here. We first check our simulations by reproducing the the work of Dehnen (2000). Using the basic parameters for the Galactic bar from Table 6.1, we can 44 Bar-induced local effects Figure 6.2: The simulated distribution of local velocity space using the basic Galactic bar parameters from Table 6.1. A similar feature to the Hercules stream is seen to arise in local velocity space, this is due to the resonance with the bar. see two features, something one can call a ‘main’ feature, at around (U ,V ) = (0,0) km s−1, and a ‘secondary’, Hercules-resembling, feature at (U ,V ) = (−15,−35) km s−1. The first reaction is that the centre of the secondary feature is offset, in comparison to the observed feature in Fig. 6.1. Increase in mass is reflected in a shift of the secondary to more negative V -velocities and an increase of the value of the OLR reflects on the secondary in a similar way. An increase in angle moves the secondary to more negative values in U , and to more positive values of V . The effects are reversed when reducing the mass, OLR, or angle. What this means is that it would require an increase in the angle of the Galactic bar in order to move the secondary to its observed position in velocity space. The change in parameters is the same, independent of which bar one chooses to model. Paper III, again, reflects a more complete analysis of the topic. The overall picture seems to be in agreement with the findings of Dehnen (2000), as we can reproduce Bar-induced local effects 45 the Hercules feature and measure a similar OLR for the bar. 6.2.2 The long bar The characteristics of the long bar, in comparison to the classical Galactic bar are its longer, but more narrow, shape, slightly lower mass, and that its present position angle is not 25◦, but 43◦. The other main difference is that it is, observationally speaking, new, as first discovered by Benjamin et al. (2005). Figure 6.3: The simulated distribution of local velocity space using the basic long bar parameters from Table 6.1. The long bar can produce a Hercules stream-like feature, just as the Galactic bar can. The basic long bar case, in Fig. 6.3, produces a primary, and a secondary feature. This time, the centre of the feature is approximately at (U ,V ) = (−25,−30) km s−1 showing an approximate agreement with the observed position of the Hercules stream. There is also another feature, albeit a very weak one. This tertiary feature, at V ∼ −90 km s−1 turns out to be caused by a bar at a high position angle. The only stream remotely in 46 Bar-induced local effects that area is the Arcturus stream (Eggen, 1971a). Williams et al. (2009) state that the Arcturus stream is at V ∼ −100 km s−1, with a wide range of U -velocities. Arifyanto & Fuchs (2006) and Helmi et al. (2006) found independently another stream, resembling the Arcturus stream, this stream is also located at V ∼ −100 km s−1. These streams could be associated with the long bar. 6.2.3 Both bars Simulating two bars at once is simple, although observationally speaking fairly odd. Cases where a galaxy has two co-existent, large scale bars seem to be non-existant. The rarity is based on searches through article databases, such as arXiv and ADS, as well as discussions with researchers in the field (Salo & Rautiainen, Private communication). At a more basic assumption, one would take the observed properties of both bars to be ac- curate. Using this approach, by adding both bars together, results in too much of an effect on local stars. Fig. 6.4 shows that there is a definite effect from using both bars, showing a large over-density at positive values of U . In this case, the bars rotate at the same frequency, or are ‘phase locked’. Paper III handles some unlocked cases, when there are two different rotation rates for the bars. In this case, additional structure is seen in the model velocity distribution that is not een in the real one. A suitable remedy would be to halve the masses of the bars. Fig. 6.5 shows that, in that specific case, velocity space starts to to look more realistic. It still has some minor problems, the Hercules feature is at slightly high V -velocities, but this can be remedied by increasing the OLR, as it moves the feature out, without adding mass. The tertiary feature becomes more weak. Since the long bar mass has been halved, its effects become less pronounced. 6.3 The Hercules stream as derived from simula- tions The cause of the generation of the Hercules feature, in both cases, is a resonance with the bar. The resonance is actually the lack of stars in a certain area, seen in the simulations as a void between the primary and Bar-induced local effects 47 Figure 6.4: The simulated distribution using both bar parameters from Table 6.1. Together, both bars produces extremely strong structures in the Solar neighbourhood velocity distribution. This is inconsistent with observations, and rules out a case where both bars are at least 1 billion years old. the secondary. The resonance means that there is a regular ‘kick’ caused by the bar that will move the stars to different orbits. They do not have a low enough V -velocity to go in and interact directly with the bar. The tertiary feature, as it turned out in further analysis, is caused by a direct interaction with the long bar (or the Galactic bar, if you shift its angle to match that of the long bar). These are stars in the simulation that reach to Galactic radii that are less than the length of the long bar, causing direct interactions. Dehnen (2000) used a quadrupole to model the bar, this is a poor approximation at small radii, and so while appropriate for modelling nearly circular orbits, highly eccentric orbits may not have been well integrated in his bar potential. For example, the Arcturus stream (Williams et al., 2009) stars have an eccentricity of ∼0.5 and pass near to the bar. A more accurate bar model, such as considered here, is required 48 Bar-induced local effects Figure 6.5: The simulated distribution using both bar parameters from Table 6.1, with the mass of both bars halved. This simulation is in much better agreement with observations of local stars. to accurately study the dynamics of higher eccentricity stars. Overall, the long bar is able to reproduce the Hercules feature as a secondary component, and also seems to add an explanation for a stream at V ∼ 100 km s−1. The Galactic bar seems to not quite accurately represent the Hercules stream. The basic Galactic bar case, exhibits a too high U - velocity for the Hercules feature. Increasing the position angle would shift the feature towards the observed location of the Hercules stream. The battle of the bars is even more open if you accept that both could co-exist. Chapter 7 Discussion and future work This thesis has aimed to sample some of the small and large scale dynamics happening in the Milky Way. Different dynamical systems are presented, unified by both their relation to each other, and through the basic technique of test particle integration. Large and small scale dynamics are equally represented: the impact of one dynamical motion on another, the impact of structure on various dynamical areas, and the statistical analysis of orbits of stars. The statistical analysis of the Sun in comparison to other nearby stars shows that the Solar orbit is a little unusual, but not significantly so. The data from Paper I, and from Section 4, show that although there might be some argument towards the typicalness, or otherwise, of the Solar orbit, there are still issues to be dealt with. The change in disc scale-length has little to no impact, but the change in the Solar motion, and the inclusion of the bar slightly changes the Solar orbit. The impact of the Solar motion on the Solar System has a strong effect on the influx of comets into the inner Solar System. In Paper IV we found that the influx of new comets due to the Galactic tide is dependent on both the radial and the vertical Solar motion. It is clear that both are important. Rickman et al. (2008) argue that perturbations by passing stars should be a more primary engine for the influx of comets into the inner Solar System. So is it even possible to see the past motion of the Sun through the Galaxy via a suitable proxy? Additionally, what role would a passage through a spiral arm, or a giant molecular cloud have? These are questions for the future. The bar in our Galaxy has a clear effect on the structure of the Galactic disc and it would be interesting to observe this. What kind of structure should we be able to see, elsewhere in the disc, due to the dense long bar? 49 50 Discussion and future work There is a clear effect based on bar geometry and bar density. It thus might be possible to constrain the mass of the bar using future surveys, such as Gaia, The Herschel Multi-tiered Extragalactic Survey (HerMES), and Large Synoptic Survey Telescope (LSST). There are computational problems to address in this type of simulations, as Paper II included only 10000 stars, a larger statistical sample would give a more complete picture on the effects the long bar has on the disc, as well as how the inner region evolves due to the bar. The effect of the bar on the Solar neighbourhood is observed via the study of dynamical streams. Paper III simulates the effects of the Galactic and long bar, separately and together, on local velocity space. As the observed Hercules Stream is a dynamical feature in local velocity space. The variation in mass, angle, and rotation rate of the bar have an impact of the simulated Hercules Feature. The simulations were performed in-plane, although the framework is set up for more complete 3D-simulations. Again, it would be computationally more expensive to perform more complete 3D- simulations, but that would require us to know better the velocities of stars at different heights and locations in our Galaxy. There are also interesting implications towards low V -velocity streams, as these can get close enough to the bar to interact directly, such streams could be the ones observed to be at V ∼ −100 km s−1 by Eggen (1971a) or Arifyanto & Fuchs (2006); Helmi et al. (2006). Large upcoming star surveys, such as Gaia and LSST could bring more light to these dynamical streams. There is still more work to be carried out to complete the model of the Milky Way. The primary missing features are spiral arms. The issues related to this are twofold. The primary is a technical challenge, to figure out what represents spiral arms accurately enough to use in modelling. The second one is observational, there seem to be a considerable diversity of opinions even on the amount of arms in our galaxy, how they move around the Galaxy, how dense the separate arms are, where the arms are, especially concerning the ones on the other side of the Galaxy. There can also be discussions on how to relate to N -body simulations. N -body offers a more realistic picture, by trying to emulate each star. The handling of the respective positions, velocities, and forces of 1011, or more, objects should be soon a reality. Computational power has increased, but even simple large-scale-universe calculations tend to take forever. Modern General Purpose Graphics Processing Units (GP-GPUs) are approaching an interesting limit, and harnessing them for calculations, for potential theory, Discussion and future work 51 or otherwise, is a completely new area of computational astrophysics. A back-of-the-envelope calculation shows that one could run ∼ 106 orbits in the barred potential, within a day, by using one GP-GPU. This would provide easy access to large number statistics, for such applications as used in Paper II. GP-GPUs could even be harnessed to simulate large catalogues at once. With advances in both observational and computational methods, there is clearly much more work to be done refining our knowledge in the structure and dynamics of the Milky Way. 52 Discussion and future work Bibliography Arifyanto M. I., Fuchs B., 2006, A&A, 449, 533 Bailer-Jones C. A. L., 2009, International Journal of Astrobiology, 8, 213 Benjamin R. A., Churchwell E., et al. 2005, ApJL, 630, L149 Bica E., Bonatto C., Barbuy B., Ortolani S., 2006, A&A, 450, 105 Binney J., 2010, MNRAS, 401, 2318 Binney J., Gerhard O., Spergel D., 1997, MNRAS, 288, 365 Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Galactic Dynamics: Second Edition, by James Binney and Scott Tremaine. ISBN 978-0-691-13026-2 (HB). Published by Princeton Uni- versity Press, Princeton, NJ USA, 2008. Bissantz N., Gerhard O., 2002, MNRAS, 330, 591 Carlson B. C., 1994, ArXiv Mathematics e-prints Chandrasekhar S., 1969, Ellipsoidal figures of equilibrium. The Silliman Foundation Lectures, New Haven: Yale University Press, 1969 Chang H., Moon H., 2005, PASJ, 57, 487 Dehnen W., 2000, AJ, 119, 800 Dehnen W., Binney J. J., 1998, MNRAS, 298, 387 Dyson F. D., 1891, Q. J. Pure Appl. Math., 25, 259 Eggen O. J., 1958, MNRAS, 118, 65 Eggen O. J., 1971a, PASP, 83, 271 53 54 BIBLIOGRAPHY Eggen O. J., 1971b, PASP, 83, 251 Eggen O. J., 1972, ApJ, 173, 63 Ferrers N. M., 1877, Q. J. Pure Appl. Math., 14, 1 Flynn C., Holmberg J., Portinari L., Fuchs B., Jahreiß H., 2006, MNRAS, 372, 1149 Flynn C., Sommer-Larsen J., Christensen P. R., 1996, MNRAS, 281, 1027 Fouchard M., Froeschle´ C., Valsecchi G., Rickman H., 2006, Celestial Me- chanics and Dynamical Astronomy, 95, 299 Fux R., 2001, A&A, 373, 511 Gies D. R., Helsel J. W., 2005, ApJ, 626, 844 Gould A., Bahcall J. N., Flynn C., 1997, ApJ, 482, 913 Groenewegen M. A. T., Blommaert J. A. D. L., 2005, A&A, 443, 143 Heisler J., Tremaine S., 1986, Icarus, 65, 13 Helmi A., Navarro J. F., Nordstro¨m B., Holmberg J., Abadi M. G., Stein- metz M., 2006, MNRAS, 365, 1309 Hogg D. W., Blanton M. R., Roweis S. T., Johnston K. V., 2005, ApJ, 629, 268 Holmberg J., Flynn C., 2000, MNRAS, 313, 209 Holmberg J., Flynn C., 2004, MNRAS, 352, 440 Holmberg J., Nordstro¨m B., Andersen J., 2007, A&A, 475, 519 Holmberg J., Nordstro¨m B., Andersen J., 2009, A&A, 501, 941 Juric´ M., Ivezic´ Zˇ., Brooks A., Lupton R. H., Schlegel D., et al. 2008, ApJ, 673, 864 Lewis J. R., Freeman K. C., 1989, AJ, 97, 139 Lo´pez-Corredoira M., Cabrera-Lavers A., Garzo´n F., Hammersley P. L., 2002, A&A, 394, 883 BIBLIOGRAPHY 55 Lo´pez-Corredoira M., Cabrera-Lavers A., Mahoney T. J., Hammersley P. L., Garzo´n F., Gonza´lez-Ferna´ndez C., 2007, AJ, 133, 154 Mark Galassi et al. 2006, GNU Scientific Library Reference Manual (2nd Ed.). Network Theory Limited Mayor M., 1974, A&A, 32, 321 McMillan P. J., Binney J. J., 2010, MNRAS, 402, 934 Mikkola S., Nurmi P., 2006, MNRAS, 371, 421 Minchev I., Nordhaus J., Quillen A. C., 2007, ApJL, 664, L31 Napier W. M., 2006, MNRAS, 366, 977 Nordstro¨m B., Mayor M., Andersen J., Holmberg J., Pont F., Jørgensen B. R., Olsen E. H., Udry S., Mowlavi N., 2004, A&A, 418, 989 Ojha D. K., 2001, MNRAS, 322, 426 Oort J. H., 1950, Bull. Astron. Inst. Netherlands, 11, 91 Pfenniger D., 1984, A&A, 134, 373 Rampino M. R., Stothers R. B., 1984, Nature, 308, 709 Rickman H., Fouchard M., Froeschle´ C., Valsecchi G. B., 2008, Celestial Mechanics and Dynamical Astronomy, 102, 111 Ruphy S., Robin A. C., Epchtein N., Copet E., Bertin E., Fouque P., Guglielmo F., 1996, A&A, 313, L21 Scho¨nrich R., Binney J., Dehnen W., 2010, MNRAS, 403, 1829 Schuster W. J., Moitinho A., Ma´rquez A., Parrao L., Covarrubias E., 2006, A&A, 445, 939 van der Kruit P. C., 1986, A&A, 157, 230 Vanhollebeke E., Groenewegen M. A. T., Girardi L., 2009, A&A, 498, 95 Weiner B. J., Sellwood J. A., 1999, ApJ, 524, 112 Wickramasinghe J. T., Napier W. M., 2008, MNRAS, 387, 153 56 BIBLIOGRAPHY Williams M. E. K., Freeman K. C., Helmi A., the RAVE collaboration 2009, in J. Andersen, J. Bland-Hawthorn, & B. Nordstro¨m ed., IAU Symposium Vol. 254 of IAU Symposium, The Arcturus Moving Group: Its Place in the Galaxy. pp 139–144 Zhao H., 1996, MNRAS, 278, 488 Original papers a rX iv :0 8 0 5 .2 9 6 2 v 2 [a s tr o -p h ] 2 2 M a y 2 0 0 8 A comprehensive comparison of the Sun to other stars: searching for self-selection effects Jose´ A. Robles1, Charles H. Lineweaver1, Daniel Grether2, Chris Flynn3, Chas A. Egan2,4, Michael B. Pracy4, Johan Holmberg5 and Esko Gardner3 ABSTRACT If the origin of life and the evolution of observers on a planet is favoured by atypical properties of a planet’s host star, we would expect our Sun to be atypical with respect to such properties. The Sun has been described by previous studies as both typical and atypical. In an effort to reduce this ambiguity and quantify how typical the Sun is, we identify eleven maximally-independent properties that have plausible correlations with habitability, and that have been observed by, or can be derived from, sufficiently large, currently available and representative stellar surveys. By comparing solar values for the eleven properties, to the resultant stellar distributions, we make the most comprehensive comparison of the Sun to other stars. The two most atypical properties of the Sun are its mass and orbit. The Sun is more massive than 95 ± 2% of nearby stars and its orbit around the Galaxy is less eccentric than 93± 1% of FGK stars within 40 parsecs. Despite these apparently atypical properties, a χ2-analysis of the Sun’s values for eleven properties, taken together, yields a solar χ2 = 8.39±0.96. If a star is chosen at random, the probability that it will have a lower value (∼ be more typical) than the Sun, with respect to the eleven properties analysed here, is only 29± 11%. These values quantify, and are consistent with, the idea that the Sun is a typical star. If we have sampled all reasonable properties associated with habitability, our result suggests that there are no special requirements for a star to host a planet with life. Subject headings: Sun: fundamental parameters — Sun: general — stars: fundamental param- eters — stars: statistics 1Planetary Science Institute, Research School of Astronomy & Astrophysics and Research School of Earth Sciences, The Australian National University, Canberra Australia; josan@mso.anu.edu.au. 2University of New South Wales, Sydney, Aus- tralia. 3Tuorla Observatory, University of Turku, Finland. 4Research School of Astronomy & Astrophysics, The Australian National University, Canberra, Aus- tralia. 5Max Planck Institute for Astronomy, Heidelberg, Germany. 1. INTRODUCTION If the properties of the Sun are consistent with the idea that the Sun was randomly se- lected from all stars, this would indicate that life needs nothing special from its host star and would support the idea that life may be common in the universe. More particularly, if there is nothing special about the Sun, we have little reason to limit our life-hunting ef- forts to planets orbiting Sun-like stars. As 1 I an example of the type of anthropic reason- ing we are using, consider the following sit- uation. Suppose uranium (a low abundance element in the Solar System and in the uni- verse) was central to the biochemistry of life on Earth. Further, suppose that a compar- ison of our Sun to other stars showed that the Sun had more uranium than any other star. How should we interpret this fact? The most reasonable way to proceed would be to try to evaluate the probability that such a co- incidence happened by chance and to deter- mine whether we are justified in reading some importance into it. Although a correlation does not necessarily imply cause, we think that a correlation between the Sun’s anoma- lous feature and life’s fundamental chemistry would be giving us important clues about the conditions necessary for life. Specifically, the search for life around other stars as envisioned by the NASA’s Terrestrial Planet Finder or ESA’s Darwin Project and as currently un- derway with SETI, would change strategy to focus on the most uranium-rich stars. An- other example: suppose the Sun had the high- est [Fe/H] of all the stars that had ever been observed. Then high [Fe/H] would be strongly implicated as a precondition for our existence, possibly by playing a crucial role in terrestrial planet formation. These are exaggerated ex- amples of the more subtle correlations that a detailed and comprehensive comparison of the Sun with other stars could reveal. Whether the Sun is a typical or atypi- cal star with respect to one or a few prop- erties has been addressed in previous stud- ies. Using an approach similar to ours (com- paring solar to stellar properties from partic- ular samples), some studies have suggested that the Sun is a typical star (Gustafsson 1998; Allende Prieto 2006), while other stud- ies have suggested that the Sun is an atypical star (Gonzalez 1999a,b; Gonzalez et al. 2001). This apparent disagreement arises from three problems: i) the language used to describe whether the Sun is, or is not typical, is often con- fusingly qualitative. For example, re- porting the Sun as “metal-rich”, can mean that the Sun is significantly more metal-rich than other stars (e.g. more metal-rich than 80% of other stars) or it can mean that the Sun is insignifi- cantly metal-rich (e.g. more metal-rich than 51% of other stars). ii) selection effects: the stellar samples cho- sen for the comparison can be biased with respect to the property of interest. iii) the inclusion (or exclusion) of stellar properties for which it is suspected or known that the Sun is atypical, will make the Sun appear more atypical (or typical). In this paper we address problem i by us- ing only quantitative measures when compar- ing the Sun’s properties to other stars. Our main interest is to move beyond the qualita- tive assessment of the Sun as either typical or atypical, and obtain a more precise quantifi- cation of the degree of the Sun’s (a)typicality. In other words, we want to answer the ques- tion ‘How typical is the Sun?’ rather than ‘Is the Sun typical or not?’ There are at least two ways to quantify how typical the Sun is. This can be done for individual parameters by determining how many stars have values be- low or above the solar value (Table 3). This can also be done by a joint analysis of multi- ple parameters (Table 2). If there are several subtle factors that have some influence over habitability, a quantitative joint analysis of the Sun’s properties may allow us to identify these factors without invoking largely specu- lative arguments linking specific properties to habitability. With respect to problem ii, most previous analyses have compared the Sun to subsets of Sun-like stars selected to be Sun-like with 2 respect to one or more parameters. In such analyses, the Sun will appear typical with respect to any parameter(s) correlated with one of the pre-selected Sun-like parameters. For example, elemental abundances [X/H] are correlated with metallicity1 [Fe/H]. The sam- ple of Edvardsson et al. (1993a) was selected to have a wide range of [Fe/H]. This pro- duced a metallicity distribution unrepresen- tative of stars in general. Recognizing this, Edvardsson et al. (1993a) conditioned on so- lar metallicity, [Fe/H] ≈ 0 and then compared solar abundances for 12 elements to the abun- dances in a group of nearby stars with solar iron abundance, solar age and solar galacto- centric radius. They found the Sun to be “a quite typical star for its metallicity, age and galactic orbit”. Similarly, Gustafsson (1998), after comparing various properties of the Sun to solar-type stars (stars of similar mass and age), concluded that the Sun seems very nor- mal for its mass and age; “The Sun, to a re- markable degree, is solar type”. The stellar samples we use for comparison with the Sun are, in our judgement, the least-biased sam- ples currently available for such a comparison. To address problem iii, in Section 2 we compare the Sun to other stars using a large number (eleven) of maximally-independent properties with plausible correlations with habitability. These properties can be observed or derived for a sufficiently large, representa- tive stellar sample (Table 1). Any property of the Sun or its environment which must be special to allow habitability would show up in our analysis. However, in contrast to previous analyses which have looked for solar anomalies with respect to individual proper- ties, we perform a joint analysis that enables us to quantify how typical the solar values are, taken as a group. In Section 3, the differ- 1Metallicity: [Fe/H] is the fractional abundance of Fe relative to hydrogen, compared to the same ratio in the Sun: [Fe/H] ≡ log(Fe/H) − log(Fe/H) ences between the solar values and the stellar samples’ medians are used to perform first a simple and then an improved version of a χ2- analysis to estimate whether the solar values are characteristic of a star selected at ran- dom from the stellar samples. The results of our joint analysis are presented in Figure 13 of Section 4. We find that the solar values, taken as a group, are consistent with the Sun being a random star. However, there are important caveats to this interpretation associated with the compromise between the number of prop- erties analyzed, and their plausibility of being correlated with habitability. In Sections 5 and 6 we discuss these caveats and summarize. We discuss the levels of correlation between our eleven properties in Appendix A. 2. Stellar Samples and Solar Values We are looking for a signal associated with a prerequisite for, or a property that favors, the origin and evolution of life (see Gustafsson 1998 for a brief discussion of this idea). If we indiscriminately include many properties with little or no plausible correlation with habit- ability, we run the risk of diluting any poten- tial signal. If we choose only a few properties based on previous knowledge that the Sun is anomalous with respect to those properties, we are making a useful quantification but we are unable to address problem iii. We choose a middle ground and try to identify as many properties as we can that have some plausible association with habitability. This strategy is most sensitive if a few unknown stellar proper- ties (among the ones being tested) contribute to the habitability of a terrestrial planet in orbit around a star. An optimal quantitative comparison of the Sun to other stars would require an unbiased, large representative stellar sample from which independent distributions, for as many prop- erties as desired, could be compared. Such 3 T a b le 1 : S a m p le s u se d to p ro d u c e th e st e ll a r d is tr ib u ti o n s p lo tt e d in F ig u re s 1 – 1 0 . F ig . P ro p e rt y R a n g e M e d ia n σ 6 8 † S o la r # st a rs S p e c tr a l d m a x S o u r c e µ 1 / 2 V a lu e ty p e [p c ] 1 M a ss [M  ] 0 .0 8 – 2 0 .3 3 0 .3 7 1 1 2 5 A 1 – M 7 7 .1 H e n ry 2 0 0 6 (R E C O N S ) 2 A g e [G y r] 0 – 1 5 5 .4 3 .2 5 4 .9 + 3 .1 − 2 .7 a 5 5 2 F 8 – K 2 2 0 0 R o c h a -P in to e t a l. 2 0 0 0 b 3 [F e / H ] − 1 .2 0 – + 0 .4 6 -0 .0 8 0 .2 0 0 4 5 3 F 7 – K 3 2 5 G re th e r & L in e w e a v e r 2 0 0 7 4 A [C / O ] − 0 .2 2 – + 0 .3 2 0 .0 7 0 .0 9 0 2 5 6 F G 1 5 0 b G 9 9 , R 0 3 , B F 0 6 4 B [M g / S i] − 0 .1 8 – + 0 .1 4 0 .0 1 0 .0 4 0 2 3 1 F G 1 5 0 c R 0 3 , B 0 5 5 v si n i [k m s− 1 ] 0 – 3 6 2 .5 1 1 .2 7 1 .2 8 d 2 7 6 F 8 – K 2 8 0 e V a le n ti & F is c h e r 2 0 0 5 6 e 0 – 1 0 .1 0 0 .0 5 0 .0 3 6 ± 0 .0 0 2 f 1 ,9 8 7 A 5 – K 2 4 0 g N o rd s tr o¨ m e t a l. 2 0 0 4 7 Z m a x [k p c ] 0 – 9 .6 0 0 .1 4 0 .1 0 0 .1 0 4 ± 0 .0 0 6 h 1 ,9 8 7 A 5 – K 2 4 0 g N o rd s tr o¨ m e t a l. 2 0 0 4 8 R G a l [k p c ] 0 – 3 0 4 .9 5 .0 3 7 .6 2 ± 0 .3 2 i — — 5 0 ,0 0 0 j B S 8 0 , G 9 6 , E 0 5 9 M g a l [M  ]k 1 0 7 – 1 0 1 2 1 0 1 0 .2 0 .4 7 1 0 1 0 .5 5 ± 0 .1 6 — — 1 0 7 l D 9 4 , C B 9 9 , L 0 0 , B J 0 1 , J 0 3 1 0 M g r o u p [M  ]k 1 0 9 – 1 0 1 3 1 0 1 1 .1 0 .4 7 1 0 1 0 .9 1 ± 0 .0 7 — — 1 0 7 E k e e t a l. 2 0 0 4 † C h a ra c te ri st ic w id th o f d is tr ib u ti o n in th e d ir e c ti o n o f th e so la r v a lu e . a W ri g h t e t a l. (2 0 0 4 ), (s e e fo o tn o te in S e c . 2 .2 ). b G 9 9 : G u st a fs so n e t a l. 1 9 9 9 , R 0 3 : R e d d y e t a l. 2 0 0 3 , B F 0 6 : B e n sb y & F e lt z in g 2 0 0 6 . c R 0 3 : R e d d y e t a l. 2 0 0 3 , B 0 5 : B e n sb y e t a l. 2 0 0 5 . d S o la r ro ta ti o n a l v e lo c it y c o rr e c te d fo r ra n d o m in c li n a ti o n (s e e S e c . 2 .5 ). e S u b -s e t o f st a r s w it h in th e m a ss ra n g e : 0 .9 M  ≤ M ≤ 1 .1 M  . f C a lc u la te d u si n g th e so la r g a la c ti c m o ti o n (D e h n e n & B in n e y 1 9 9 8 ) a n d th e G a la c ti c p o te n ti a l (s e e S e c . 2 .6 ). g S u b -s e t o f v o lu m e c o m p le te A 5 – K 2 st a r s w it h in 4 0 p c . h In te g ra te d so la r o rb it in th e G a la c ti c p o te n ti a l o f F ly n n e t a l. (1 9 9 6 ) (s e e S e c . 2 .6 ). i E is e n h a u e r e t a l. 2 0 0 5 . j B S 8 0 : B a h c a ll & S o n e ir a 1 9 8 0 , G 9 6 : G o u ld e t a l. 1 9 9 6 , E 0 5 : E is e n h a u e r e t a l. 2 0 0 5 . k S te ll a r m a ss , n o t to ta l b a ry o n ic m a ss , n o r to ta l m a ss . l D 9 4 : D ri v e r e t a l. 1 9 9 4 , C B 9 9 : C o u rt e a u & v a n d e n B e rg h 1 9 9 9 , L 0 0 : L o v e d a y 2 0 0 0 , B J 0 1 : B e ll & d e J o n g 2 0 0 1 , J 0 3 : J a rr e tt e t a l. 2 0 0 3 . a distribution for each property of interest would allow a straightforward analysis and outcome: the Sun is within the n% of stars around the centroid of the N -dimensional dis- tribution. However, observational and sample selection effects prevent the assembly of such an ideal stellar sample. In this study, we compare the Sun to other stars with respect to the following eleven ba- sic physical properties: (1) mass, (2) age, (3) metallicity [Fe/H], (4) carbon-to-oxygen ratio [C/O], (5) magnesium-to-silicon ratio [Mg/Si], (6) rotational velocity v sin i, (7) ec- centricity of the star’s galactic orbit e, (8) maximum height to which the star rises above the galactic plane Zmax, (9) mean galacto- centric radius RGal, (10) the mass of the star’s host galaxy Mgal, (11) the mass of the star’s host group of galaxies Mgroup. These eleven properties span a wide range of stel- lar and galactic factors that may be associ- ated with habitability. We briefly discuss how each parameter might have a plausible cor- relation with habitability. For each property we have tried to assemble a large, represen- tative sample of stars whose selection crite- ria is minimally biased with respect to that property. For each property the percentage of stars with values lower and higher than the solar value are computed. For proper- ties (9), (10) and (11), the uncertainties in the percentages are determined from the un- certainties of the distributions. For the rest of the properties, nominal uncertainties ∆, on the percentages were calculated assum- ing a binomial distribution (e.g. Meyer 1975): ∆ = (nlow × nhigh/Ntot) 1/2 where nlow (nhigh) is the fraction of stars with a lower (higher) value than the Sun and Ntot is the total num- ber of stars in the sample. The solar value is indicated with the symbol “” in all figures. We compare the Sun and its environment to other stars and their environments. The analysis of these larger environmental con- texts provides information about properties that otherwise could not be directly measured. For example, suppose the metallicity of the Sun were normal with respect to stars in the solar neighborhood but that these stars as a group, had an anomalously high metallicity with respect to the average metallicity of stars in the Universe. This fact would strongly sug- gest that habitability is associated with high metallicity, but our comparison with only lo- cal stars would not pick this up. In the ab- sence of an [Fe/H] distribution for all stars in the Universe, we use galactic mass as a con- venient proxy for any such property that cor- relates with galaxy mass. 2.1. Mass Mass is probably the single most impor- tant characteristic of a star. For a main se- quence star, mass determines luminosity, ef- fective temperature, main sequence life-time and the dimensions, UV insolation and tem- poral stability of the circumstellar habitable zone (Kasting et al. 1993). Low mass stars are intrinsically dim. Thus a complete sample of stars can only be ob- tained out to a distance of ∼ 7 parsecs (≈ 23 lightyears). Figure 1 compares the mass of the Sun to the stellar mass distribution of the 125 nearest main sequence stars within 7.1 pc, as compiled by the RECONS consortium (Henry 2006). Over-plotted is the stellar Initial Mass Function (IMF) (Kroupa 2002, Eqs. 4 & 5, Table 1) normalised to 125 stars more mas- sive than the brown dwarf limit of 0.08 M. Since the IMF appears to be fairly univer- sal (Kroupa & Weidner 2005), these nearby comparison stars are representative of a much larger sample of stars. There is good agree- ment between the histogram and the IMF — the Sun is more massive than 95 ± 2% of the nearest stars, and more massive than 94±2% of the stars in the Kroupa (2002) IMF. 5 Fig. 1.— Mass histogram of the 125 near- est stars (Henry 2006, RECONS). The me- dian (µ1/2 = 0.33 M) of the distribution is indicated by the vertical grey line. The 68% and 95% bands around the median are indi- cated respectively by the vertical dark grey and light grey bands. We also use these con- ventions in Figs. 2–11. The solid curve and hashed area around it represents the Initial Mass Function (IMF) and its associated un- certainty (Kroupa 2002). The Sun, indicated by “”, is more massive than 95±2% of these stars. Fourteen brown dwarfs and nine white dwarfs within 7.1 pc were not included in this sample. Including them yields 94% — the same result obtained from the IMF. Our 95%± 2% result should be compared with the 91% reported by Gonzalez (1999b). The Sun’s mass is the most anomalous of the properties studied here. 2.2. Age If the evolution of observers like ourselves takes on average many billions of years, we might expect the Sun to be anomalously old (Carter 1983). Accurate estimation of stel- lar ages is difficult. For large stellar surveys (> a few hundred stars), the most commonly used age indicators are based on isochrone fit- ting and/or chromospheric activity (R′HK in- dex). Rocha-Pinto et al. (2000b) have esti- mated a Star Formation Rate (SFR), or equiv- alently, an age distribution for the local Galac- tic disk from chromospheric ages of 552 late- type (F8–K2) dwarf stars in the mass range 0.8 M ≤ M ≤ 1.4 M at distances d ≤ 200 pc (Rocha-Pinto et al. 2000a). They applied scale-height corrections, stellar evolution cor- rections and volume incompleteness correc- tions that converted the observed age distri- bution into the total number of stars born at any given time. Hernandez et al. (2000) and Bertelli & Nasi (2001) have made esti- mates of the star formation rate in the solar neighborhood and favour a smoother distri- bution (fewer bursts) than Rocha-Pinto et al. (2000b). In Figure 2 we compare the chromo- spheric age of the Sun (τ = 4.9 ± 3.0 Gyr, Wright et al. 2004) 2 to the stellar age distribution representing the Galactic SFR (Rocha-Pinto et al. 2000b). The median of this distribution is 5.4 Gyr. The Sun is younger than 53 ± 2% of the stars in the thin disk of our Galaxy. Over-plotted is the cosmic SFR derived by Hopkins & Beacom (2006). According to this distribution with a median µ1/2 = 9.15 Gyrs, the Sun was born after 86± 5% of the stars that have ever been born. The Galactic and cosmic SFRs are differ- ent because the cosmic SFR was dominated by bulges and elliptical galaxies in which the largest fraction of stellar mass in the Universe resides. Bulges and elliptical galaxies (early- 2To ensure that the Sun’s age is determined in the same way as the stellar ages to which it is being compared, we adopt the chromospheric solar age τ = 4.9 ± 3.0 Gyr over the more accurate meteoritic age τ = 4.57± 0.002 Gyr (Alle`gre et al. 1995). 6 type galaxies) formed their stars early and quickly and then ran out of gas. The disks of spiral galaxies, like our Milky Way, seem to have undergone irregular bursts of star for- mation over a longer period of time as they interacted with their satellite galaxies. The volume limited (dmax = 40 pc) sub- set from Nordstro¨m et al. (2004) contains isochrone ages for 1126 A5–K2 stars. The median of this sub-set is 5.9 Gyr and the Sun is younger than 55 ± 2% of the stars. The similarity of this isochrone age result to the chromospheric age result is not obvious since the agreement between these two age tech- niques is rather poor. This mismatch can be seen in Fig. 15D, Reid et al. (2007), and in Fig. 8 of Feltzing et al. (2001). Fig. 2.— The Galactic stellar age dis- tribution (median µ1/2 = 5.4 Gyr) from Rocha-Pinto et al. (2000b). The Sun is younger than 53± 2% of the stars in the disk of our Galaxy. The grey curve is the cosmic Star Formation Rate (SFR) with its associ- ated uncertainty (Hopkins & Beacom 2006), according to which the Sun is younger than 86± 5% of the stars in the Universe. 2.3. Metallicity Iron is the most frequently measured el- ement in nearby stars. Metallicity [Fe/H], is known to be a proxy for the fraction of a star’s mass that is not hydrogen or he- lium. In the Sun and possibly in the Uni- verse, the dominant contributors to this mass fraction in order of abundance are: O(44%), C(18%), Fe(10%), Ne(8%), Si(6%), Mg(5%), N(5%), S(3%) (Asplund et al. 2005; Truran & Heger 2005). The corresponding abundances by number are: O(48%), C(26%), Ne(7%), N(6%), Mg (4%), Si(4%), Fe(3%), S(2%). Importantly for this analysis, this short list contains the dominant elements in the composition of terrestrial planets (O, Fe, Si and Mg) and life (C, O, N and S). Over the last few decades, much effort has gone into determining abundances in nearby stars for a wide range of elements. Stellar ele- mental abundances for element X are usually normalised to the solar abundance of the same element using a logarithmic abundance scale: [X/H] ≡ log(X/H) − log(X/H). Hence all solar elemental abundances [X/H], are de- fined as zero. Spectroscopic abundance anal- yses are usually made differential relative to the Sun by analysing the solar spectrum (re- flected by the Moon, asteroids or the telescope dome) in the same way as the spectrum of other stars. In this approach, biases intro- duced by the assumption of local thermody- namic equilibrium (LTE), largely cancel out for Sun-like stars (Edvardsson et al. 1993b). A comparison between solar and stellar iron abundances is a common feature of most abundance surveys and most have con- cluded that the Sun is metal-rich compared to other stars (Gustafsson 1998; Gonzalez 1999a,b). However, for our purposes, the appropriateness of these comparisons de- pends on the selection criteria of the stel- lar sample to which the Sun has been com- pared. Stellar metallicity analyses such 7 as Edvardsson et al. (1993a); Reddy et al. (2003); Nordstro¨m et al. (2004); Valenti & Fischer (2005) have stellar samples selected with dif- ferent purposes in mind, e.g., Edvardsson et al. (1993a) aimed to constrain the chemical evo- lution of the Galaxy and their sample is bi- ased towards low metallicity (average [Fe/H]= −0.25). The sample of Valenti & Fischer 2005 (average [Fe/H]= −0.01), was selected as a planet candidate list and contains some bias towards high metallicity (see Grether & Lineweaver 2007). To assess how typical the Sun is, Gustafsson (1998) limited the sample of Edvardsson et al. (1993a) to stars with galac- tocentric radii within 0.5 kpc of the solar galactocentric radius, and to ages between 4 and 6 Gyrs. The distribution of stars given by this criteria has an average [Fe/H]= −0.09. Grether & Lineweaver (2006, 2007) com- piled a sample of 453 Sun-like stars within 25 pc. These stars were selected from the Hip- parcos catalogue, which is essentially com- plete to 25 pc for stars within the spec- tral type range F7–K3 and absolute mag- nitude of MV = 8.5 (Reid 2002). Metallic- ities for this sample were assembled from a wide range of spectroscopic and photometric surveys. In Figure 3, we compare the Sun to the Grether & Lineweaver (2007) sample, which has a median [Fe/H]= −0.08. To our knowledge this is the most complete and least- biased stellar spectroscopic metallicity dis- tribution. The Sun is more metal-rich than 65± 2% of these stars. This result should be compared with Favata et al. (1997) who constructed a volume- limited (dmax = 25 pc) sample of 91 G and K dwarfs ranging in color index (B−V ) between 0.5− 0.8 (Favata et al. 1996). Their distribu- tion has a median [Fe/H]= −0.05 and com- pared to this sample, the Sun is more metal rich than 56 ± 5% of the stars. Fuhrmann (2008) compared the Sun to a volume com- plete (dmax = 25 pc) sample of about 185 thin-disk mid-F-type to early K-type stars down to MV = 6.0. He finds a mean [Fe/H] = −0.02± 0.18. This mean [Fe/H] is lowered by 0.01 dex if the 43 double-lined spectroscopic binaries in his sample are included. His re- sults are consistent with ours. Fig. 3.— Stellar metallicity histogram of the 453 FGK Hipparcos stars within 25 pc (Grether & Lineweaver 2007). The median µ1/2 = −0.08. The Sun is more metal-rich than 65± 2% of the stars. 2.4. Elemental ratios [C/O] and [Mg/Si] The elemental abundance ratios of a host star have a major impact on its proto- planetary disk chemistry and the chemical compositions of its planets. Oxygen and carbon make up ∼ 62% of the Solar Sys- tem’s non-hydrogen-non-helium mass content (Z = 0.0122, Asplund et al. 2005). Car- bon and oxygen abundances are among the hardest to determine. This is due to high temperature sensitivity and non-LTE effects in their permitted lines (e.g. C I λ6588, O I λ7773), and to the presence of blends in the forbidden lines ([C I] λ8727, [O I] 8 λ6300). See Allende Prieto et al. (2001) and Bensby & Feltzing (2006) for details on C and O abundance derivations. Carbon pairs up with oxygen to form car- bon monoxide. In stars with a C/O ratio larger than one, most of the oxygen con- denses into CO which is largely driven out of the incipient circumstellar habitable zone by the stellar wind. In this oxygen-depleted scenario, planets formed within the snow- line are formed in reducing environments and are mostly composed of carbon compounds, e.g. silicon carbide (Kuchner & Seager 2005). Thus, the C/O ratio could be strongly associ- ated with habitability. As most heavy element abundances rela- tive to hydrogen (e.g. [O/H], [C/H], [N/H]) are correlated with [Fe/H], they were not in- cluded in our analysis. After the overall level of metallicity (represented by [Fe/H]), and af- ter the ratio of the two most abundant metals, [C/O], the magnesium to silicon ratio [Mg/Si] is the most important ratio of the next most abundant elements (excluding the noble gas Ne). For example [Mg/Si] sets the ratio of olivine to pyroxene which determines the abil- ity of a silicate mantle to retain water (Hugh O’Neill, private communication). Stellar elemental abundance ratios are defined as [X1/X2] = [X1/H] − [X2/H]. Hence, systematic errors associated with the determination of absolute solar abundances cancel for abundances relative-to-solar. We compile [C/O] and [Mg/Si] ratios from sam- ples with the largest number of stars and high- est signal-to-noise stellar spectra: [C/O]: 256 stars from Gustafsson et al. 1999; Reddy et al. 2003; Bensby & Feltzing 2006 [Mg/Si]: 231 stars from Reddy et al. 2003; Bensby et al. 2005 Due to their selection criteria, these sam- ples are biased towards low metallicity and therefore cannot be used to create a represen- tative [Fe/H] distribution. Because a correla- tion exists between the [C/O] and [Mg/Si] ra- tios and [Fe/H] (e.g. Gustafsson et al. 1999), the samples we use have a relatively nar- row range of [Fe/H] to reduce the influence of the correlation. Therefore, these small correlations can be neglected in this study — see the bottom panels of Fig. 4 where [Fe/H] versus [C/O] as well as [Fe/H] ver- sus [Mg/Si] are plotted. The top panels show the corresponding stellar distribution histograms. The Sun’s [C/O] ratio is lower than 81 ± 3% of the stars. This is consis- tent with Gonzalez (1999b) who suggested — based on data from Edvardsson et al. (1993a) and Gustafsson et al. (1999) — that the Sun has a low [C/O] ratio relative to Sun-like stars at similar galactocentric radii. See however, Ramı´rez et al. (2007) who find that the Sun is oxygen poor compared to solar metallicity stars. The Sun’s [Mg/Si] ratio is lower than 66± 3% of the stars. The [C/O] and [Mg/Si] ratios are also largely independent of each other (see Fig. 14 in Appendix A). 2.5. Rotational velocity Stellar rotational velocities are related to the specific angular momentum of a proto- planetary disk and possibly to the magnetic field strength of the star during planet forma- tion, and to protoplanetary disk turbulence and mixing. An unusually low stellar rota- tional velocity may be associated with the presence of planets (Soderblom 1983). One or several of these factors could be related to habitability. There is a known correlation between mass and v sin i at higher stellar masses (e.g. see Fig. 18.21 of Gray 2005, p. 485). In order to minimise the effect of this correlation (and maximize independence between parameters), 9 Fig. 4.— A: Comparison of the Sun’s carbon-to-oxygen ratio ([C/O] ≡ 0) to the [C/O] ratios of 256 stars compiled from Gustafsson et al. (1999); Reddy et al. (2003) and Bensby & Feltzing (2006). The Sun’s [C/O] ratio is lower than 81±3% of the stars in this sample which has a median µ1/2 = 0.07. B: Comparison of the Sun’s magnesium-to-silicon ratio ([Mg/Si] ≡ 0), to [Mg/Si] values from 231 stars from Reddy et al. (2003) and Bensby et al. (2005). The Sun’s [Mg/Si] ratio is lower than 66 ± 3% of the stars in this sample with median µ1/2 = 0.01. The bottom panels C & D show the small correlations of these distributions with [Fe/H]. These small correlations can be neglected for this study. 10 we assembled a sample containing 276 stars within the mass range 0.9–1.1 M (F8–K2) from Valenti & Fischer (2005). The selec- tion criteria of the Valenti & Fischer (2005) stars introduces some bias against more ac- tive stars. We compared the high v sin i tail of our Valenti & Fischer (2005) sample with the high v sin i tail of a sub-sample from Nordstro¨m et al. (2004). We estimate that for our Valenti & Fischer (2005) sample, the bias introduced by the selection criteria is lower than ∼ 5%. The v sin i values in Valenti & Fischer (2005) are obtained by fix- ing the macroturbulence for the stars of a given color, without modeling the stars indi- vidually. If the macroturbulence value was underestimated for T > 5800 K, the result- ing v sin i values (especially when v sin i is near zero) would be overestimated (Sec. 4 of Valenti & Fischer 2005). The inclination of the stellar rotational axis to the line of sight is usually unknown so the observable is v sin i. Using the so- lar spectrum reflected by the asteroid Vesta, Valenti & Fischer (2005) derived a solar v sin i = 1.63 km s−1. For the purposes of this analysis we use the mean value that would be derived for the Sun, when viewed from a random in- clination: v sin i = 1.63(π/4) km s −1 ≈ 1.28 km s−1. The Sun rotates more slowly than 83± 7% of the stars in our Valenti & Fischer (2005) sample (Fig. 5). This is in agreement with Soderblom (1983, 1985) who reported that the Sun is within one standard deviation of stars of its mass and age. 2.6. Galactic orbital parameters The Galactic velocity components of a star (U ,V ,W ) with respect to the local standard of rest (LSR) may be used to compute a star’s orbit in the Galaxy. How typical or atypical is the solar orbit compared to the orbits of other nearby stars in the Galaxy? The orbit may be Fig. 5.— Rotational velocity histogram for 276 F8–K2 (0.9 ≤ M ≤ 1.1 M) stars (Valenti & Fischer 2005). The Sun (v sin i = 1.28 km s−1) rotates more slowly than 83±7% of the stars. There is one star to the right of the plot with v sin i = 36 km s−1. related to habitability because more eccentric orbits bring a star closer to the Galactic cen- ter where there is a larger danger to life from supernovae explosions, cosmic gamma and X- ray radiation and any factors associated with higher stellar densities (Gonzalez et al. 2001; Lineweaver et al. 2004). For a standard model of the Galactic po- tential, Nordstro¨m et al. (2004) computed or- bital paramters for the Sun, and for a large sample (∼ 16700) of A5–K2 stars. Their adopted components of the solar velocity relative to the local standard of rest were (U, V,W ) = (10.0±0.4, 5.25±0.62, 7.17±0.38) km s−1 (Dehnen & Binney 1998). For each of the 1,987 stars within 40 pc in the Nordstro¨m et al. (2004) catalog, an in- ner and outer radii Rmin and Rmax were com- puted. This yielded the orbital eccentricity e ≡ (Rmax −Rmin)/(Rmin + Rmax). The solar 11 eccentricity was computed using the compo- nents of the solar motion (Dehnen & Binney 1998) relative to the local standard of rest in the Galactic potential of Flynn et al. (1996). The bottom panel of Figure 6 shows the cor- relation between Galactic orbital eccentricity e and the magnitude of the galactic orbital ve- locities with respect to the local standard of rest: vLSR ≡ (U 2+V 2+W 2)1/2. Eccentricity e and vLSR are strongly correlated. We include e, not vLSR, in the analysis since e is less cor- related with the maximum height above the Galactic plane Zmax, than is vLSR. This is shown in Fig. 16 in Appendix A. The Sun’s eccentricity was determined with the same relation as the stellar eccentrici- ties. The uncertainty in our estimate of solar eccentricity came from propagating the un- certainty in the adopted solar motion. We find e = 0.036 ± 0.002 (consistent with the e = 0.043 ± 0.016 found by Metzger et al. 1998). The Sun has a more circular orbit than 93±1% of the A5–K2 stars within 40 pc (with median eccentricity µ1/2 = 0.1). This is the second most anomalous of the eleven solar properties we consider here. The frequency of the passage of a star through the thin disk could be associated with Galactic gravitational tidal perturbations of Oort cloud objects that might increase the impact rate on potentially habitable planets. This is correlated with the maximum height, Zmax, to which the stars rise above the Galac- tic plane. Figure 7 shows the stellar distri- bution of Zmax for the stars shown in Fig- ure 6. We find that 59 ± 3% of the A5–K2 stars within 40 pc of the Sun reach higher above the Galactic plane than the Sun does (Zmax, = 0.104 ± 0.006 kpc). The solar Zmax, was derived by integrating the solar or- bit in the Galactic potential. The uncertainty on W , produces the uncertainty on Zmax, and hence the ±3% uncertainty on 59%. Our results for eccentricity and Zmax are consistent Fig. 6.— Top panel: eccentricity distribu- tion for the 1,987 stars at d ≤ 40 pc from Nordstro¨m et al. (2004). The Sun has a more circular orbit than 93±1% of the A5–K2 stars within 40 pc. After mass, eccentricity is the second most anomalous parameter. Bottom panel: Correlation between vLSR and eccen- tricity for the same stars presented in the top panel. Since these properties are highly corre- lated we select only one for the analysis. The large grey point with error bars represents the median and the 68% widths of the two one- dimensional distributions. As in Fig. 4 the contours correspond to 38%, 68%, 82% and 95%. 12 Fig. 7.— The distribution of maximum heights above the Galactic plane for the Nordstro¨m et al. (2004) sample. 59 ± 3% of nearby A5–K2 stars (dmax =40 pc) reach higher above the Galactic plane than the Sun reaches. There are 22 stars evenly distributed over Zmax between 1.5 and 9.6 kpc. Their ex- clusion from the comparison reduces the 59% result by less than 1%. with those obtained using Hogg et al. (2005) LSR values: (U, V,W ) = (10.1 ± 0.5, 4.0 ± 0.8, 6.7 ± 0.2). Using the Hogg et al. LSR values, 92 ± 1% of A5–K2 stars within 40 pc have higher eccentricities than the Sun and 62 ± 4% of A5–K2 stars within 40 pc have larger Zmax values. How does the Sun’s distance from the cen- ter of the Milky Way compare to the distances of other stars from the center of the Milky Way? In Fig. 8 we show the distribution of the mean radial distances of stars from the Galac- tic center, based on the star count model of Bahcall & Soneira (1980). To represent the entire Galactic stellar population we include the disk (thin + thick) and spheroidal (bulge + halo) components. Using the current Solar distance from the center (R0 = 7.62 ± 0.32 kpc, Eisenhauer et al. 2005) and a disk scale length h = 3.0 ± 0.4 kpc (Gould et al. 1996), we estimate that the Sun lies farther from the Galactic center than 72+8−5% of the stars in the Galaxy. The uncertainty on the result comes from the 68% bounds of the total distribution, which come from the scale length uncertainty (±0.4 kpc). Fig. 8.— Mean stellar galactocentric radius distribution dN/dRGal. The solid curve rep- resents the sum of the disk (dashed line) and spheroidal (dotted line) stellar components. The 68% uncertainty of the total distribution is shown by the cross-hatched area. The Sun is farther from the Galactic center than 72+8 −5% of the stars in the Galaxy. 2.7. Host galaxy mass The mass of a star’s host galaxy may be correlated with parameters that have an in- fluence on habitability. For example, galaxy mass affects the overall metallicity distribu- tion that a star would find around itself — an 13 effect that would not show up in Fig. 3, which only shows the local metallicity distribution. The Milky Way is more massive than ∼ 99% of all galaxies — the precise fraction depends on the lower mass-limit chosen for an object to be classified as a galaxy, and the behaviour of the low-mass end of the galaxy mass function (Silk 2007). We are re- ferring here to the stellar mass, not the to- tal baryonic mass or the total mass. De- spite the Milky Way’s large mass compared to other galaxies, if most stars in the Uni- verse resided in even more massive galaxies, the Milky Way would be a rather low mass galaxy for a star to belong to. To estimate the fraction of all stars in galaxies of a given mass, we first estimate the distribution of galaxy masses by taking the K-band luminos- ity function of Loveday (2000) (K-band most closely reflects stellar mass since it is less sen- sitive than other bands to differences in stel- lar populations) and weighting it by luminos- ity. We convert this to stellar mass assum- ing a constant stellar-mass-to-light ratio of 0.5 (Bell & de Jong 2001). This function, plotted in Fig. 9, shows the amount of stellar mass contributed by galaxies of a given mass — or assuming identical stellar populations — the fraction of stars residing in galaxies of a given stellar mass. We estimate the K-band luminosity of the Milky Way by converting the published V- band magnitude of Courteau & van den Bergh (1999) to the K-band assuming the mean color of an Sbc spiral galaxy from the 2MASS Large Galaxy Atlas (Jarrett et al. 2003) and apply- ing the color conversion from (Driver et al. 1994). We then convert this to stellar mass using the same stellar-mass-to-light ratio used above, i.e., 0.5. In this way we estimate the stellar-mass content of the Milky Way to be 1010.55±0.16 = 3.6+1.5−1.1 × 10 10M (see also Flynn et al. 2006). Comparing this to the stellar masses of other galaxies (Fig. 9), Fig. 9.— Fraction of all stars that live in galaxies of a given mass, dN/dM (solid curve). The mass of the Sun’s galaxy is in- dicated by the “”. This distribution repre- sents the amount of stellar mass contributed by galaxies of a given mass. Approximately 77+11−14% of stars live in galaxies less massive than ours. The cross-hatched band shows the 1σ uncertainty associated with the un- certainty in the two Schechter function pa- rameters, α and L∗ (Loveday 2000; Schechter 1976). The dashed line shows the unweighted luminosity function (the number of galaxies per luminosity interval dNgals/dM) according to which the Milky Way is more massive than ∼ 99% of galaxies. we find that 77+11 −14% of stars reside in galaxies less massive than the Milky Way. 2.8. Host group mass The mass of a star’s host galactic group or galactic cluster may be correlated with parameters that have an influence on hab- itability. For example, group mass is cor- related with the density of the galactic en- 14 vironment (number of galaxies/Mpc3) which could, like galactocentric radius, be associ- ated with the dangers of high stellar densi- ties: “The presence of a giant elliptical at a distance of 50 kpc would have disrupted the Milky Way Galaxy, so that human beings (and hence astronomers) probably would not have come into existence.” (van den Bergh 2000). Our Local Group of galaxies seems rather typical (van den Bergh 2000) but we would like to quantify this. Proceeding simi- larly to our analysis of galaxy mass in Sect. 2.7, we ask: What fraction of stars live in galactic groups less massive than our Local Group? Figure 10 shows the luminosity- weighted (∼ stellar-mass-weighted) number density of galactic groups. The number dis- tribution and luminosity distribution of galac- tic groups is taken from the Two-degree Field Galaxy Redshift Survey Percolation-Inferred Galaxy Group (2PIGG) catalogue (Eke et al. 2004). It spans the range from weak groups to rich galaxy clusters. We estimated the stellar masses of the 2PIGG groups and Local Group galaxies (Courteau & van den Bergh 1999) by con- verting from the B-band assuming a constant stellar-mass-to-light ratio of 1.5 (Bell & de Jong 2001). This gives an estimated stellar mass of the local group of 1010.91±0.07 = 8.1+1.4−1.2 × 1010M. Figure 10 indicates that our Local Group is a typical galactic grouping for a star to be part of. Approximately 58±5% of stars live in galactic groups more massive than our Local Group. With respect to the mass of its galaxy and the mass of its galactic group, the Sun is a fairly typical star in the Universe. 3. Joint Analysis of 11 Solar Proper- ties 3.1. Solar χ2-analysis We would like to know if the solar prop- erties, taken as a group, are consistent with Fig. 10.— The dashed histogram shows the luminosity function of galactic groups (number of groups per interval of B-band lu- minosity). The solid histogram shows the luminosity-weighted group luminosity func- tion (approximately the fraction of stars which inhabit a group of given stellar mass). The horizontal axis has been converted to stel- lar mass assuming a constant B-band stellar- mass-to-light ratio of 1.5 (Bell & de Jong 2001). The “” shows the estimated mass of the Local Group (Courteau & van den Bergh 1999) and lies just below the median (vertical grey line). noise, i.e., are they consistent with the values of a star selected at random from our stel- lar distributions. We take a χ2 approach to answering this question. First we estimate the solar χ2, by adding in quadrature, for all eleven properties, the differences between the solar values and the median stellar values. We find: χ2 = N=11∑ i=1 (x,i − µ1/2,i) 2 σ268,i = 7.88+0.08−0.30 (1) 15 where i is the property index, N = 11 is the number of properties we are considering, µ1/2,i is the median of the ith stellar distribution and σ68,i is the difference between the median and the upper or lower 68% zone, depending on whether the solar value x,i is above or below the median. The uncertainty on χ2 is obtained using the uncertainties of x,i. Equation (1) can be improved upon by tak- ing into account: i) the non-Gaussian shapes of the stellar distributions and ii) the larger uncertainties of the medians of smaller sam- ples (our smallest sample is ∼ 100 stars). We employ a bootstrap analysis (Efron 1979) to randomly resample data (with re- placement) and derive a more accurate esti- mate of χ2. Because the bootstrap is a non- parametric method, the distributions need not be Gaussian. We obtain χ2 = 8.39 ± 0.96. Figure 11 shows the resulting solar chi-squared distri- bution. The median of this distribution is our adopted solar chi-squared value. Dividing our adopted solar chi-squared by the number of degrees of freedom gives our adopted reduced solar chi-squared value: χ2/11 = 0.76± 0.09 (2) The standard conversion of this into a probability of finding a lower chi-squared value (assuming normally distributed inde- pendent variables) yields: P (< χ2 = 8.39|N = 11) = 0.32 ± 0.09. (3) 3.2. Estimate of P (< χ2) To quantify how typical the Sun is with re- spect to our 11 properties, we compare the so- lar χ2(= 8.39) to the distribution of χ 2 values obtained from the other stars in the samples. We perform a Monte Carlo simulation (Metropolis & Ulam 1949) to calculate an es- timate of each star’s chi-squared value (χ2). Fig. 11.— Bootstrapped solar chi-squared distribution. The median of the distribution (white “”) is χ2 = 8.39± 0.96. This should be compared to the solar χ2 value from Eq. 1: 7.88+0.08−0.30 which is over-plotted (grey “” on dotted line). The histogram shown in Figure 12 is the resulting Monte Carlo stellar χ2 distribu- tion. Three standard chi-squared distribu- tions have been over-plotted for comparison (N = 10, 11, 12). The probability of finding a star with chi-squared lower than or equal to solar is: PMC(≤ χ 2  = 8.39|N = 11) = 0.29±0.11 (4) The Monte Carloed χ2 distribution has a simi- lar shape to the standard chi-squared distribu- tion function for N = 11, and thus both yield similar probabilities: PMC(≤ χ 2) = 0.29 ∼ P (≤ χ2) = 0.32 (Eqs. 3 and 4). The more appropriate Monte Carlo distribution has a longer tail, produced by the longer super- Gaussian tails of the stellar distributions. Table 2 summarizes our analysis for the So- lar χ2 values and the probabilities P (< χ 2 ). 16 Our simple χ2 = 7.88 estimate increased to 8.39 and the uncertainty increased by a factor of ∼ 3 after non-Gaussian and sam- ple size effects were included as additional sources of uncertainty. Our improved anal- ysis yields PMC(≤ χ 2 ), with a longer tail and brings the probability down from 0.32 ± 0.09 to 0.29 ± 0.11. If this value were close to 1, almost all other stars would have lower chi- squared values and we would have good reason to suspect that the Sun is not a typical star. However, this preliminary low value of 0.29 indicates that if a star is chosen at random, the probability that it will be more typical (∼ have a lower χ2 value) than the Sun (with re- spect to the eleven properties analysed here), is only 29± 11%. The details of our improved estimates of χ2 and P (< χ 2 ) can be found in the Appendix B. 4. Results Figure 13 shows four different representa- tions of our results. Panel (A) compares the solar values to each stellar distribution’s me- dian and 68% and 95% zones. The Sun lies be- yond the 68% zone for three properties: mass (95%), eccentricity (93%) and rotational ve- locity (88%). No solar property lies beyond the 95% zone. The histogram in panel (B) is the distribution of solar values in units of standard deviations: zi = x,i − µ1/2,i σ68,i (5) For each stellar property i, the Sun has a larger value than ni% of the stars. If the Sun were a randomly selected star, we would expect the percentages ni% to be scattered roughly evenly between 0% and 100%. When the ni% values are lined up in decreasing or- der (panel C), we expect them to be near the line given by: ni,expected% = ( 1− (i− 1/2) N ) × 100% (6) Fig. 12.— Stellar chi-squared distribution from our Monte Carlo simulation. PMC(< χ2 = 8.39) = 0.29 ± 0.11 (represented by the grey shade) is calculated integrating from χ2 = 0 to χ2 = χ2. For comparison, three χ 2 distribution-curves are over-plotted with 10, 11, and 12 degrees of freedom. The standard probability from the N = 11 curve yields: P (< χ2 = 8.39|N = 11) = 0.32 ± 0.09. and plotted in Panel C. Any anomalies would show up as points ‘’ significantly distant from the line. Panel (D) compares the percentages ni% of stars having sub-solar values (shown in Panel C) with the solar values expressed in units of standard deviations from each distribution’s median (shown in Panel B). If the stellar dis- tributions were perfect Gaussians, the trans- lation from zi to ni would be given by the cumulative Gaussian distribution (black line in Panel D). That the points lie along this line demonstrates that the approximation of our distributions as Gaussians is reasonable. Table 3 lists percentages ni% of stars for each property (as shown in Fig. 13). In the 17 Table 2 Summary of χ2 and P (< χ2) results. Analysis χ2 χ 2 /11 P (< χ 2 |N = 11) PMC(< χ 2 |N = 11) simple 7.88+0.08−0.30 (Eq. 1) 0.72 +0.01 −0.03 0.28 +0.01 −0.03 (Eq. B1) — improved 8.39± 0.96 0.76 ± 0.09 (Eq. 2) 0.32 ± 0.09 (Eq. 3) 0.29± 0.11 (Eq. 4) lower half of the table we list properties not included in this analysis because of correla- tions with properties that are included. Individual stellar uncertainties make the observed characteristic widths (σ68, column 5 of Table 1) larger than the widths of the intrinsic distributions. This broadening effect makes the Sun appear more typical than it re- ally is when σ68 and the individual stellar un- certainties (σ) are of similar size and the in- dividual stellar uncertainties are much larger than the solar uncertainty (σ). We estimate that our results are not significantly affected by this broadening effect. Our resulting probability of finding a star with a χ2 lower or equal to the solar value of 29± 11% (Eq. 4), is consistent with the prob- ability we would obtain if stellar multiplicity were included in our study. Using the vol- ume limited sample used for stellar mass in Section 2.1 (125 A1–M7 stars within 7.1 pc) the probability that a randomly selected star will be single is 52.8±4.5%, which means that ∼ half of stars are single while the other half have one or more companions. Including this in our bootstrap analysis and Monte Carlo simulations (see Appendix B.1) marginally in- creases the probability in Eq. 4 to 33 ± 11%. If the multiplicity data for 246 G dwarfs from Duquennoy & Mayor (1991) is used instead — the probability that a randomly selected G dwarf will be single is 37.8 ± 2.9% — then the probability in Eq. 4 would increase to 34±11%. The inclusion of stellar multiplicity marginally increases our reported probability. In Figures 6 and 7 of Radick et al. (1998), the Sun’s short-term variability as a function of average chromospheric activity, appears ∼ 1σ low, compared to a distribution of 35 F3–K7 Sun-like stars (Lockwood et al. 1997). Lockwood et al. (2007) suggest that the Sun’s small total irradiance variation compared to stars with similar mean chromospheric activ- ity, may be due to their limited sample and the lack of solar observations out of the Sun’s equatorial plane. We do not include short or long term variability (chromospheric or pho- tometric) in Table 3 because of the small size of the Lockwood et al. (2007) sample. We also do not include the chromospheric index R′HK (see Table 3, bottom panel) as one of our 11 properties because of its correlation with the chromospheric ages of our sample. 18 Fig. 13.— Various representations of our main results. A: Solar values of eleven properties com- pared to the distribution for each property Each distribution’s median value is indicated by a small filled circle. The dark and light grey shades represent the 68% and 95% zones respectively. B: Histogram of the number of properties as a function of the number of standard deviations the solar value is from the median of that property. The grey curve is a Gaussian probability distribution normalised to 11 parameters. C: Percentage ni% of stars with sub-solar values as a function of property. The average signal expected from a random star is shown by the solid line (see Sec. 4). D: Percentage ni% of stars with sub-solar values as a function of the number of standard devia- tions the solar value is from the median of that property. The solid curve is a cumulative Gaussian distribution — if every sample were a Gaussian distribution, every solar dot would sit exactly on the line. Just as in (C), the dashed lines encompass the 68% and 95% zones. Similar to the results from Figure 12, these four panels indicate that the Sun is a typical star. 19 Table 3 Summary of How the Sun Compares to Other Stars (see Fig. 13) Parameter Fig. ni% Level of Anomaly Mass 1 95± 2% of nearby stars are less massive than the Sun. Age 2 53± 2% of stars in the thin disk of the Galaxy are older than the Sun. [Fe/H] 3 65± 2% of nearby stars are more iron-poor than the Sun. [C/O] 4A 81± 3% of nearby stars have a higher C/O ratio than the Sun. [Mg/Si] 4B 66± 3% of nearby stars have a higher Mg/Si ratio than the Sun. v sin i 5 83± 7% of nearby Sun-like-mass stars rotate faster than the Sun. e 6 93± 1% of nearby stars have larger galactic orbital eccentricities than the Sun. Zmax 7 59± 3% of nearby stars reach farther from the Galactic plane than the Sun. RGal 8 72 +8 −5% of stars in the Galaxy are closer to the galactic center than the Sun. Mgal 9 77 +11 −14% of stars in the Universe are in galaxies less massive than the Milky Way. Mgroup 10 58± 5% of stars in the Universe are in groups more massive than the local group. Properties not included in the analysis because they are correlated with the selected 11 parameters Mass: IMFStellar 1 94± 2% of nearby stars are less massive than the Sun. Age: SFRCosmic 2 86± 5% of stars in the Universe are older than the Sun. Agea — 55± 2% of nearby Sun-like-mass stars are older than the Sun. [Fe/H]b — 56± 5% of nearby stars are more iron-poor than the Sun. v sin ic — 92± 5% of nearby Sun-like-mass stars rotate faster than the Sun. logR′HK d — 51± 2% of nearby FGKM stars are more chromospherically active. [O/Fe] — 75± 3% of nearby stars have a lower O/Fe ratio than the Sun. Rmin — 91± 1% of nearby stars get closer to the Galactic center. vLSR — 93± 1% of nearby stars have smaller velocity with respect to the LSR. |U | — 75± 1% of nearby stars have larger absolute radial velocity. |V | — 82± 1% of nearby stars have larger absolute tangential velocity. |W | — 58± 1% of nearby stars have larger absolute vertical velocity. a1126 stars (A5–K2) from Nordstro¨m et al. (2004). b91 stars (GK) from Favata et al. (1997). c590 stars (F8–K2) from Nordstro¨m et al. (2004). d866 stars (FGKM) from Wright et al. (2004). 20 5. Discussion and Interpretation The probability PMC(≤ χ 2 ) = 0.29 ± 0.11 classifies the Sun as a typical star. How ro- bust is this result? The probability of finding a star with a chi-squared lower than or equal to χ2, depends on the properties selected for the analysis (see problem iii of Section I). For example, if we had chosen to consider only mass and eccentricity data, this analy- sis would yield PMC(χ 2 ≤ χ2) = 0.94 ± 0.4, i.e., the Sun would appear mildly (∼ 2σ) anomalous. If on the other hand, we had cho- sen to remove mass and eccentricity from the analysis, we would obtain PMC(χ 2 ≤ χ2) = 0.07 ± 0.04, which is anomalously low. The most common cause of such a result is the over-estimation of error bars. The next most common cause is the preselection of properties known to have ni% ∼ 50%. Gustafsson (1998) discussed the atypically large solar mass, and proposed an anthropic explanation — the Sun’s high mass is prob- ably related to our own existence. He sug- gested that the solar mass could hardly have been greater than ∼ 1.3M since the main se- quence lifetime of a 1.3M star is ∼ 5 billion years (Clayton 1983). He also discussed how the dependence of the width of the circum- stellar habitable zone on the host star’s mass probably favours host stars within the mass range 0.8–1.3M. Our property selection criteria is to have the largest number of maximally independent properties that have a plausible correlation with habitability and, ones for which a rep- resentative stellar sample could be assembled. Our joint analysis does not weight any pa- rameter more heavily than any other. If the only properties associated with habitability are mass and eccentricity then we have diluted a ∼ 2σ signal that would be consistent with Gustafsson’s proposed anthropic explanation. Our analysis points in another direction. If mass and eccentricity were the only properties associated with habitability, then the solar values for the remaining 9 properties would be consistent with noise. However, a joint analy- sis of just the remaining 9 properties produces a χ2,9 = 3.6 ± 0.4 and the anomalously low probability: P (≤ χ2,9) = 0.07 ± 0.04, which suggests that the 9 properties are unlikely to be the properties of a star selected at random with respect to these properties. The χ2 fit of the 11 points in Panel C of Fig. 13 to the diagonal line yields a fit that is substantially better then the fit of the re- maining 9 properties to Eq. 7 with N = 9. In other words, the joint analyis suggests that although mass and eccentricity are the most anomalous solar properties, it is unlikely that they are associated with habitability, because without them, it is unlikely that the remaining solar properties are just noise. Thus, the Sun, despite its mildly (∼ 2σ) anomalous mass and eccentricity, can be considered a typical, ran- domly selected star. There may be stellar properties crucial for life that were not tested here. If we have left out the most important properties, with re- spect to which the Sun is atypical, then our Sun-is-typical conclusion will not be valid. If we have sampled all properties associated with habitability, our Sun-is-typical result suggests that there are no special require- ments on a star for it to be able to host a planet with life. 6. Conclusions We have compared the Sun to representa- tive stellar samples for eleven properties. Our main results are: • Stellar mass and Galactic orbital eccen- tricity are the most anomalous proper- ties. The Sun is more massive than 95± 2% of nearby stars and has a Galactic orbital eccentricity lower than 93 ± 1% 21 FGK stars within 40 pc. • Our joint bootstrap analysis yields a so- lar chi-squared χ2 = 8.39 ± 0.96 and a solar reduced chi-squared χ2/11 = 0.76 ± 0.09. The probability of find- ing a star with a chi-squared lower than or equal to solar PMC(≤ χ 2  = 8.39 ± 0.96) = 0.29 ± 0.11. To our knowledge, this is the most compre- hensive and quantitative comparison of the Sun with other stars. We find that taking all eleven properties together, the Sun is a typi- cal star. This finding is largely in agreement with Gustafsson (1998), however our results undermine the proposition that an anthropic explanation is needed for the comparatively large mass of the Sun. Further work could encompass the inclu- sion of other properties potentially associ- ated with habitability. Another improvement would come when larger stellar samples be- come available for which all properties could be derived, instead of using different samples for different properties as was done here. In addition, research in the molecular evolution that led to the origin of life may, in the future, be able to provide more clues as to which stel- lar properties might be associated with our existence on Earth, orbiting the Sun. Acknowledgments: We would like to thank Charles Jenkins for clarifying discussions of statistics, particularly on how to include stel- lar multiplicity, and Martin Asplund and Jorge Mele´ndez for discussions of elemental abundances. JAR acknowledges an RSAA PhD research scholarship. MP acknowledges the financial support of the Australian Re- search Council. EG acknowledges the finan- cial support of the Finnish Cultural Founda- tion. 22 A. Property-correlations The χ2 formalism and the use of the χ2-distribution to obtain P (< χ2|N), — improved using Monte-Carlo simulations in Section 3.2 to obtain PMC(≤ χ 2 ) — assumes that each parameter is independent of the others. In selecting our 11 properties we have selected properties which are maximally independent based on plotting property 1 vs property 2 for the same stars. We show seven such plots in this Appendix. If there are correlations between the analysed properties, then the number of degrees of freedom N could drop from 11 to ∼ 10.5 (see Fig. 12). Some properties have been excluded from the analysis due to a correlation with another property in the analysis. A.1. Elemental ratios Fig. 14.— Carbon to oxygen ratio [C/O] versus magnesium to silicon ratio [Mg/Si] of 176 FG stars with abundances for these elements (Reddy et al. 2003). In Figure 4 (bottom panels) we showed that the [C/O] and [Mg/Si] distributions are largely independent of [Fe/H]. Here we show that these distributions are also largely independent of each other. Note that in this comparison we only use the data from Reddy et al. (2003), since it is the largest available sample with C, O, Mg and Si abundances. A.2. Mass, age and rotational velocity In Figure 15 we show four correlation plots for mass, chromospheric age, rotational velocity and v sin i. We use the stars common to both Wright et al. (2004) and Valenti & Fischer (2005) for which these observables are available. 23 Fig. 15.— Correlation plots between various properties. For all four panels we use the stars common to both Wright et al. (2004) and Valenti & Fischer (2005). Panel (A): mass vs rotational velocity v sin i for 713 FGK stars. This panel shows the degree of correlation between mass and v sin i. See Gray (2005) for a stronger correlation between these two variables when a larger mass range and more active stars are kept in the sample. To minimize the effect of this correlation on our analysis, we restrict the range of mass in Fig. 5 to 0.9 to 1.1M. Panel (B): chromospheric age versus v sin i for 641 FGK stars. The lack of correlation between chromospheric determined ages and rotational velocities is shown. Panel (C): no strong correlation between mass and chromospheric age for 639 FGK stars. Panel (D): the ages of 637 stars determined by the chromospheric method versus their ages from the isochrone method. 24 A.3. Galactic orbital parameters The Galactic orbital eccentricity (e) and the magnitude of the galactic orbital velocities with respect to the local standard of rest (vLSR) are strongly correlated (see Fig. 6 in Sec. 2.6). We selected e instead of vLSR because of its near independence of the maximum height above the galactic plane (Zmax). Fig. 16.— Left panel: Galactic orbital eccentricity e versus Zmax for 1987 FGK stars within 40 pc (Nordstro¨m et al. 2004). The orbital eccentricity is not correlated with Zmax. Right panel: vLSR versus Zmax for the same stars. Because vLSR is more strongly correlated with Zmax than eccentricity, eccentricity has been selected for the joint analysis instead of vLSR. As in Fig. 4, the contours correspond to 38%, 68%, 82% and 95%. B. Improved Estimates of χ2 and P (< χ 2 ) In Section 3.2, with 11 degrees of freedom, the reduced chi-squared from Equation 4 is χ2 / 11 = 0.72+0.01−0.03. Since χ 2  / 11 < 1, the Sun’s properties are consistent with the Sun being a randomly selected star. To improve on this preliminary analysis (but with a similar conclusion), as mentioned in Section 3.2, we employ a bootstrap analysis (Efron 1979) to randomly resample data (with replacement) and derive a more accurate estimate of χ2. Because the bootstrap is a non-parametric method, the distributions need not be Gaussian. For every iteration, each parameter’s stellar distribution is randomly resampled and a χ2 value is calculated using Eq. (1). The uncertainties σ,i of the solar values x,i are also included in the bootstrap method: for every iteration, the Solar value for each parameter is replaced in Eq. (1) by a randomly selected value from a normal distribution with median µ1/2,i = x,i and standard 25 deviation σ,i. The process was iterated 100,000 times, although the resulting distribution varies very little once the number of iterations reaches ∼ 10, 000. The median of this distribution and the error on the median yields our improved value for the reduced χ2 (Fig. 11). The uncertainty of the median of each re-sampled distribution varies inversely proportionally to the square root of the number of stars in the distribution, ∆µ1/2,i ∝ 1/ √ N,i. In other words, median values are less certain for smaller samples and this uncertainty is included in our improved estimate of χ2, and its uncertainty. We find the probability of finding a star with a χ2 value lower than the solar χ 2 , for N = 11 degrees of freedom in the standard way (Press et al. 1992) and obtain: P (< χ2 = 7.88 +0.08 −0.30|11) = 0.28 +0.01 −0.03 (B1) To improve our estimate of the probability of finding a star with lower chi-squared value than the Sun, we perform a Monte Carlo simulation (Metropolis & Ulam 1949) to calculate an estimate of each star’s chi-squared value (χ2). For every iteration, we randomly select a star from each stellar distribution. We then calculate its χ2 value by replacing the solar value x,i with that star’s value x,i in Eq. (1). This process was repeated 100,000 times to create our Monte Carloed stellar chi-squared distribution. Stars were randomly selected with replacement, thus the simulated χ2 distribution accounts for small number statistics and non-Gaussian distributions. The probability of finding a star with chi-squared lower than or equal to solar is PMC = 0.29 ± 0.11. The results of our analysis for the Solar χ2 values and the probabilities P (< χ 2 ) are summarized in Table 2 B.1. Addition of a discrete parameter In Section 4 we discuss the addition of stellar multiplicity to our analysis. Since stellar multiplic- ity cannot easily be approximated by a one-sided Gaussian (particularly because the Sun is on the edge of the distribution, i.e., it is of multiplicity one), we modified our Monte Carlo procedure to include this discrete parameter. The likelihood of observing a particular χ2 for the 11 parameters is exp ( − 1 2 11∑ i=1 χ2i ) . (B2) We take the probability p(1) of a star being a single star, to be 53.8 ± 4.5%, obtained from our sample of nearby stars (Sec. 2.1). The likelihood L of observing a particular χ2 and p(1) is the product L = p(1) exp ( − 1 2 11∑ i=1 χ2i ) . (B3) Taking logarithms we can then compute the distribution of the statistic S, where S = ln p(1)− 1 2 11∑ i=1 χ2i . (B4) The distribution of S allows us to obtain the results for multiplicity reported at the end of Section 4. 26 REFERENCES Alle`gre, C. J., Manhe`s, G., & Go¨pel, C. 1995, Geochim. Cosmochim. Acta, 59, 1445 Allende Prieto, C. 2006, ArXiv Astrophysics e-prints Allende Prieto, C., Lambert, D. L., & As- plund, M. 2001, ApJ, 556, L63 Asplund, M., Grevesse, N., & Sauval, A. J. 2005, in ASP Conf. Ser. 336: Cosmic Abun- dances as Records of Stellar Evolution and Nucleosynthesis, 25 Bahcall, J. N., & Soneira, R. M. 1980, ApJS, 44, 73 Bell, E. F., & de Jong, R. S. 2001, ApJ, 550, 212 Bensby, T., & Feltzing, S. 2006, MNRAS, 367, 1181 Bensby, T., Feltzing, S., Lundstro¨m, I., & Ilyin, I. 2005, A&A, 433, 185 Bertelli, G., & Nasi, E. 2001, AJ, 121, 1013 Carter, B. 1983, Philos. Trans.R. Soc. Lon- don, A, 310 Clayton, D. D. 1983, Principles of stellar evo- lution and nucleosynthesis (Chicago: Uni- versity of Chicago Press, 1983) Courteau, S., & van den Bergh, S. 1999, AJ, 118, 337 Dehnen, W., & Binney, J. J. 1998, MNRAS, 298, 387 Driver, S. P., Phillipps, S., Davies, J. I., Mor- gan, I., & Disney, M. J. 1994, MNRAS, 268, 393 Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485 Edvardsson, B., Andersen, J., Gustafsson, B., Lambert, D. L., Nissen, P. E., & Tomkin, J. 1993a, A&A, 275, 101 —. 1993b, A&AS, 102, 603 Efron, B. 1979, The Annals of Statistics, 7, 1 Eisenhauer, F. et al. 2005, ApJ, 628, 246 Eke, V. R. et al. 2004, MNRAS, 355, 769 Favata, F., Micela, G., & Sciortino, S. 1996, A&A, 311, 951 —. 1997, A&A, 323, 809 Feltzing, S., Holmberg, J., & Hurley, J. R. 2001, A&A, 377, 911 Flynn, C., Holmberg, J., Portinari, L., Fuchs, B., & Jahreiß, H. 2006, MNRAS, 372, 1149 Flynn, C., Sommer-Larsen, J., & Christensen, P. R. 1996, MNRAS, 281, 1027 Fuhrmann, K. 2008, MNRAS, 384, 173 Gonzalez, G. 1999a, MNRAS, 308, 447 —. 1999b, Astronomy and Geophysics, 40, 25 Gonzalez, G., Brownlee, D., & Ward, P. 2001, Icarus, 152, 185 Gould, A., Bahcall, J. N., & Flynn, C. 1996, ApJ, 465, 759 Gray, D. F. 2005, The Observation and Analysis of Stellar Photospheres (The Observation and Analysis of Stellar Photo- spheres, 3rd Edition, by D.F. Gray. ISBN 0521851866. http://www.cambridge.org/us/ /catalogue/catalogue.asp?isbn=0521851866. Cam- bridge, UK: Cambridge University Press, 2005.) Grether, D., & Lineweaver, C. H. 2006, ApJ, 640, 1051 —. 2007, ApJ, 669, 1220 27 Gustafsson, B. 1998, Space Science Reviews, 85, 419 Gustafsson, B., Karlsson, T., Olsson, E., Ed- vardsson, B., & Ryde, N. 1999, A&A, 342, 426 Henry, T. J. 2006, RECONS database Hernandez, X., Valls-Gabaud, D., & Gilmore, G. 2000, MNRAS, 316, 605 Hogg, D. W., Blanton, M. R., Roweis, S. T., & Johnston, K. V. 2005, ApJ, 629, 268 Hopkins, A. M., & Beacom, J. F. 2006, ApJ, 651, 142 Jarrett, T. H., Chester, T., Cutri, R., Schnei- der, S. E., & Huchra, J. P. 2003, AJ, 125, 525 Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108 Kroupa, P. 2002, Science, 295, 82 Kroupa, P., & Weidner, C. 2005, in ASSL Vol. 327: The Initial Mass Function 50 Years Later, ed. E. Corbelli, F. Palla, & H. Zin- necker, 175 Kuchner, M. J., & Seager, S. 2005, ArXiv As- trophysics e-prints Lineweaver, C. H., Fenner, Y., & Gibson, B. K. 2004, Science, 303, 59 Lockwood, G. W., Skiff, B. A., Henry, G. W., Henry, S., Radick, R. R., Baliunas, S. L., Donahue, R. A., & Soon, W. 2007, ApJS, 171, 260 Lockwood, G. W., Skiff, B. A., & Radick, R. R. 1997, ApJ, 485, 789 Loveday, J. 2000, MNRAS, 312, 557 Metropolis, N., & Ulam, S. 1949, Journal of the American Statistical Association, 44, 335 Metzger, M. R., Caldwell, J. A. R., & Schechter, P. L. 1998, AJ, 115, 635 Meyer, S. L. 1975, Data Analysis for Scien- tists and Engineers (Data Analysis for Sci- entists and Engineers, by Stuart L. Meyer p. 186. ISBN 0471599956. NY, USA: John Wiley & Sons, 1975.) Nordstro¨m, B. et al. 2004, A&A, 418, 989 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in FORTRAN. The art of scientific computing (Cambridge: University Press, —c1992, 2nd ed.) Radick, R. R., Lockwood, G. W., Skiff, B. A., & Baliunas, S. L. 1998, ApJS, 118, 239 Ramı´rez, I., Allende Prieto, C., & Lambert, D. L. 2007, A&A, 465, 271 Reddy, B. E., Tomkin, J., Lambert, D. L., & Allende Prieto, C. 2003, MNRAS, 340, 304 Reid, I. N. 2002, PASP, 114, 306 Reid, I. N., Turner, E. L., Turnbull, M. C., Mountain, M., & Valenti, J. A. 2007, ApJ, 665, 767 Rocha-Pinto, H. J., Maciel, W. J., Scalo, J., & Flynn, C. 2000a, A&A, 358, 850 Rocha-Pinto, H. J., Scalo, J., Maciel, W. J., & Flynn, C. 2000b, ApJ, 531, L115 Schechter, P. 1976, ApJ, 203, 297 Silk, J. 2007, Astronomy and Geophysics, 48, 30 Soderblom, D. R. 1983, ApJS, 53, 1 —. 1985, AJ, 90, 2103 Truran, Jr., J. W., & Heger, A. 2005, Origin of the Elements (Meteorites, Comets and Planets: Treatise on Geochemistry, Volume 1) 28 Valenti, J. A., & Fischer, D. A. 2005, ApJS, 159, 141 van den Bergh, S. 2000, The Galaxies of the Local Group (The galaxies of the Local Group, by Sidney Van den Bergh. Pub- lished by Cambridge, UK: Cambridge University Press, 2000 Cambridge Astro- physics Series Series, vol no: 35, ISBN: 0521651816.) Wright, J. T., Marcy, G. W., Butler, R. P., & Vogt, S. S. 2004, ApJS, 152, 261 This 2-column preprint was prepared with the AAS LATEX macros v5.2. 29 The Galaxy Disk in Cosmological Context Proceedings IAU Symposium No. 254, 2008 J. Andersen (Chief Editor), J. Bland-Hawthorn & B. Nordstro¨m c© 2008 International Astronomical Union Dynamical effects of the long bar in the Milky Way Esko Gardner1, Kimmo A. Innanen2 and Chris Flynn1 1Tuorla Observatory, Department of Physics and Astronomy, University of Turku Va¨isala¨ntie 20, FI-21500 Piikkio¨, Finland e-mail: esgard@utu.fi 2Dept. of Physics and Astronomy York University, Toronto, Canada Abstract. We examine the dynamical effects on disk stars of a “long bar” in the Milky Way by inserting a triaxial rotating bar into an axisymmetric disk+bulge+dark halo potential and integrating 3-D orbits of 104 tracer stars over a period of 2 Gyr. The long bar has been detected via “clump giants” in the IR by Lo´pez-Corredoira et al. (2007), and is estimated to have semi-major axes of (3.9× 0.6× 0.1) kpc and a mass of 6× 109 M. We find such a structure has a slight impact on the inner disk-system, moving tracers near to the bar into the bar-region, as well as into the bulge. These effects are under continuing study. Keywords. Galaxy: kinematics and dynamics, Galaxy: center, Galaxy: disk 1. Introduction Our initial interest was to study a barred galaxy potential, and especially the resonance effects which might lead to vertical ejection into the halo of globular clusters formed in the disk as suggested by Innanen (2007). We took a triaxial ellipsoid, a Ferrers’ potential, to model the bar (see Pfenniger (1984) for a similar application of the potential as a bar). We decided to look at the effects of the “long bar” described by Lo´pez-Corredoira et al. (2007), taking their parameters directly, and applying it to our system. The dynamical effects of the bar are quite interesting, and it should be possible to constrain the bar parameters with analysis of these effects. We also present a brief discussion on the scale-length of the Galactic disk, and a revised version of the Flynn et al. (1996) Galactic disk model with a shorter scale-length of 3 kpc rather than 4.4 kpc, to be more consistent with recent observations. 2. Scale-length of the Milky Way disk Our motivation in constructing a disk model with a revised disc scale-length stems from recent infrared observations, which show a considerably shorter scale-length for the Galactic disk than optical observations. The disk scale-length (hR) adopted in Flynn et al. (1996) was 4.4 kpc, based on the kinematic determination (using velocity dispersions of old disk red giants) by Lewis & Freeman (1989), who obtained a scale-length of hR = 4.4 ± 0.3 kpc. It has been long known that the scale-length measured in the infrared is considerably shorter than in the optical, Pioneer 10 data giving hR = 5.5 ± 1.0 kpc already in the 1980s (van der Kruit (1986)). More recently near-IR data in J and K with DENIS give hR = 2.3±0.1 kpc and 2MASS of K giants gives hR = 3.34±0.1 kpc (Lo´pez- Corredoira et al. (2002), Ruphy et al. (1996), Ojha (2001)). These latter surveys probe 1 II 2 Gardner et al. the scale-length of emitted light, mainly by giant stars, rather than probing directly the scale-length of the mass, which is what we are actually interested for modeling the disk’s potential. M dwarfs do trace the mass directly: star-counts with HST in the 1990s of disk M dwarfs indicate a disk scale-length of 3.0± 0.4 kpc (Gould et al. (1997)). On the other hand, local kinematics of disk stars, which should also be mostly sensitive to the mass rather than the light distribution, give a range of possible scale-lengths 1.7−2.9 kpc (Bienayme & Sechaud (1997)). Very recently, F to K type dwarfs, for which distances and metallicities have been determined in the huge numbers in 6500 square degrees of sky and probing to distances of up to 20 kpc above and below the disk near the Sun using SDSS data, give a scale-length of 2.6± 0.5 kpc (Juric´ et al. (2008)). 3. The disk model The Flynn et al. (1996) model of the Galactic disk has a hR = 4.4 kpc, and is built from three superimposed Miyamoto-Nagai potentials (with different linear scales) to obtain the exponentially falling density profile. In light of recent observations, we would like a shorter disk scale-length, and modified the Flynn et al. (1996) model to a disk with hR = 3.0 kpc. The local volume density local surface density of the disk are as in the 1996 model (i.e. ρ0 ≈ 0.1M/pc 3 and Σ ≈ 50M/pc 2), and are consistent with local measurements (Holmberg & Flynn (2000)). We adopt a solar Galactocentric distance of R = 8 kpc. Figure 1 shows the surface density profile of the disk model, with the dashed line indicating a scale-length of 3 kpc. This is a good fit over a wide range of Galactocentric radii (2 to 15 kpc). Note that the disk truncates strongly at approximately 18 kpc. Figure 1. The surface density of the disk component as a function of Galactocentric radius. The dashed line corresponds to an exponential density falloff, 3 kpc, which is a good fit to the model over a wide range of radii. Note that the density truncates strongly at 18 kpc. 4. Inserting the long bar We insert a triaxial bar potential into the central regions of the model. We chose a uniform density Ferrers’ potential (see eg. Chandrasekhar (1969), or Binney & Termaine). The linear semi-major axes (3.9:0.6:0.1 kpc), and mass (6 × 109M) were chosen from Dynamical effects of the long bar in the Milky Way 3 Lo´pez-Corredoira et al. (2007), who used “red clump” stars detected in a 10 degree wide region observed in the near IR along the Galactic plane. Note that the long bar is very thin vertically, a property which leads to some interesting effects on orbits of disk stars which come near it. The bar is set to rotate around the z-axis with with a pattern speed of 50 km/s/kpc, giving it a speed of 200 km/s at its tip. This is the speed suggested by OLR (Dehnen (2000)), although we have tried a range of pattern speeds from 0 to 75 km/s/kpc, without substantially altering the conclusions of the paper. The initial phase angle of the bar is as observed by Lo´pez-Corredoira et al. (2007), i.e. 43 degrees. An important point to note is that we aim to investigate the effects of the “long bar”, rather than the smaller “COBE bar” (Bissantz & Gerhard (2002)), which has dimensions (3.5:1.4:1) kpc, and a position angle of 22 degrees is rather dissimilar to the long bar. 5. Disk simulation We set up an exponential distribution of tracer stars with the same scale-length as the disk. We chose U, V , and W velocities from a normal distribution with a standard deviation of 10 km/s around the local standard of rest (i.e. the circular velocity) in the total potential, so that the stars are very nearly on exact circular orbits. Thus the disk is initially inherently cold, and aids in following the secular development of the long bar’s kinematic effects. We integrated the orbits of 104 objects in steps of 104 years, with a total time of 2 Gyr, in both the barred and unbarred systems, using the long bar with the parameters advocated by Lo´pez-Corredoira et al. (2007). We also examined the effects of a thicker, less dense long bar, with an axial ratio of (3.9:2:1) kpc. 6. Results and discussion We present here a preliminary look at the results. In Figures 2, and 3, we show the final distribution of the objects, after 2 Gyr of orbital integrations, for the cases of “no bar”, a long bar of high density (6 M/pc 3) and a long bar with a much lower density (0.2 M/pc 3). The initial (orange), and final radial distributions of the objects are shown in Figure 2. The vertical distribution is shown in Figure 3. In both, the “no bar” case (cyan), shows little secular evolution, indicating that the disk has been set up in an sta- ble manner. In the denser long bar (magenta), the bar is injecting the inner regions with disk stars from approximately 1 kpc beyond the tip of the bar. Within the bar region, significant modification of the distribution profile takes place, as one would expect. This is seen as a knee in the distribution profile at logR = 3.0− 3.6. The less dense, thick bar (blue), has much smaller effects, with a small knee forming around 1 kpc, which coincides with the spherical area inside the bar. The bars have almost no effect on the z-distribution of the disk, there is a slight shift towards higher values in the thin case (magenta), these are mostly tracers that orbitally evolve towards the center of the system, ie. into the bulge. There is no significant ejection to large heights above the disk. Acknowledgements: We thank Martin Lo´pez-Corredoira for his useful comments on the manuscript. EG acknowledges the financial support of the Finnish Cultural Foundation, and the Finnish Graduate School in Astronomy and Space Physics. CF and KI are very grateful to the Academy of Finland for financial support. 4 Gardner et al. Figure 2. Histograms of Galactocentric radii of the tracer stars, for the initial distribution (orange), and after all three simulations have run (2 Gyr). The vertical line marks the position of the Sun at 8 kpc from the Galactic center. The “no bar” case (cyan) shows no secular evolution, indicating the initial setup in the disk is stable. Cases with the bar turned on show some secular evolution. Analysis shows that the long thin bar (magenta) has slight dynamical effects on disk stars which pass near it, injecting the inner disk, and bulge, with stars that have orbits near the tip of the bar. The thicker bar exibhits a feature around 1 kpc. <0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Log(z) 0 500 1000 1500 2000 2500 C o u n t Starting No Bar Thin Bar Thick Bar Figure 3. Histograms of the vertical height of the tracer stars, initially, and at the end of 2 Gyrs integration, for the “no bar” case (cyan), the long bar and the thicker long bar cases (magenta, blue). No secular evolution is seen in the distributions with the bar, although there is a slight trend to higher z-values in the thin bar case. References Bienayme, O. & Sechaud, N. 1997, A&A, 323, 781 Binney, J. & Tremaine, S. 2008, Galactic Dynamics 2 ed., Princeton University Press Bissantz, N. & Gerhard, O. 2002, MNRAS, 330, 591 Chandrasekhar, S. 1969, Ellipsoidal Figures of Equilibrium, Yale University Press Dehnen, W. 2000, AJ, 119, 800 Dynamical effects of the long bar in the Milky Way 5 Flynn, C., Sommer-Larsen, J., & Christensen, P.R. 1996, MNRAS, 281, 1027 Gould, A., Bachall, J. N. & Flynn, C. 1997, ApJ, 482, 913 Holmberg, J. and Flynn, C. 2000, MNRAS,313,209 Innanen, K. A. 2007, CeMDA, 99, 161 Juric´, M., Ivezic´, Zˇ., Brooks, A., Lupton, R. H., Schlegel, D., Finkbeiner, D., Padmanabhan, N., Bond, N., Sesar, B., Rockosi, C. M. and Knapp, G. R., Gunn, J. E., Sumi, T., Schneider, D. P. and Barentine, J. C., Brewington, H. J., Brinkmann, J., Fukugita, M., Harvanek, M., Kleinman, S. J., Krzesinski, J., Long, D., Neilsen, Jr., E. H., Nitta, A., Snedden, S. A. & York, D. G. 2008, ApJ, 673, 864 van der Kruit, P. C. 1986, A&A, 157, 230 Lewis, J. R. & Freeman, K. C. 1989, AJ, 97, 139 Lo´pez-Corredoira, M., Cabrera-Lavers, A., Garzo´n, F.& Hammersley, P. L. 2002 A&A, 394, 883 Lo´pez-Corredoira, M., Cabrera-Lavers, A., Mahoney, T. J., Hammersley, P. L., Garzo´n, F. & Gonza´lez-Ferna´ndez, C. 2007, AJ, 133, 154 Ojha, D. K. 2001, MNRAS, 322, 426 Pfenniger, D. 1984, A&A, 134, 373 Ruphy, S., Robin, A. C., Epchtein, N., Copet, E., Bertin, E., Fouque, P. & Guglielmo, F. 1996, A&A, 313, 21 Walterbos, R. A. M. & Kennicutt, Jr., R. C. 1988, A&A, 198, 61 Mon. Not. R. Astron. Soc. 000, 1–9 (2010) Printed 29 June 2010 (MN LATEX style file v2.2) Probing the Galaxy’s bars via the Hercules stream Esko Gardner1 and Chris Flynn1 1Tuorla Observatory, Department of Physics and Astronomy, University of Turku, Va¨isa¨la¨ntie 20,FI-21500, Piikkio¨, Finland April 2010 ABSTRACT It has been suggested that a resonance between a rotating bar and stars in the solar neighbourhood can produce the so called ‘Hercules stream’. Recently, a second bar may have been identified in the Galactic centre, the so called ‘long bar’, which is longer and much flatter than the traditional Galactic bar, and has a similar mass. We looked at the dynamical effects of both bars, separately and together, on orbits of stars integrated backwards from local position and velocities, and a model of the Galactic potential which includes the bars directly. Both bars can produce Hercules like features, and allow us to measure the rotation rate of the bar(s). We measure a pattern speed, for both bars, of 1.87 ± 0.02 times the local circular frequency. This is on par with previous measurements for the Galactic bar, although we do adopt a slightly different Solar motion. Finally, we identify a new kinematic feature in local velocity space, caused by the long bar, which is tempting to identify with the high velocity ‘Arcturus’ stream. Key words: Galaxy: centre – Galaxy: kinematics and dynamics – solar neighbour- hood – Galaxy: structure. 1 INTRODUCTION The inner regions of the Milky Way contain a non- axisymmetric structure sometimes referred to as the triax- ial (or boxy) bulge and/or the Galactic bar. Primary ev- idence for this structure are the COBE/Diffuse Infrared Background Experiment (DIRBE) near infrared (NIR) lu- minosity maps of the inner Galaxy (Binney et al. 1997; Bissantz & Gerhard 2002), but it has also been mapped via star counts using Two Micron All Sky Survey (2MASS) (Lo´pez-Corredoira et al. 2005), 2MASS and Optical Grav- itational Lensing Experiment II (OGLE-II) (Vanhollebeke et al. 2009). Typical parameters for the triaxial bulge/bar are axis ratios of 1.0:0.4:0.3, that it lies at an angle of be- tween 10◦ and 40◦ to the line of sight to the Galactic centre (see e.g. table 1 of Vanhollebeke et al. 2009) and that its mass is of order 1 × 1010 M  , (Zhao 1996; Weiner & Sell- wood 1999). It has recently been proposed (Benjamin et al. 2005) that the inner Galaxy also contains a ‘long bar’, discovered as an over-density in star count data at 4.5 µm taken in the Galactic Legacy Infrared Mid-Plane Survey Extraordi- naire (GLIMPSE) survey of the Galactic centre made with Spitzer. Benjamin et al. (2005) measure the half-length of this bar as 4.4 ± 0.5 kpc, and lying at a position angle of 44± 10◦ to the line of sight to the Galactic centre. Only the  E-mail:esgard@utu.fi closer end of the bar was detected, as the far end would lie at apparent magnitudes fainter than the completeness limit of the survey. The long bar has been confirmed by Lo´pez- Corredoira et al. (2007), who detect it via red clump stars identified from the 2MASS survey. They find a half-length of 3.9 kpc and a position angle of 43 ± 7◦, very similar to Benjamin et al. (2005). The long bar is much longer than any previously proposed bar/bulge in the Galactic centre, for which the long axis is thought to lie in the range 2.4 (e.g. Dwek et al. 1995) to 3.5 kpc (e.g. Bissantz & Gerhard 2002). Most interestingly, the long bar seems to lie at a dif- ferent angle to the traditional bar/bulge (hereafter ‘Galactic bar’), for which estimates typically lie in the range 25 ± 10 degrees (Vanhollebeke et al. 2009). Most recently, Cabrera-Lavers et al. (2007) have used wider field 2MASS data to isolate red clump stars over a wide range of Galactic longitude, showing that the long bar, the triaxial bulge and the disc can all be seen as separate structures (at least on our side of the Galaxy) with the im- portant result that the long bar is much flatter than the bulge. The long bar is found to have a vertical scale-height of around 100 pc whereas the bulge has a vertical scale- height of 500 pc. Furthermore, the red clump stars used to trace the long bar can be used to estimate its total mass (assuming that the mass to light ratio of red clump stars is the same in the long bar as it is in the bulge) as 6 × 109 M  , of the same order of the mass of the Galactic bar 1 × 1010 M  , (Zhao 1996; Weiner & Sellwood 1999). c© 2010 RAS III 2 Esko Gardner and Chris Flynn Figure 1. Local velocity space, data includes all the stars with U and V velocities from the Geneva-Copenhagen survey (Holmberg et al. 2007) and Schuster et al. (2006). The Hercules stream is the prominent feature area around (−25,−30) km s−1. The velocities have been corrected to the solar velocity from Scho¨nrich et al. (2010) Dehnen (2000) has shown that a bar in the inner galaxy has interesting dynamical consequences for the orbits of stars even in the solar neighbourhood, 8 kpc from the centre, due to a resonance of stellar orbits with the Outer Lindblad Resonance (OLR) of the bar. The manifestation of this res- onance in our local space is the ζ Herculis stream (Eggen 1971b), or the ‘Hercules stream’ as it is usually called. In this paper, we use a similar method to Dehnen (2000), to reexamine the dynamical effects on local disc stars of these non-axisymmetric features in the Galactic centre. We have examined both the long bar (see Lo´pez-Corredoira et al. 2007) and the Galactic bar (see Bissantz & Ger- hard 2002), separately, and as a dual-bar system. We adopt R  = 8 kpc throughout the paper. This facilitates compar- ison with Dehnen (2000), who adopted the same value. We have made a dynamical study of the effects of the bar(s) on local velocity space, by the means of potential theory. We will present the statistical evidence produced by various variations of the long bar-system, and the Galactic bar-system, as well as, some interesting combined models. 2 THE HERCULES STREAM In Fig. 1 we show local velocity space in (U, V ), where U is the velocity towards the Galactic centre and V is the ve- locity in the direction of Galactic rotation. The figure has been prepared using 12939 stars from the Geneva Copen- hagen Survey (Holmberg et al. 2007) of abundances, ages and kinematics of nearby stars, and has been supplemented by 1535 stars from the the high velocity catalogue of Schus- ter et al. (2006). The bin size is 2 km s−1, and the contours are in 5 percent steps relative to the maximum density. The velocities have been corrected for the solar motion taken from Scho¨nrich et al. (2010) of (U  , V  ) = (11.1, 12.2) km s−1. The Hercules stream is the feature at approximately (U, V ) = (−25,−30) km s−1. This is the feature which Dehnen (2000) identifies as having a dynamical origin as a resonance with the bar. Our simulations in this paper of a Galactic bar, long bar (and both) are primarily constrained by this feature. The bulk of the stars in Fig. 1 belong to the compli- cated feature above the Hercules stream. The classically identified features in this region are the Pleiades, Sirius and Hyades streams, which occur at around (0,−10), (20, 20), and (−25,−5) km s−1 respectively. Further streams than those mentioned have been identified by a number of groups, but these mainly concern higher velocity stars and cannot be pointed to easily in Fig 1. One of the very interesting fea- tures at higher velocities is the Arcturus stream at V ≈ 100 km s−1 (for a wide range of U velocities) identified by Eggen (1971a) and discussed recently by Williams et al. (2009), who have identified it in the RAdial Velocity Experiment (RAVE) survey. The Arcturus stream may have a dynami- cal origin, as we show using the simulations in this paper. 3 MODELLING THE DISC AND BAR We model the effects of a bar on local disc stars by setting up a potential describing the disc, bulge, dark halo and bar, and integrating a library of orbits passing through the Solar neighbourhood, similarly to Dehnen (2000). We begin with the axisymmetric Galactic potential of Flynn et al. (1996). This model contains a disc, bulge and dark halo. The scale-length of disc matter (hR) in the model is 4.4 kpc, which these days appears rather long, as recent studies of the disc scale-length give results closer to 3 kpc. Analysis of NIR data in J and K with DEep Near Infrared Survey of the Southern Sky (DENIS) give hR = 2.3 ± 0.1 kpc (Ruphy et al. 1996), while K giants seen in 2MASS give hR = 3.3 ± 0.1 kpc (Lo´pez-Corredoira et al. 2002). These surveys tend to probe the scale-length of luminosity in the disc, rather than the mass, which is traced rather better via main sequence dwarfs. Hubble Space Telescope (HST) star counts of M dwarfs yield a disc scale-length of hR = 3.0±0.4 (Gould et al. 1997). F and K dwarfs, traced in huge num- bers in the 6500 deg2 SDSS survey, yield a disc scale-length of hR = 2.6 ± 0.5 (Juric´ et al. 2008). For a more complete review of the disc scale-length, see Gardner et al. (2009). We chose a disc scale-length of 3 kpc to reflect the com- bined observational constraints. This was straightforward and only involved modifying the Miyamoto discs from which the exponential-like disc is built up. The surface density (Σ  ) of the new disc at the Sun’s position is 50 M  pc−2 and the local density ρ  is 0.10 M  pc−3, consistent with observational constraints (Holmberg & Flynn 2000). The surface density of the disc as a function of Galactocentric radius is shown in Fig. 2. The new disc model has a 3 kpc scale-length, with a truncation at about 18 kpc. The full parameters of the new disc model are shown in Table 1, along with the other parameters of the model (which are unchanged but are reproduced for clarity). The equations for the axisymmetric potential are: c© 2010 RAS, MNRAS 000, 1–9 Probing the Galaxy’s bars via the Hercules stream 3 Figure 2. The surface density of the disc component (solid line) as a function of Galactocentric radius. The dashed line corre- sponds to an exponential density falloff, 3 kpc, which is a good fit to the model over a wide range of radii. Note that the density truncates strongly at 18 kpc. Table 1. Parameters of the model. Parameter Value Unit Vh 220 km s −1 r0 8.5 kpc MC1 3 10 9 M  rC1 2.7 kpc MC2 16 10 9 M  rC2 0.42 kpc Md1 77.04 10 9 M  ad1 5.81 kpc Md2 −68.48 10 9 M  ad2 17.43 kpc Md3 26.75 10 9 M  ad3 34.84 kpc b 0.3 kpc Φ = ΦH +ΦC +ΦD +Φbar (1) ΦH = 1 2 V 2 h ln(r 2 + r20) (2) ΦC = − GMC1√ r2 + r2C1 − GMC2√ r2 + r2C2 (3) ΦDn = −GMDn√ (R2 + (adn + √ (z2 + b2))2) , (4) where ΦD = Σ 3 n=1ΦDn , r = x 2 + y2 + z2, and R = x2 + y2. We have used a different method to modelling the bar as was adopted in Dehnen (2000) and Minchev et al. (2007). Instead of modelling the bar as a local quadrupole perturba- tion to the potential, we model the bar directly in the inner galaxy. This means that orbits of stars entering into the in- ner disc regions and even passing through the bar are much more accurately modelled, even if it is computationally far more expensive. Table 2. Default simulation parameters. Parameter Value Unit Galactic Bar Angle 25 ◦ Dimensions 3.5:1.4:1.0 kpc Mass 10 109 M  Pattern Speed 55.9 km s−1 kpc−1 Local standard of rest 239 km s−1 OLR 1.87 Long Bar Angle 43 ◦ Dimensions 3.9:0.6:0.1 kpc Mass 6 109 M  Pattern Speed 54.9 km s−1 kpc−1 Local standard of rest 235 km s−1 OLR 1.87 The bar is modelled as a Ferrers n = 2 potential, as laid out by Pfenniger (1984), and as such, is a triaxial ellipsoid. The equations for the bar are too complex to present here (for full details see Pfenniger 1984). The bar rotates as a rigid body. A bar is specified by its three axis lengths, its mass, its angular (or pattern) speed and its angle to the line connect- ing the Galactic centre and the Sun. We have two standard bars, with parameters taken from Lo´pez-Corredoira et al. (2007) and from Bissantz & Gerhard (2002). The dimensions of the ‘long bar’ are (3.9:0.6:0.1) kpc, with a mass of 6 109 M  and an angle of 43◦. The ‘Galactic bar’ has dimensions of (3.5:1.4:1.0) kpc, a mass of 1010 M  and an angle of 25◦. These masses correspond to a Dehnen (2000) α value of 0.0037, for the long bar case, and 0.0040, for the Galactic bar case. Note that in the simulations we do not vary the dimensions of the respective bars, only the mass and position angle. An important parameter is the pattern speed of the bar. Following Dehnen (2000) and Minchev et al. (2007), we parametrise the bar’s pattern speed by the OLR, i.e. by the ratio of the period of a star with velocity of the local standard of rest (LSR) in a given Galactic model to the period of rotation of the bar. The definition of the LSR is straightforward in axisym- metric models of the Galaxy, but needs to be redefined slightly once we have added a bar to the model. In the pres- ence of a bar, a star can no longer follow a circular orbit at the solar radius, but rather oscillates around a mean radius. For typical bar parameters, the oscillation is up to 10 km s−1, compared to a mean velocity of 220 km s−1 in typical models. For each model run in our simulations, we deter- mine, experimentally, at what velocity a star, which starts at the Solar position, would have a mean orbital radius of 8 kpc. Integration is carried out for 10 Gyr, so that the effect of the angle of the bar is averaged out. This is defined as the LSR for that model. 4 SIMULATIONS AND METHODS The basic simulation approach is to model the appearance of local velocity space using the method described by Dehnen c© 2010 RAS, MNRAS 000, 1–9 4 Esko Gardner and Chris Flynn (2000), a ‘backward restricted N -body method’. The idea is to trace a library of orbits starting from the solar neighbour- hood backwards over some fixed time (typically 1 Gyr), and use the known density distribution and velocity distribution of the disc to reconstruct the velocity distribution of stars in the solar neighbourhood. We introduce two refinements, firstly, a more realistic representation of the bar potential and secondly, we use ob- servational constraints on the velocity dispersion and asym- metric drift of stars in the Galactic disc over a wide range of Galactocentric radius obtained by Lewis & Freeman (1989). In practise, this latter refinement is not very different from the approach adopted by Dehnen (2000), (simple exponen- tial functions of stellar density and radial velocity dispersion for the old stellar disc) and indeed we find that we repro- duce his results quite well when we adopt a bar similar to his. Our simulations are in two dimensions (2D), as in Dehnen (2000). The initial conditions are set by a velocity-grid, of 100 by 100 elements in the UV -velocity plane, covering U = −100 to 100 km s−1 , and V= −150 km s−1) to V = 50 km s−1), in steps of 2 km s−1. A tracer with the UV -velocity selected from the grid is then set into motion, following the tracer backwards in time over a period of 1 Gyr, tracking its position every 1 Myr. The initial position is selected as the solar position (8 kpc), and both the z-coordinate and the W -velocity are set to zero. The initial velocity in the rotational direction V is, of course, relative to the LSR for particular model being tested. Statistically, we use a density and velocity-based system, where a weight value corresponds to a 8 kpc- centric exponential distribution, and a Gaussian velocity- distribution, with standard deviations as described in Lewis & Freeman (1989). The V -velocity is corrected back to a ‘true’ V -velocity using the LSR of the position and the asymmetric drift (Lewis & Freeman 1989). 5 RESULTS We examine in this section the effects of both the ‘Galactic bar’ and the ‘long bar’ on local velocity space. We vary the bar mass, value of the OLR (i.e. the pattern speed of the bar), and the angle of the bar. Finally, we look at the effects of including both bars at the same time. 5.1 Effect of the mass of the Galactic bar We begin by running the standard model of the ‘Galactic bar’ but varying its mass over a range of 1 to 10 × 109 M  , showing the results in Fig. 3. The OLR for this case is 1.87. The variations of mass within one order of magnitude have prominent effects. As found by Dehnen (2000) we find that a Hercules-like feature arises in local velocity space in the presence of a bar. We will henceforth refer to this as a secondary feature in velocity space, with the primary fea- ture (above it) containing most of the stars. The effect of increasing the mass of the bar is similar to that found by Dehnen (2000), with the strength of the secondary feature increasing (and the depth of the resonance gap between the secondary feature and primary feature also growing). The secondary component becomes more ellipsoidal, as mass is increased, and, as expected, becomes more prominent. Also, the position of the secondary feature moves to more negative V velocities as the mass is increased. While it is clear that something like the Hercules feature can be generated by this bar, the detailed match is not yet perfect. In Fig. 1, the Her- cules feature is found at (U, V ) = (−25,−30) km s−1. In the simulations, the generated feature approaches this location for the higher mass models, but at the expense of generating a sharp ‘tail’ of stars at positive U values spreading outward from the primary feature. This tail of stars is also seen in the Dehnen (2000) and Minchev et al. (2007) simulations, and increases with the mass of the adopted bar. Such a tail, if present at all in the observational data (Fig. 1) is quite weak. Apart from this tail, the generated feature is fairly consistent with various estimates of bar mass where they are consistently around and over 1010 M  (see e.g. Zhao 1996; Weiner & Sellwood 1999). 5.2 Long bar We next look at the ‘long bar’ case. We adopt the standard model for its dimensions and initial position angle and vary only the mass. These models all give a Hercules-like sec- ondary feature, similar to the Galactic bar case. The effect of varying its mass on local velocity space is shown in Fig. 4. One can see that mass takes a significant role in the division of the primary and secondary feature, where the secondary is stronger with increasing mass. Again, these plots superfi- cially resemble the real local velocity space (Fig. 1) but not in detail. The position of the secondary feature tends to be at too high V velocities and the tail of stars at high U values also appears for the higher mass models. One very interesting ‘tertiary’ feature appears with the long bar, located at V = −80 to −100 km s−1. The tertiary feature becomes more prominent as mass increases. This fea- ture is not seen in the Galactic bar simulations (Fig. 3). We will return to this feature in section 5.4. 5.3 Pattern speed of the bar As shown by Dehnen (2000) and Minchev et al. (2007), the position of the Hercules feature is quite dependent on the adopted pattern speed of the bar, parametrised by the OLR. The feature can therefore be used to measure the pattern speed of the bar. Dehnen (2000) finds that the feature can be fit by an OLR of 1.85 ± 0.15, using similar simulations to ours. Minchev et al. (2007) find a similar and more tightly constrained value of the OLR of 1.87 ± 0.04, using similar simulations as here but also from constraints implied by the locally measured Oort constant C. Our standard model for both the Galactic and long bars assumes an OLR of 1.87. For the Galactic bar case, the effect of shifting from an OLR value of 1.87 to 1.90 is shown in Fig 5. The secondary feature shifts by about 5 km s−1 to more negative V velocities. A more detailed view of this is shown for the long bar case in Fig. 6, where we have varied the value of the OLR ranging from 1.83 to 1.95. The secondary feature shifts quite markedly away from the primary feature over this range of assumed OLR values. We derive a best fitting OLR of 1.87±0.02, assuming that the position of the secondary feature is at V = −30± 5 km s−1. c© 2010 RAS, MNRAS 000, 1–9 Probing the Galaxy’s bars via the Hercules stream 5 Figure 3. Effect on local velocity space of varying the mass of the Galactic bar, from 109 to 1010 M  in increments of 109 M  . For other simulation parameters see Table 2. Our determination of the OLR value for the bar is in agreement with Dehnen (2000) and Minchev et al. (2007) within the error bars. We note that our value is about the same as theirs, even though we adopt the solar motion of Scho¨nrich et al. (2010), which shifts this feature by 7 km s−1 to more positive V velocities in the observational plane (the Hercules feature is centred on about V = −35 km s −1 with the solar motion adopted by Dehnen (2000), rather than at Figure 4. Variations of velocity space with different values of for long bar mass. Masses are 4, 5, 6 (middle), 7, and 8 109M  . For other simulation parameters see Table 2. Figure 5. Variation of velocity space with a Galactic bar model, OLR values of 1.87 and 1.90. V = −30 km s −1 as in our Fig. 1. Had we adopted the older solar motion, we would obtain an OLR of 1.90± 0.03. 5.4 Bar angle Dehnen (2000) used the Hercules stream to probe the angle of the bar, deriving an angle of 25◦. This is to be compared to the observational evidence for the (Galactic bar) angle, which lies in the range 20◦ to 40◦. We probe the effect of bar angle as well, for both the Galactic bar and long bar cases. As found by Dehnen (2000), the effect is to change of the position of the secondary feature to towards more negative values of U and more positive values of V with increasing bar angle, and also to change the strength of the feature, as c© 2010 RAS, MNRAS 000, 1–9 6 Esko Gardner and Chris Flynn Figure 6. Effect on local velocity space with different values of for the OLR for the long bar case. Values are 1.83, 1.85, 1.87, 1.90, 1.93, and 1.95, from left to right, top to bottom. There is a clear shift in the position of the secondary feature to lower V velocity as the OLR increases, allowing a determination to be made of the pattern speed of the bar. seen in Fig. 7 and Fig. 8. However, these changes are rather subtle and we consider them too weak to constrain the bar angle. The tertiary feature is clearly dependent on bar angle, both for the long bar case and the Galactic bar case. The feature appears in both cases when the initial bar angle is greater than 40◦. The orbits of stars in this feature are highly eccentric and take them to within 2 kpc of the Galactic cen- tre, ideal for direct interaction with a bar. This accounts for why bar angle is important to the generation of the tertiary feature. There is a stream in the Solar neighbourhood at about this velocity, V ∼ −100 km s−1, the Arcturus stream (Eggen 1971a). The stream has recently been recovered in the RAVE survey by Williams et al. (2009). Williams et al. (2009) ar- gue against an accretion origin for this stream, because its stars have a wide range of metallicity, so that the accreted system would have to have an unreasonably high mass, per- haps as large as the Large Magellanic Cloud (LMC). It is interesting that our simulations show that a dynamical ori- gin for the stream is plausible. Streams with a dynamical origin are expected to have a wide range of metallicity, as wide as the stars of its parent population. If so, the more plausible bar for originating this stream is the long bar, be- cause its measured bar angle is 43◦. The Galactic bar could Figure 7. Variation of Galactic bar angle, 15◦, 25◦, 35◦, and 45◦ Figure 8. Variation of long bar angle, 23◦, 33◦, 43◦, and 53◦. also produce this stream, but it would have to also be at an angle of about 45◦, whereas most studies find that it lies much closer to the line of sight to the Galactic centre, at about 20◦ (Vanhollebeke et al. 2009). 5.5 Integration time We have chosen an integration time of 1 Gyr, but this choice is somewhat arbitrary. Fig. 9 shows the effects of adopting shorter and longer integration times of 500 Myr and 2 Gyr for the long bar. After 500 Myr, the secondary feature has started form- ing but is not yet prominent. This represents about 4 bar rotations, and it is not surprising that it takes some time for the resonance with the bar to develop the trough be- tween the primary and secondary features. After 2 Gyr we c© 2010 RAS, MNRAS 000, 1–9 Probing the Galaxy’s bars via the Hercules stream 7 Figure 9. Variation in integration time, 500 Myr, and 2 Gyr, using the long bar. can see the stronger splitting of the primary and secondary feature, with the resonance-area clearing, as well as the pri- mary feature tailing out to more positive U . Simulations as long as this assume of course that the bar has existed for 2 Gyr with the same parameters, in particular that the rota- tion rate and mass have not changed during this time. Our choice of 1 Gyr is meant to allow structure enough time to develop, while not taxing the credibility of bar properties being very long lived in galaxies like the Milky Way. Fux (2001) has pointed out a potential difficulty with the back- wards integration technique used here. He finds that as one integrates further backwards in time, increasingly spurious structures develop in the UV -plane. We find no such be- haviour over approximately 16 periods of the bar (Fig. 9). We are confident in our 1 Gyr timescale, for interpreting the simulations. 5.6 Dual bar models Our simulations have shown that both the Galactic bar, and the newly discovered ‘long bar’ can generate features in ve- locity space which are superficially similar to those seen in the Solar neighbourhood, in particular the Hercules stream. We now investigate the effect of including both bars. We show a selection of dual bar models in Fig. 10. We begin with both bars having the same OLR - i.e. the two bars have the same pattern speed and rotate together. The up- per panels of Fig. 10 show the local effect on phase space for OLR values of 1.87 and 1.90. The overwhelming impression here is that the main feature in phase space is driven to large values of U , inconsistent with observations. The main cause for this is simply the large amount of mass which now resides in the Galactic centre, in these cases 1.6 1010 M  . Similar behaviour is seen even for cases with a single bar, when the mass of the bar exceeds about 1010 M  . More than this much mass and the bar(s) starts to generate an unaccept- able amount of activity on the Solar neighbourhood (U, V ) distribution, typically pushing stars to too large values of U . Two decoupled scenarios were tested, with one of the bars at an OLR of 1.50 and the other at 1.87. These are shown in the middle panels of Fig. 10. On the right, the Galactic bar has an OLR of 1.50, while the long bar has an OLR of 1.50 on the left. Again, the total mass in the two bars is generating significant structure in local veloc- ity space, with the same problem of too many stars being shifted into a tail of high U velocities. In addition, the ter- tiary feature at V ∼ −80 becomes very prominent, as the chances of an interaction with either bar is now higher. We Figure 10. Dual bar models, phase-locked 1.87 and 1.90 (top), unlocked, long at 1.87, galactic at 1.5 (middle left), long at 1.5 galactic at 1.87 (middle right), half-mass locked at 1.87 (bottom left), half-mass unlocked long at 1.87, galactic at 1.5 (bottom right). tested this by running single-bar models with an OLR of 1.5, which results in no secondary or tertiary feature for the Galactic bar case, and in the long bar case it causes no secondary and a very weak tertiary, such as seen in other single-bar simulations of the long bar. Both the upper and middle panels show more structure than really exists in the Solar neighbourhood, primarily because too much mass is assigned to the bar structures in the inner Galaxy. To alle- viate this problem, we tested a scenario where the masses of both bars are simply halved. This is shown in the lower panels of Fig. 10. On the left, the bars both have an OLR of 1.87, while on the right, the bars are unlocked, with the Galactic bar at an OLR of 1.5 and the long bar with an OLR of 1.87. Reducing the mass brings the simulations into much better agreement with the observations, although the secondary feature is at a too high V velocity. This can be remedied by shifting to a higher value of 1.90 for the OLR. 6 DISCUSSION & CONCLUSIONS We have performed simulations of the effects of bars in the inner Galaxy on the distribution of stars in velocity space near the Sun. Our simulations are similar to those of Dehnen (2000) and Minchev et al. (2007), in which a large library of orbits of stars passing through the Solar neighbourhood are c© 2010 RAS, MNRAS 000, 1–9 8 Esko Gardner and Chris Flynn performed in a model of the Galactic potential including a bar. These earlier studies showed that the Hercules stellar stream in the solar neighbourhood could be the result of a resonance of stellar orbits with a fast moving bar (i.e. with a co-rotation with the Galactic disc not much beyond the end of the bar). Both Dehnen (2000) and Minchev et al. (2007) find that the Hercules stream can be used to constrain the rotation rate and angle of the bar. Our motivation to reexamine this is the recent discovery that the Galaxy may contain a second bar, in addition to the traditional Galactic bar, by Benjamin et al. (2005) and Lo´pez-Corredoira et al. (2007). These authors have found a ‘long bar’ in the Galactic centre. It is quite different than the Galactic bar, being highly flattened (semi-major axes of 3.9, 0.6 and 0.1 kpc as compared to 3.5, 1.4 and 1.0 kpc). Furthermore, it is aligned at an angle of about 43 degrees to the line of sight to the Galactic centre, as opposed to about 20 degrees for the Galactic bar. For the Galactic bar we adopt a mass of 1010 M  (Zhao 1996; Weiner & Sellwood 1999) and for the long bar 6 × 109 M  (Lo´pez-Corredoira et al. 2007). The long bar has a mass quite comparable to the traditional bar. We have simulated the effects of the long bar, as well as the Galactic bar, and both bars, on the velocities of stars in the Solar neighbourhood. Rather than simulating bars as a quadrupole perturbation in the local Galactic potential, we simulate the bars with a Ferrers-potential (which gener- ates a triaxial ellipsoidal density distribution). This is com- putationally more expensive, but allows for more accurate modelling of a bar as a whole, especially for stars passing through, or trapped in, a bar. The simulations are performed in 2-D, although the model itself allows 3-D simulations. We confirm the basic picture that the Hercules stream can be generated by a resonance between local stars and the bar. Both the long bar and the Galactic bar produce Hercules-like features in the Solar neighbourhood. The po- sition of the Hercules stream in the simulations is found to depend on the mass of the bar, the rotation rate of the bar, and rather weakly, on the angle of the bar. The geometry of the bars is not very important either, since both bars generate Hercules-like streams. We measure the rotation rate of the bars from the posi- tion of the Hercules stream, finding that both bars can pro- duce acceptable fits to the (U, V ) velocities of nearby stars. We measure the rotation rate of both bars via the Outer Lindblad Resonance value they produce (i.e. by the ratio of the period of a star with velocity of the Local Standard of Rest in a given Galactic model to the period of rotation of the bar). We measure values of 1.87 ± 0.02 for the OLR of both the Galactic bar and the long bar. This value gives a pattern speed of 55.9 and 54.9 km s−1 kpc−1, respectively, for the Galactic bar and long bar. This puts the correspond- ing co-rotation, for the long bar, as defined by its potential model, at the approximate tip of the bar, 3.9 kpc. For the Galactic bar, our model produces no co-rotation radius. It should be noted however that this value is sensitive to the assumed Solar motion. We assume the Solar motion recently derived by Scho¨nrich et al. (2010). Had we adopted the So- lar motion used by Dehnen (2000), we would have derived a value of 1.90 ± 0.03. Both of these derived values are in reasonable agreement with Dehnen (2000)’s determination of 1.85 ± 0.15 and Minchev et al. (2007)’s determination of 1.87± 0.04. After looking at the effects of each bar alone, we tried putting both bars into the modelling. The effects on local velocity space are quite dramatic, as we have now added quite a lot of mass in an non-axisymmetric component to the Galactic centre. In particular, these simulations produce a striking tail of stars in the Solar neighbourhood at high U velocities, inconsistent with observations. The simple ex- pedient of halving the mass of both bars produces much better fits to the data. On the basis of these experiments then, if there are two bars in the Galactic centre, we won- der whether the mass estimates of each bar include material from the other bar, and are perhaps overestimates. This of- fers a challenge to observational studies of the bar masses. Finally, we have found that an additional feature in lo- cal velocity space can arise under some conditions, especially for the long bar case. A weak stream of stars can form at about V ∼ −100 and U < 0 km s−1. We believe this is due to a direct interaction of such stars with the bar, rather than a resonance with the bar. It is tempting to associate this feature with the known Arcturus stream at V∼ −100 km s−1 (Williams et al. 2009) in the Solar neighbourhood. Thus, a dynamical origin for the stream is possible, whereas people have mainly discussed accretion origins (of a satellite galaxy) to date. We point out that our method has shortcomings. Firstly, the simulations are in 2-D, and while this is almost certainly acceptable for typical disc stars with near circular orbits, full 3-D simulations may be needed for features like that which we associate with the Arcturus stream. In ad- dition, there are no spiral arms or other non-axisymmetric components, besides the bars, in the modelling, and these may be responsible for additional complexity in local veloc- ity space which we cannot model. We are not affected by the problems of long-term backwards integration as mentioned in Fux (2001). Compared to earlier studies, we try to more accurately represent the Galaxy’s observed properties, such as local densities and bar geometry. The model also lets us look at direct interaction, as well as resonances, this is es- pecially important for high velocity stars (e.g. those in the Arcturus stream), where they have high enough eccentrici- ties to enter the inner parts of the Galaxy. The model also makes it possible, in the future, to fully study the 3D impact of the bar(s). ACKNOWLEDGEMENTS The authors wish to thank Professor Kimmo Innanen for starting out this project, as well as Professor Daniel Pfen- niger for providing his code for the bar potential. We also wish to thank Burkhard Fuchs and the referee for insight- ful comments. EG acknowledges the financial support of the Finnish Cultural Fund and the Finnish Graduate School in Astronomy and Space Physics. CF acknowledges financial support by the Academy of Finland. REFERENCES Benjamin R. A., Churchwell E., et al. 2005, ApJL, 630, c© 2010 RAS, MNRAS 000, 1–9 Probing the Galaxy’s bars via the Hercules stream 9 L149 Binney J., Gerhard O., Spergel D., 1997, MNRAS, 288, 365 Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Galactic Dynamics: Second Edition, by James Binney and Scott Tremaine. ISBN 978-0-691-13026-2 (HB). Published by Princeton University Press, Prince- ton, NJ USA, 2008. Bissantz N., Gerhard O., 2002, MNRAS, 330, 591 Cabrera-Lavers A., Hammersley P. L., Gonza´lez-Ferna´ndez C., Lo´pez-Corredoira M., Garzo´n F., Mahoney T. J., 2007, A&A, 465, 825 Carlson B. C., 1994, ArXiv Mathematics e-prints Chandrasekhar S., 1969, Ellipsoidal figures of equilibrium. The Silliman Foundation Lectures, New Haven: Yale Uni- versity Press, 1969 Dehnen W., 2000, AJ, 119, 800 Dwek E., Arendt R. G., Hauser M. G., Kelsall T., Lisse C. M., Moseley S. H., Silverberg R. F., Sodroski T. J., Weiland J. L., 1995, ApJ, 445, 716 Eggen O. J., 1971a, PASP, 83, 271 Eggen O. J., 1971b, PASP, 83, 251 Flynn C., Sommer-Larsen J., Christensen P. R., 1996, MN- RAS, 281, 1027 Fux R., 2001, A&A, 373, 511 Gardner E., Innanen K. A., Flynn C., 2009, in J. Andersen, J. Bland-Hawthorn, & B. Nordstro¨m ed., IAU Symposium Vol. 254 of IAU Symposium, Dynamical effects of the long bar in the Milky Way. p. 21P Gould A., Bahcall J. N., Flynn C., 1997, ApJ, 482, 913 Holmberg J., Flynn C., 2000, MNRAS, 313, 209 Holmberg J., Nordstro¨m B., Andersen J., 2007, A&A, 475, 519 Juric´ M., Ivezic´ Zˇ., Brooks A., Lupton R. H., Schlegel D., et al. 2008, ApJ, 673, 864 Lewis J. R., Freeman K. C., 1989, AJ, 97, 139 Lo´pez-Corredoira M., Cabrera-Lavers A., Garzo´n F., Ham- mersley P. L., 2002, A&A, 394, 883 Lo´pez-Corredoira M., Cabrera-Lavers A., Gerhard O. E., 2005, A&A, 439, 107 Lo´pez-Corredoira M., Cabrera-Lavers A., Mahoney T. J., Hammersley P. L., Garzo´n F., Gonza´lez-Ferna´ndez C., 2007, AJ, 133, 154 Minchev I., Nordhaus J., Quillen A. C., 2007, ApJL, 664, L31 Pfenniger D., 1984, A&A, 134, 373 Ruphy S., Robin A. C., Epchtein N., Copet E., Bertin E., Fouque P., Guglielmo F., 1996, A&A, 313, L21 Scho¨nrich R., Binney J., Dehnen W., 2010, MNRAS, 403, 1829 Schuster W. J., Moitinho A., Ma´rquez A., Parrao L., Co- varrubias E., 2006, A&A, 445, 939 Vanhollebeke E., Groenewegen M. A. T., Girardi L., 2009, A&A, 498, 95 Weiner B. J., Sellwood J. A., 1999, ApJ, 524, 112 Williams M. E. K., Freeman K. C., Helmi A., the RAVE collaboration 2009, in J. Andersen, J. Bland-Hawthorn, & B. Nordstro¨m ed., IAU Symposium Vol. 254 of IAU Symposium, The Arcturus Moving Group: Its Place in the Galaxy. pp 139–144 Zhao H., 1996, MNRAS, 278, 488 c© 2010 RAS, MNRAS 000, 1–9 Mon. Not. R. Astron. Soc. 000, 1–8 (2010) Printed 24 June 2010 (MN LATEX style file v2.2) Tidal influence of a variable Galactic potential on the orbits of Oort cloud comets E. Gardner1, P. Nurmi1, C. Flynn2 and S. Mikkola1 1Tuorla Observatory, Department of Physics and Astronomy, University of Turku, Va¨isa¨la¨ntie 20, FI-21500, Piikkio¨, Finland 2Finnish Centre for Astronomy with ESO (FINCA), University of Turku, Va¨isa¨la¨ntie 20, FI-21500, Piikkio¨, Finland Accepted –. Received –; in original form – ABSTRACT The long-term dynamics of Oort cloud comets are studied under the influence of both the radial and the vertical components of the Galactic tidal field. Sporadic dynamical perturbation processes are ignored, such as passing stars, since we aim to study the influence of just the axisymmetric Galactic tidal field on the cometary motion and how it changes in time. We use a model of the Galaxy with a disc, bulge and dark halo, and a local disc density, and disc scale length constrained to fit the best available observational constraints. By integrating a few million of cometary orbits over 1 Gyr we calculate the time variable flux of Oort cloud comets that enter the inner Solar System, for the cases of a constant Galactic tidal field, and a realistically varying tidal field which is a function of the Sun’s orbit. We find that the periodicity in the cometary flux is complicated and quasi-periodic. The amplitude of the variations in the flux are of order 30 %. The radial motion of the Sun is the chief cause of this behaviour, and should be taken into account when the Galactic influence on the Oort cloud comets is studied. Key words: methods: numerical – celestial mechanics – comets: general – Solar system : general – stars: kinematics and dynamics. 1 INTRODUCTION Strong observational evidence supports the idea that the in- ner Solar system is subject to a steady flux of ‘new’ comets which originate from the ‘Oort cloud’ (Oort 1950). The semi- major axes of these comets are thought to evolve under the influence of external forces such as the Galactic tidal field and passing stars. The comets evolve from a typical aorig  3 × 10 4 AU to either hyperbolic orbits or to larger binding energies, depending on the orbital evolution dur- ing the cometary encounters with the Solar System planets. During one orbital evolution, comets typically experience perturbations that change the comet’s orbit into a short- period or hyperbolic orbit, which leads to a subsequent ejec- tion into interstellar space (Nurmi 2001). The Galactic potential exerts a force on Oort cloud comets, and is important for the steady state flux of comets with a > 20000 AU, e.g. Byl (1983); Heisler & Tremaine (1986); Matese & Whitman (1989). The Galactic tidal force has a dominant component that is perpendicular to the Galactic plane; the radial component is 10 times weaker than the perpendicular component Heisler & Tremaine (1986). For this reason, many studies (e.g. Matese et al. 1995; Wick-  E-mail: esgard@utu.fi ramasinghe & Napier 2008) have assumed that the radial Galactic tidal field component is negligible. Other studies have included the radial component, for more accurate mod- elling of cometary motion, such as in Matese & Whitmire (1996); Brasser (2001). The radial component of the tide has been found to have an effect on the long-term evolution of the comets’ perihelia, on the distribution of the longitudes of the perihelion (Matese & Whitmire 1996), and the origin of chaos in the cometary motion (Breiter et al. 2008). In recent large scale simulations, Rickman et al. (2008) found that a fundamental role is played by perturbations due to pass- ing stars on comets, contrary to the investigations during the previous two decades, starting with Heisler & Tremaine (1986). The stellar perturbations, do of course act together with the Galactic tide. The topic of this paper is the effect on comets due to the Galactic tide alone. We use a realistic model of the lo- cal Galaxy, which is well constrained by observations, which contains a disc, bulge and halo. We follow the orbit of the Sun in this model using recent constraints on the Solar mo- tion. This motion allows us to compute the change in the ver- tical and radial components of the Galactic tide with time, and the change in the tidal force due to the radial motion of the Sun is fully accounted for. Qualitatively, the tidal effect on the cometary orbits can c© 2010 RAS IV 2 E. Gardner, P. Nurmi, C. Flynn and S. Mikkola be evaluated by studying the change in angular momentum averaged over one orbit. The Galactic tidal force periodically changes the angular momentum (J = √ GM  a(1− e2)) of the Oort cloud comets (Ferna´ndez & Ip 1991). The angular momentum of the comets changes as (Heisler & Tremaine 1986): dJ dt = − 5πGρ0 G2M2  E 2(E2 − J2)[1− (J2z /J 2)] sin 2ωg . (1) Here, ωg is the Galactic argument of perihelion, Jz is the z-component of angular momentum, and is perpendic- ular to the Galactic plane, E is the total energy, ρ0 is the local mass density, and G is the gravitational constant. The periodically changing angular momentum causes variations mainly in perihelion distance q, for comets in near-parabolic orbits, since q ≈ J2/(2GM  ). The typical assumptions in the cometary flux calculations due to the Galactic tidal force suppose that the local tidal field is axisymmetric, perpendic- ular to the mid-plane, and adiabatically changing (Matese & Whitman 1992). The aim of many of these studies is to correlate the mo- tion of the Sun in the Galaxy with phenomena on the Earth, such as mass extinctions of species, the cratering record and climate change. An extensive review of this topic has been made recently by Bailer-Jones (2009). For example, early studies have shown that the ages of well dated impact craters on Earth are not distributed randomly, but that there is a possible 28 Myr (Alvarez & Muller 1984) or 30 ± 1 Myr periodicity in crater ages over the past 250 Myr (Rampino & Stothers 1984). Since then, several authors have claimed that there is a significant pe- riodic signal present, but the periods differ quite a lot from study to study. The signal is the most prominent for 40 large, well-dated craters that are up to 250 Myr old (Napier 2006), but the period is difficult to measure, with estimates of between 24-26 Myr (Napier 2006), 30 Myr (Napier 2006; Stothers 2006), 36 Myr (Napier 2006; Stothers 2006), 38 Myr (Yabushita 2004; Wickramasinghe & Napier 2008) and 42 Myr (Napier 2006). In some studies, the reliability of the signal is questioned altogether, based on the inaccuracy of the age estimates of the impact craters, possible biases caused by rounding the ages of craters, and the small number of craters (Grieve & Pesonen 1996; Jetsu & Pelt 2000). By critically reviewing many studies that have tried to connect the Solar motion and periodicity in terrestrial phe- nomena such as biodiversity, impact cratering and climate change, Bailer-Jones (2009) has concluded that there is little evidence to support these connections. By studying the arti- ficial cratering data Lyytinen et al. (2009) came to the same conclusion, that the reliable detection of any periodicity is currently impossible with the existing cratering data. In this study, we statistically analyse how the Galactic tidal force changes cometary orbits over 1 Gyr, using nu- merical simulations. The 1 Gyr time-scale is long enough to observe changes in the Galactic tide due to both the ra- dial and vertical motion of the Sun. A simple axisymmetric Galactic potential is adopted. To our knowledge, there has been no study to date, in which the effect of both radial and vertical components, in a time varying Galactic poten- tial (via the variation in mass density ρ), has been analysed in detail. Our purpose is to study the statistical effects of Table 1. Modified and additional potential parameters. The pa- rameters not mentioned here stay the same as in Gardner & Flynn (2010). Property value Unit r0 10 kpc b 0.45 kpc MD1 66.06×10 9 M  MD2 −59.05 ×10 9 M  MD3 22.97×10 9 M  bg 0.12 kpc Mg1 18.63×10 9 M  Mg2 −16.66×10 9 M  Mg3 6.48×10 9 M  the complete Galactic potential to the Oort cloud comets in detail, especially concentrating on the comets that enter the Solar System (q < 30 AU). In particular, we analyse the differences in cometary motion for when the tidal field is constant, and when it varies as the Sun moves in a realistic orbit in a fairly realistic Galactic potential. We find that it is important to include the radial motion of the Sun in the calculations, since the local density varies significantly as the Sun moves towards and away from the Galactic centre. 2 METHODS The method of simulation requires two, traditionally sepa- rate, components. The first is to simulate the motion of the Sun around the Galaxy. The second involves the evolution of the orbits of comets in the Oort cloud. We integrate the mo- tion of the Sun in an axisymmetric Galactic potential. The method of integration of the comets is described in Mikkola & Nurmi (2006). We do not include the random perturba- tions caused by the planets, since the aim is to identify the tidal effects of the Galactic potential on the flux of comets reaching the inner Solar System. 2.1 The Galactic potential The Galactic potential consists of a disc, bulge, and dark halo, and is described in Gardner & Flynn (2010). With a notable exception in the disc-model, as well as a slight modification of the dark halo. We noticed that the vertical density profile of the disc, from Gardner & Flynn (2010), does not accurately reproduce the observational profile from Holmberg & Flynn (2004). We proceeded to modify our disc- model by changing a few of the model’s parameters, and adding three more Miyamoto-Nagai potentials to the model, to emulate a very thin layer of gas in the disc. The equation for the new disc component is: Φg = 3∑ n=1 −GMgn√ R2 + [adn + √ (z2 + b2g)]2 , (2) where G is the gravitational constant, R is the Galacto- centric radius, and z is the height. The modified parameters are in Table 1. The vertical density profile of the model can be seen in Fig. 2. c© 2010 RAS, MNRAS 000, 1–8 Dynamics of Oort cloud comets 3 Figure 1. The surface density of the disc component as a func- tion of Galactocentric radius. The dashed line corresponds to an exponential density falloff of 3 kpc, which is a good fit to the model over a wide range of radii. Note that the density truncates strongly at 18 kpc. The mass density at the Sun is 0.11 M  /pc3, consistent with observational constraints (0.10 ± 0.01 M  /pc3 Holm- berg & Flynn (2000) and 0.105 ± 0.005 M  /pc3 Korchagin et al. (2003)). This corresponds to a nominal Tz=83 × 10 6 year period for small amplitude simple harmonic motion in the vertical direction. The surface density of disc mat- ter in the model is 54.9 M  /pc2, compared with a mea- sured disc surface density of 56 ± 6M  /pc2 Holmberg & Flynn (2004). The adopted current position of the Sun in the model, (R, z)  , is (8,0) kpc. The local circular velocity of the model is 221 km s−1. 2.1.1 Density and rotation curve Fig. 1 shows the surface density of the disc with radius (R), the change in density with height (z) at the Sun (R = 8kpc) (Fig. 2), and the rotation curve (Fig. 3), as these are the two factors that have the most impact on the orbits of the comets. The disc has a scale-length of 3 kpc, and a local scale-height of 0.24 kpc, consistent with recent measure- ments by Juric´ et al. (2008) using the Sloan Digital Sky Survey (SDSS). The potential is axisymmetric, so we do not model ef- fects such as molecular clouds, spiral arms, and bubbles in the interstellar matter and passing stars. Our interest here is on the global Galactic effects of the Solar motion on the comets. 2.1.2 Motion of the Sun The orbit of the Sun, using the Solar motion measured by Scho¨nrich et al. (2010), where (U, V,W )=(11.1, 12.24, 7.25) Figure 2. Vertical density of the model, at the Sun (R = 8 kpc). The dotted line represents the baryonic contribution in the model, the dashed line the dark matter contribution, and the dashed- dotted line the total density of the model at a certain height (z). The solid line corresponds to an exponential fit to the baryonic component of the model of 0.25 kpc. 0 50 100 150 200 250 300 0 10 20 30 40 50 V c ( k m /s ) R (kpc) Dark Halo Bulge Disc Total Figure 3. Rotation curve for the model Milky Way (dotted), and the different contributions of the disc (dashed), bulge (long dashed), and dark halo (solid). km s−1, and assuming (R, z)  to be (8,0) kpc, is shown in the (R, z)-plane in Fig. 4. Properties of this orbit, such as eccentricity, radial and vertical period and maximum z height, are shown in Table 2. The eccentricity of the Solar orbit (e) is defined as: (Rmax −Rmin)/(Rmax +Rmin). Table 2. Eccentricity, e, maximum vertical height, zmax, radial oscillation period, TR, and vertical oscillation period, Tz , of the Sun in the adopted potential. The adopted solar motion is that of Scho¨nrich et al. (2010). Property value Unit e 0.059 ± 0.003 zmax 0.102 ± 0.006 kpc TR 149 ± 1 Myr Tz 85 ± 4 Myr c© 2010 RAS, MNRAS 000, 1–8 4 E. Gardner, P. Nurmi, C. Flynn and S. Mikkola -0.15 -0.1 -0.05 0 0.05 0.1 0.15 7.8 8 8.2 8.4 8.6 8.8 9 z (k p c) R (kpc) Figure 4. Motion of the Sun in the model potential for 1 Gyr, in radial R and vertical z components. The initial position of the Sun is (R, z)0 = (8, 0) kpc. 2.2 Variation of the tidal parameters caused by the motion of the Sun Our integration method for cometary motion is based on the computation of secular motion of the comet in quadratic perturbation by the Galactic potential. The method used is presented in Mikkola & Nurmi (2006). As the comet orbits the Sun under the Galactic potential it experiences a force F per unit mass: F = − GM  r3 r −G1xx−G2yy −G3zz, (3) where G1 = −(A − B)(3A + B), G2 = (A − B) 2, and G3 = 4πGρ(R, z)−2(B 2 −A2). Here r is the Sun-Comet dis- tance, ρ(R, z) is the local mass density, and G is the gravita- tional constant (Heisler & Tremaine 1986). A and B are the Oort constants, and are obtained from the Galactic model. Usually the radial components (x,y) are neglected, so that G1 = G2 = 0. We will examine the effect on comets in two particular cases. Firstly we study what we call the ’dynamic’ case, in which the Galactic tide changes realistically along the Solar orbit, for the case that the Sun oscillates both vertically and radially in the potential. Secondly, we assume that there is a ’constant’ tidal field with time (i.e. the tidal field does not change as the Sun moves around the Galaxy). In the constant case the Sun moves on a flat, circular orbit. We integrate the Solar orbit in the ’dynamic’ case for 1 Gyr, sampling the values of the G-parameters every 100 kyr: these are shown in Fig. 5. The changes in G3 are dominated by the changes in local density during the orbit, since the changes in A and B over the orbit are not very large. Fig. 5 (top and centre panel) shows how the changing values of the Oort constants, A and B, affect the values of G1 and G2. Fig. 6 shows the combined effects of the radial (R) and vertical motion (z) of the orbit, on G3. It is clear that G3 increases in two situations: when the radial position is the closest to the Galactic centre, and when the vertical motion crosses the mid-plane. Due to the slightly eccentric motion of the orbit, it is the radial component of the motion which dominates the changes in local density (ρ(R, z)), rather than the vertical motion. This is partly because we adopt the Solar motions of Scho¨nrich et al. (2010), which produces a -7.8 x 10-16 -7.6 x 10-16 -7.4 x 10-16 -7.2 x 10-16 -7.0 x 10-16 -6.8 x 10-16 -6.6 x 10-16 -6.4 x 10-16 0 200 400 600 800 1000 G 1 ( 1 /y r2 ) Time (Myr) 6.4 x 10-16 6.6 x 10-16 6.8 x 10-16 7.0 x 10-16 7.2 x 10-16 7.4 x 10-16 7.6 x 10-16 7.8 x 10-16 8.0 x 10-16 8.2 x 10-16 0 20 40 60 80 100 G 2 ( 1 /y r2 ) Time (Myr) 3.0 x 10-15 3.5 x 10-15 4.0 x 10-15 4.5 x 10-15 5.0 x 10-15 5.5 x 10-15 6.0 x 10-15 6.5 x 10-15 0 20 40 60 80 100 G 3 ( 1 /y r2 ) Time (Myr) Figure 5. Variation of the parameters G1 (top panel), G2 (centre panel), and G3 (bottom panel) along the Solar orbit in the ’dy- namic’ case (i.e. including full vertical and radial motions) over 1 Gyr. mildly eccentric orbit for the Sun (e = 0.059 ± 0.003): it oscillates between Galactocentric radii of 7.9 and 8.9 kpc. As such, the evolution of G3 depends on the radial motion and the vertical motion, both being of equal magnitude (Fig. 6). In the second case studied, the Sun is set on a perfectly circular and flat orbit, so that the local mass density does not change with time. For the ’constant’ case, the values of Gi are the mean values from the ’dynamical’ case, and have c© 2010 RAS, MNRAS 000, 1–8 Dynamics of Oort cloud comets 5 7.8 8 8.2 8.4 8.6 8.8 9 0 1 2 3 4 5 6 7 8 9 10 R ( kp c) -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0 1 2 3 4 5 6 7 8 9 10 Z ( kp c) 3.2 3.6 4 4.4 4.8 5.2 5.6 6 6.4 0 2 4 6 8 10 G 3 ( 1 e -1 5 / yr 2 ) Time (Myr) Figure 6. Galactocentric radius R, vertical height z and G3 for the ’dynamic’ case of the Solar motion over 1 Gyr in the Galactic adopted model. Upper panel: the radial motion. Middle panel: vertical motion. Lower panel: G3. the values G1 = −7.00897 × 10 −16, G2 = 7.27613 × 10 −16, and G3 = 4.55749× 10 −15yr−2. 3 SIMULATIONS OF THE COMETS Observations support the idea that the in situ flux of new comets is roughly linear with respect to heliocentric dis- tance, up to the distance of Jupiter (Hughes 2001). Comets with highly eccentric orbits with randomised directions of motion have a q-distribution that is flat (O¨pik 1966). If a long period comet has a perihelion distance between ∼ 10−15 AU, it is quickly ejected to interstellar space or per- turbed so that it becomes a short-period comet (Wiegert & Tremaine 1999). From the perspective of the Oort cloud, the comets are removed from the Oort cloud and from the loss cone to the orbital parameter distribution (Ferna´ndez 1981). The loss cone (lc) is the population of comets that have orbits that will allow them to penetrate the planetary sys- tem, making it possible to observe them (Hills 1981). Plan- etary perturbations move comets efficiently from q values less than qlc  15 AU into either hyperbolic orbits or into orbits that are more tightly bound to the solar system (Hills 1981). In the steady state situation, the new comets are dis- tributed uniformly to the perihelion distances of q  qlc, for a > 30000 AU, while q > qlc comets come also from the in- ner region. For this reason, in all of our simulations we have assumed that initial perihelion distances are distributed uni- formly outside the loss cone. The inclination distribution of the outer Oort cloud of comets is isotropic (uniform in cos i). All the other angular elements are uniformly distributed be- tween [0− 2π]. 3.1 Structure of the Oort cloud An important issue in studying the structure of the Oort cloud, is what energy distribution to adopt for the comets. We assume that the density of comets between 3000 AU and 50000 AU in the Oort cloud is proportional to 1/r3.5±0.5, so that the number of comets N is dN ∝ 1/rαdr, where α = 1.5± 0.5 (Duncan et al. 1987). The conclusions of this paper are not particularly sensitive to the adopted value of α. The existence of an inner Oort cloud has been speculated upon in many studies although there is no direct evidence for it (Hills 1981). An inner Oort cloud is the extension of the Kuiper belt, filling the gap between the Kuiper belt and the outer Oort cloud. Semi-major axes in the inner Oort cloud are typically between 50−15000 AU (Leto et al. 2008). The steady state flux from the inner Oort cloud cannot be uni- formly distributed in perihelion distance, since the planetary perturbations move comets efficiently from q values less than qlc  15 AU to either hyperbolic orbits or into orbits that are more tightly bound to the solar system (Hills 1981). In the steady state situation, the new comets come uniformly to perihelion distances q  qlc when a > 30000 AU, while the q > qlc comets also come from the inner region. 3.2 Simulation parameters of the comets Due to this complicated picture, we study different semi- major axes in separate simulations, and evaluate the effi- ciency of tidal injection in each simulation. We chose initial sample conditions for the Oort cloud comets, setting the semi-major axis to be 10000, 20000, 30000, 40000, 50000, and 60000 AU. We choose the comet’s eccentricity (e) ran- domly, so that the resulting values of q would be uniformly distributed from 35 to the value of the semi-major axis (a). All the other parameters, cos(i), Ω, ω, and the initial ec- centric anomaly were all chosen randomly in appropriate intervals. The number of comets in each simulation is 106 at each of the sampled semi-major axes. For the computation of the numbers of comets reaching the inner Solar System, the results from each semi-major axis are normalised to the adopted number density law. 3.3 Simulations of the Galactic tide To study the effect of the Galactic potential on cometary motion, we chose two hypothetical solar systems: a ‘con- stant’ background density, where the Sun would be on a pure circular orbit with no vertical motion, and a realistic ’dynamic’ Solar orbit. In all systems, we analyse the flux of Oort cloud comets into the Solar System. This means that we consider a comet to have been detected in the inner Solar System when its q is within 30 AU, and it has a heliocentric radius of less than 1000 AU. The last criterion is important, as the osculating elements of comets can evolve to have q  30 AU, far away from the Sun, and evolve to more than 30 AU, without ever entering the inner Solar System. Comets which have been detected are removed from the simulation. We do not replace them with new comets. 4 RESULTS Our aim is to determine the effect of the Galactic tide on Oort cloud comets by separately adopting a constant local density and a varying local density. Examining the flux of comets into the inner Solar System (Table 3), we find that there is no significant difference between the two models. Most of the detected comets come from the middle (30−40 kAU) range of semi-major axes. Likewise, the resulting dis- tribution of orbital elements, for comets coming into the c© 2010 RAS, MNRAS 000, 1–8 6 E. Gardner, P. Nurmi, C. Flynn and S. Mikkola Table 3. Relative fraction of comets entering the inner Solar System, from each of the simulated semi-major axis values, for the ’constant’ and ’dynamic’ cases. Semi-major axis (AU) constant dynamic 10000 0.0729 0.0732 20000 0.1803 0.1804 30000 0.2405 0.2405 40000 0.2381 0.2369 50000 0.1504 0.1507 60000 0.1179 0.1184 Solar System, is the same in both the constant and dynamic cases. Fig. 7 shows example distributions for the simulation cases where the semi-major axis is 30000 AU. Similarly, the q-distribution of incoming comets shows no significant dif- ference between the two cases, as seen in Fig. 8. Finally, the total flux of comets entering the inner-Solar system is not significantly different between the two cases. For the ’constant’ case, from a source pool of 1.8 million comets in the simulated Oort cloud, we find that approximately 120 comets per Myr reach the inner Solar System (which here means comets with q < 30 AU). Assuming that there are about 1012 comets in the Oort cloud (Wiegert & Tremaine 1999), this corresponds to about 70 comets per year with q < 30 AU. This is a bit higher than the observed num- ber of new comets with q < 30 AU, which is about 20 per year, assuming a cometary flux of 0.65 ± 0.18 yr−1 AU−1 (Ferna´ndez 2010). 4.1 Time evolution of the cometary flux There is a clear difference between the two models when we look at the temporal evolution of the cometary flux. The top panel of Fig. 9 shows, for the ’dynamic’ case (and after a relaxation time of about 100 Myr) that the mean cometary flux (shown as a 10 Myr moving average) is well correlated with the changes in G3 (dashed line). The resulting fluxes from the simulation have been weighted according to their relative number-density (a−1.5), the resulting total amount of comets in the simulation, by using this weighting system is 1.8 ×106. From Section 2.2, there are two main causes for the evolution of the flux, the major being the radial motion, the minor being the vertical motion. The top panel of Fig 9 shows the major trend clearly following the radial motion. The vertical motion is also followed, as is clear in the 10 Myr moving average (which smooths out some of the Pois- son noise in the individual 1 Myr samples). We also ran a separate simulation with a circular orbit, and with the verti- cal component intact, the resulting fluxes followed perfectly the vertical motion, as was expected due to the changes in G3 being solely contributed by the vertical component. In the case where the values for Gi have been kept constant (i.e. corresponding to a completely circular orbit with no vertical motion), there is no evidence of any evolution, as seen in the bottom panel of Fig. 9. Simple Fourier-analysis of the dynamic flux (in the top panel of Fig. 9) finds two distinct periods in the cometary flux. The strongest signal is produced by a pe- riod of 143−167 Myr, and an equal signal from a period of 9 27 45 63 81 99 117 135 153 171 i (degrees) 0 2000 4000 6000 8000 10000 C o u n t 18 54 90 126 162 198 234 270 306 342 Ω (degrees) 0 1000 2000 3000 4000 5000 6000 7000 8000 C o u n t 18 54 90 126 162 198 234 270 306 342 ω (degrees) 0 5000 10000 15000 20000 C o u n t Figure 7. Distribution of the orbital elements i, Ω, and ω for a = 30000. The solid column represents constant case, and the shaded column the dynamic case. There is no significant difference between the two distributions. 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 q (AU) 0 2000 4000 6000 8000 10000 12000 14000 16000 C o u n t Figure 8. Distribution of the perihelion distance q for the in- coming comets, for a = 30000. The solid column represents the constant case, and the shaded column the dynamic case. There is no significant difference between the two distributions. c© 2010 RAS, MNRAS 000, 1–8 Dynamics of Oort cloud comets 7 0 200 400 600 800 1000 Time (Myr) 0×10 0 2×10 −5 4×10 −5 6×10 −5 8×10 −5 1×10 −4 R e l a t i v e F l u x 0 200 400 600 800 1000 Time (Myr) 0×10 0 2×10 −5 4×10 −5 6×10 −5 8×10 −5 1×10 −4 R e l a t i v e F l u x Figure 9. The time evolution of the relative flux of comets in the dynamic case (top panel), and the constant case (bottom panel). The solid line represents the relative flux of comets in 1 Myr bins, the white line a 10 Myr moving average. The dashed line in the top panel represents the evolution of the G3-parameter. 41−45 Myr. The former corresponds to the radial period (152 Myr) and the latter to the half-period of the vertical period (43 Myr). The vertical signal is found to lie in the range 41−45 Myr, and is quasi-periodic, as has been seen earlier by Matese et al. (2001). 4.2 Comparison with earlier results Matese et al. (1995) derived various periods for the verti- cal Solar oscillation, depending on the adopted model of the disc. For example, their ‘No-Dark-Disk Model’ has a crossing-period very close to ours (43 Myr), while their ‘Best-fit Model’ produces a much lower period of 33 Myr. They assumed a W -velocity of 7.5 km s−1, compared to ours of 7.25 km s−1, and does not cause much qualitative difference in the vertical period of the orbit of the Sun. The lower period comes from assuming a considerably higher mid-plane density (ρ ≈ 0.13 M  /pc3) than ours. Matese et al. (2001) found that the radial period of Solar motion modulates the vertical period. However, they used very high local densities of matter in the disc, so that the vertical oscillations dominated the radial oscillations. This meant that the flux of comets into the inner Solar System with each passage through the disc was greatly amplified compared to our simulations. We find that the cometary flux varies with an amplitude of about 20%, whereas Matese et al. (2001) find flux variations of about a factor of two. The much smaller amplitude in the signal which we advocate, even taking into account the radial oscillations, would make finding a period in the scant cratering record very difficult. Fouchard et al. (2006) calculated that a local mass den- sity of ρ = 0.1 M  /pc3 corresponds to a cometary flux of around 104 for their first 500 Myr interval for a source pop- ulation of 106 comets, where observed comets have q  15 AU. Assuming that there are 1012 comets in the Oort cloud. This produces a flux of ∼20 comets/yr, which is a factor of two less than our estimate. The Galactic tide case in Rickman et al. (2008) gives a flux of 100 comets per 50 Myr, from a source popula- tion of 106 comets, where observed comets have evolved from q > 15AU to q < 5 AU. This corresponds to a flux of 2 comets/yr, again assuming an Oort cloud with 1012 comets. Using similar analysis we find that we get a flux of 290 comets per 50 Myr, corresponding to a factor of three larger than Rickman et al. (2008). A review by Bailer-Jones (2009) correlates terrestrial events with cometary signals. One of the more interesting ones in the review is the 140 ± 15 Myr period found by Ro- hde & Muller (2005) in the number of known marine animal genera as a function of time. While this could correspond to the Sun’s radial period, Bailer-Jones (2009) considers that it has not been significantly detected, the main problem being that the entire time-span of the data covers no more than three oscillations. Many other periods, proxies, and studies are mentioned in Bailer-Jones (2009), although the conclu- sion is that there is no proven impact on biodiversity as a result of the orbital motion of the Sun. 5 DISCUSSION AND CONCLUSIONS We have studied the long-term dynamics of Oort cloud comets under the influence of both the radial and the ver- tical components of the Galactic tidal field. Other perturb- ing forces on the comets, such as passing stars or passage through spiral arms are ignored, since we aim to study the influence of just the axisymmetric Galactic tidal field on the cometary motion. We use an axisymmetric model of the Galaxy, and a re- cently revised value for the Solar motion by Scho¨nrich et al. (2010). This leads to vertical oscillations of the Sun with an amplitude of about 100 pc, and radial oscillations over about 1 kpc. The changing tidal forces on the Oort cloud are computed as the Sun orbits for 1 Gyr in this potential, and the flux of comets entering the inner Solar System is computed in simulations. As expected, the cometary flux is strongly coupled to the G3-parameter in the tidal forces, which is dominated by the local mass density seen along the Solar orbit. Both the radial and vertical motions of the Sun can be seen in the cometary flux, although the amplitude of the variations is small, implying that detecting such a signal from the small number of age-dated craters would be very difficult. This agrees with the recent review of the detectability of the Solar motion in terrestrial proxies (Bailer-Jones 2009). As G3 is directly coupled to local density, it is easily af- fected by the non-axisymmetric components in the motion of the Sun in the Galaxy. This implies that spiral arms, a Galactic bar, giant molecular clouds, or any other intermit- tently encountered structure should have an effect on the flux. c© 2010 RAS, MNRAS 000, 1–8 8 E. Gardner, P. Nurmi, C. Flynn and S. Mikkola ACKNOWLEDGEMENTS PN wants to thank the Academy of Finland for the financial support in this work. EG acknowledges the support of the Finnish Graduate School in Astronomy and Space Physics. REFERENCES Alvarez W., Muller R. A., 1984, Nature, 308, 718 Bailer-Jones C. A. L., 2009, International Journal of As- trobiology, 8, 213 Brasser R., 2001, MNRAS, 324, 1109 Breiter S., Fouchard M., Ratajczak R., 2008, MNRAS, 383, 200 Byl J., 1983, Moon and Planets, 29, 121 Duncan M., Quinn T., Tremaine S., 1987, AJ, 94, 1330 Ferna´ndez J. A., 1981, A&A, 96, 26 Ferna´ndez J. A., 2010, in J. A. Ferna´ndez, D. Lazzaro, D. Prialnik, & R. Schulz ed., IAU Symposium Vol. 263 of IAU Symposium, The discovery rate of new comets in the age of large surveys. Trends, statistics, and an updated evaluation of the comet flux. pp 76–80 Ferna´ndez J. A., Ip W., 1991, in R. L. Newburn Jr., M. Neugebauer, & J. Rahe ed., IAU Colloq. 116: Comets in the post-Halley era Vol. 167 of Astrophysics and Space Science Library, Statistical and evolutionary aspects of cometary orbits. Kluwer Academic Publishers, Doordecht, the Netherlands, 1991., pp 487–535 Fouchard M., Froeschle´ C., Valsecchi G., Rickman H., 2006, Celestial Mechanics and Dynamical Astronomy, 95, 299 Gardner E., Flynn C., 2010, MNRAS, 405, 545 Grieve R. A. F., Pesonen L. J., 1996, Earth Moon and Planets, 72, 357 Heisler J., Tremaine S., 1986, Icarus, 65, 13 Hills J. G., 1981, AJ, 86, 1730 Holmberg J., Flynn C., 2000, MNRAS, 313, 209 Holmberg J., Flynn C., 2004, MNRAS, 352, 440 Hughes D. W., 2001, MNRAS, 326, 515 Jetsu L., Pelt J., 2000, A&A, 353, 409 Juric´ M., Ivezic´ Zˇ., Brooks A., Lupton R. H., Schlegel D., et al. 2008, ApJ, 673, 864 Korchagin V. I., Girard T. M., Borkova T. V., Dinescu D. I., van Altena W. F., 2003, AJ, 126, 2896 Leto G., Jakub´ık M., Paulech T., Neslusˇan L., Dybczyn´ski P. A., 2008, MNRAS, 391, 1350 Lyytinen J., Jetsu L., Kajatkari P., Porceddu S., 2009, A&A, 499, 601 Matese J. J., Innanen K. A., Valtonen M. J., 2001, in M. Y. Marov & H. Rickman ed., Astrophysics and Space Science Library Vol. 261 of Astrophysics and Space Sci- ence Library, Variable Oort cloud flux due to the Galactic tide. pp 91–102 Matese J. J., Whitman P. G., 1989, Icarus, 82, 389 Matese J. J., Whitman P. G., 1992, Celestial Mechanics and Dynamical Astronomy, 54, 13 Matese J. J., Whitman P. G., Innanen K. A., Valtonen M. J., 1995, Icarus, 116, 255 Matese J. J., Whitmire D., 1996, ApJL, 472, L41+ Mikkola S., Nurmi P., 2006, MNRAS, 371, 421 Napier W. M., 2006, MNRAS, 366, 977 Nurmi P., 2001, MNRAS, 323, 911 Oort J. H., 1950, BAN, 11, 91 O¨pik E. J., 1966, Irish Astronomical Journal, 7, 141 Rampino M. R., Stothers R. B., 1984, Nature, 308, 709 Rickman H., Fouchard M., Froeschle´ C., Valsecchi G. B., 2008, Celestial Mechanics and Dynamical Astronomy, 102, 111 Rohde R. A., Muller R. A., 2005, Nat, 434, 208 Scho¨nrich R., Binney J., Dehnen W., 2010, MNRAS, 403, 1829 Stothers R. B., 2006, MNRAS, 365, 178 Wickramasinghe J. T., Napier W. M., 2008, MNRAS, 387, 153 Wiegert P., Tremaine S., 1999, Icarus, 137, 84 Yabushita S., 2004, MNRAS, 355, 51 c© 2010 RAS, MNRAS 000, 1–8