Table Of Content

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed2February2008 (MNLATEXstylefilev2.2) The stellar velocity distribution in the solar neighbourhood ⋆ Richard S. De Simone, Xiaoan Wu and Scott Tremaine Princeton University Observatory, PeytonHall, Princeton, NJ 08544-1001, USA; e-mail: desimone [email protected], [email protected], [email protected] 4 0 2February2008 0 2 n ABSTRACT a J We explore the heating of the velocity distribution in the solar neighbourhood by 9 stochastic spiral waves. Our investigation is based on direct numerical integration of 2 initiallycirculartest-particleorbitsintheshearedsheet.Weconfirmthe conclusionof otherinvestigatorsthatheatingbyspiralstructurecanexplaintheprincipalfeaturesof 2 theage-velocitydispersionrelationandotherparametersofthevelocitydistributionin v 6 the solar neighbourhood.In addition, we find that heating by strong transientspirals 0 cannaturallyexplain the presence ofsmall-scalestructure in the velocitydistribution 9 (“movinggroups”).Heating by spiralstructure also explains why the stars in a single 0 velocity-space moving group have a wide range of ages, a result which is difficult to 1 understand in the traditional model that these structures result from inhomogeneous 3 star formation. Thus we suggest that old moving groups arise from irregularities in 0 the Galactic potential rather than irregularities in the star-formationrate. / h Key words: solar neighbourhood – Galaxy: kinematics and dynamics – Galaxy: p fundamental parameters – stars: kinematics – galaxies:kinematics and dynamics - o r t s a v: 1 INTRODUCTION σi2j ≡(vi−vi)(vj−vj); (2) i The velocity distribution function (df) of stars in thesolar X where X f(v)X(v)dv/ f(v)dv, and thematrix σ2 is neighbourhood provides unique insights into the Galactic ≡ ij r potential field, the dynamical history of the Galactic disk, theinverse oRf thematrix αijR. a In a steady-state, axisymmetric galaxy (i) the mean and the relationships between kinematics, age, and metal- velocity relative to the LSR is tangential, so v = v = licity for disk stars. x z 0; (ii) the tensors α and σ2 are both diagonal; (iii) Let us define the Local Standard of Rest (LSR) to be ij ij the ratio σ2 /σ2 is determined by the local gravitational afictitiouspoint that coincides with theSunat thepresent xx yy force and its radial gradient (e.g. Chandrasekhar 1942, instant and travels in a circular orbit around the Galactic Binney & Tremaine 1987; see also eq.14 below). centre.WeintroducearotatingCartesiancoordinatesystem The velocity dispersion of stars in the solar neigh- with origin at the LSR, x-axis pointing radially outward, bourhood increases with age, probably because the disk is y-axis pointing in the direction of Galactic rotation, and z- axispointingtothesouthGalacticpole.Thedff(v)isde- “heated” by one or more dynamical mechanisms. However, theinterpretationoftheobservedage-velocitydispersionre- finedsothatf(v)dv isthenumberofstarsperunitvolume avr avr withvelocitiesintherange[v,v+dv].Thestandardempir- lation( )isuncertain:(i)Oneschoolmodelsthe as a smooth power law, σ (t) tp, and interprets this be- ical model for the past century has been the Schwarzschild xx ∝ df haviorasevidenceforacontinuousheatingmechanism.Es- (1907) , which is a triaxial Gaussian of theform timates of the exponent p in the literature span the range 3 0.2–0.5. (ii) Some authors argue that the dispersion rises f(v)∝exp(cid:20)− 21 αij(vi−vi)(vj−vj)(cid:21), (1) steeplywith agefor stars <5Gyrold,andthereafterisrel- iX,j=1 atively flat (Carlberg et al∼. 1985; Gómez et al. 1997), per- where v (v ,v ,v ). The two lowest moments of the df haps because the continuous heating mechanism saturates ≡ x y z once the dispersion is large enough. (iii) A third model is are the mean velocity, v = (v ,v ,v ), and the velocity- x y z that the dispersion does not increase smoothly with age. dispersion tensor, Forexample,Freeman(1991)andQuillen & Garnett(2001) argue that there are three discrete age groups (t<3Gyr, 3Gyr<t<10Gyr, t>10Gyr) with different disp∼ersions, ⋆ Towhomcorrespondenceshouldbeaddressed. but w∼ithi∼n each grou∼p there is no evidence for a correla- 2 R. De Simone, X. Wu & S. Tremaine tion between dispersion and age. Such groups might arise gular momentum, since this determines their mean angular because the continuous heating saturates after only 3 Gyr, velocity). The existence and membership of these groups andthehigherdispersion oftheoldest starsis duetoadis- was controversial until the Hipparcos satellite measured re- crete event such as a merger. liable distances and proper motions for a large, homoge- A wide variety of mechanisms for disk heating has neousdatabaseof nearbystars, and verifiedthepresenceof df been discussed (see Lacey 1991 for a review): (i) Spitzer rich substructure in the velocity of both young and old & Schwarzschild (1951, 1953) suggested that massive gas stars,includinganumberoffeaturesthatcoincidewithmov- cloudscouldgravitationallyscatterstars,leadingtoasteady ing groups already identified by Eggen (e.g. Dehnen 1998, increase in velocity dispersion with age, and thereby pre- Chereul et al. 1998,1999) dictedtheexistenceofgiantmolecularclouds(gmcs).How- In this paper we explore a quite different explanation ever, heating by gmcs alone fails to explain several ob- for substructure in the velocity df. We suggest that sub- servations (Lacey 1984; Villumsen 1985; Jenkins& Binney structurearisesnaturallyfromthesamespiralgravitational 1990; Lacey 1991; Jenkins 1992): the predicted ratio fluctuationsthatexcitethegrowthofthevelocitydispersion. σ /σ of vertical to radial dispersion may be too high, Inotherwords,substructureiscausedbyhomogeneousstar zz xx roughly 0.72 compared to the observed value of 0.5 for formationinanirregularpotential,asopposedtoinhomoge- old stars (but see Ida, Kokubo& Makino 1993 for an neousstarformationinaregularpotentialinthetraditional opposing view); the predicted exponent in the avr is model. somewhat too low, p<0.25; and the masses and num- Weinvestigatethishypothesisbysimulatingtheevolu- ber density of gmcs∼determined from CO observations tionofthevelocitydfinducedbytransientspiralstructure are too low to explain the observed dispersion, proba- in a simple two-dimensional model of the local Galaxy. We bly by a factor of five or so but with substantial uncer- restrictourselvestotwodimensionsbecausespiralstructure tainty. (ii) Transient spiral waves lead to potential fluctu- does not excite vertical motions efficiently,and because the df ations in the disk that excite the random motions of disk velocity appearstobewell-mixedintheverticaldirection stars(Barbanis & Woltjer1967;Sellwood & Carlberg 1984; (Dehnen1998). Carlberg & Sellwood 1985). However, spiral waves only ex- FollowingEggen,weshallusetheterm“movinggroup” cite the horizontal (x and y) velocity components effec- to denote substructure in the velocity df of old stars (>1 tively, since their characteristic spatial and temporal scales Gyr) at a given position. Unfortunately, the same term∼is are much larger than the amplitude or period of verti- sometimes also applied to OB associations, which are spa- cal oscillations for young stellar populations. (iii) These tially localized concentrations of much younger stars (e.g., considerations lead naturally to a hybrid model, in which deZeeuw et al. 1999). spiral waves excite non-circular velocities in the plane, Section3describesoursimplifieddynamicalmodel.The and the velocities are then redistributed between horizon- results of our simulations are analyzed and compared to tal and vertical motion through gmc scattering (Carlberg observed data in 4. Section 5 examines briefly the tradi- § 1987). The hybrid model has been investigated thoroughly tional hypothesis that substructure arises from inhomoge- by Jenkins & Binney (1990) and Jenkins (1992) using the neous star formation. Section 6 contains a brief discussion Fokker-Planck equation. They find that they can success- of the closely related process of radial migration of stars, fully reproduce the observed axis ratio σ /σ , the expo- and 7 contains concluding remarks. avr zz xx § nentpinthe andtheradialdispersionofoldstars.(iv) Otherpossibleheatingmechanisms,allofwhichrelytosome extentonhypotheticalorpoorlyunderstoodcomponentsor features of the Galaxy, include scattering by massive com- 2 PROPERTIES OF SPIRAL STRUCTURE pact halo objects or halo substructure (Lacey & Ostriker 1985), mergers with dwarf galaxies (Tóth & Ostriker 1992; OurmodeldependsonseveralpropertiesoftheGalaxy’sspi- Walker, Mihos & Hernquist1996;Huang& Carlberg1997), ralstructure,suchasthenumberofarms,thearmstrength, or the outer Lindblad resonance from the Galactic bar andthepitchangle. Weareconcerned withspiral structure (Kalnajs1991;Dehnen1999,2000;Fux2001;Quillen2003). inthedisksurfacedensity(ratherthan,say,inthedistribu- df TheSchwarzschild (1)onlyapproximatestheveloc- tion of youngstars or gas). Theproperties of this structure ity distribution on the largest scales in velocity space. On in external galaxies are revealed by near-infrared images, smaller scales, there is substructure, which is most promi- which are dominated by the stars that contribute most of nent in the youngest stars but present in stars of all ages. the mass (Rix & Rieke 1993, but see also Rhoads 1998 and df Discussion of substructure in the velocity dates back James & Seigar 1999 for qualifications). to Kapteyn’s (1905) model of “two star streams”, and for Rix & Zaritsky (1995) examined the K-band spiral decades Eggen advocated the case for substructure in the structureof18face-ongalaxies.Theyfoundthatalmosthalf formof“movinggroups”inthesolarneighbourhood(Eggen hadstrongtwo-arm(m=2)spirals,witharm-interarmcon- 1996 and references therein). trasts I /I 2. If the arms are sinusoidal, with frac- max min ≃ Eggenandothershaveusuallyexplainedmovinggroups tionalamplitudeǫrelativetotheaxisymmetricbackground, as the result of inhomogeneous star formation in the disk: thearm-interarm contrastisI /I =(1+ǫ)/(1 ǫ),so max min − in this model, stars in a moving group are born at a com- a contrast of 2 corresponds to ǫ 0.3. ≃ mon place and time, and then disperse into a stream that Block & Puerari (1999) examined K-band spiral struc- happens to intersect the solar neighbourhood. This model ture in 19 spirals. They found pitch angles ranging from 8◦ predictsthatstars in amovinggroup should havethesame to 49◦ with a median value of 22◦, and m = 2 amplitudes age, metallicity, and azimuthal velocity (i.e. the same an- ranging from ǫ=0.03 to ǫ=0.5, with a median ǫ=0.1. Velocity distribution in solar neighbourhood 3 Seigar & James (1998) conducted a similar survey of The solutions of theequations of motion are 45 galaxies. They found that the dominant Fourier mode x = x +acos(κt+φ), usually has m= 1 (36%) or m= 2 (31%), and the median g pitch angle of the spiral structure was 8◦, with little or no y = y (t) 2Ωasin(κt+φ), correlation with Hubble type. They measure the strength g − κ of the spiral arms in terms of the “equivalent angle”; their yg(t) = yg0 2Axgt, (6) − medianequivalentangleof14◦ correspondstoanamplitude where [x ,y (t)] are the coordinates of the guiding centre, g g ǫ 0.07 for a sinusoidal m=2 spiral. ≃ a is the epicycle amplitude, φ is a phase constant, and κ is Elmegreen et al. (1999) stress that the amplitude of theradial or epicycle frequency, m = 2 spiral structure in the near-infrared depends on whethertheopticalspiralarmsareclassifiedasflocculentor κ=Ω 4+2dlnΩ 1/2 or κ2 =4Ω(Ω A). (7) grand-design (Elmegreen 1998). Grand-design spirals have (cid:16) dlnR(cid:17)R0 − arm-interarmcontrastsof1.5–6,correspondingtoamplitude The energy and angular momentum are related to the ǫ = 0.2–0.7, while flocculent galaxies have contrast <1.7, guiding-centreradius and epicycle amplitudeby corresponding toǫ<0.25. ∼ At optical wav∼ebands, Ma et al. (1999) find that the E =2A(A−Ω)x2g+ 12κ2a2, H =2(Ω−A)xg, (8) meanpitchanglefor51Sbcgalaxies(thesameHubbletype and the guiding-centre radius is related to the phase-space as the Galaxy) is 15◦. coordinates by ThemeasurementofspiralstructureinourownGalaxy is more difficult than in face-on external galaxies. Drimmel y˙+2Ωx x = . (9) (2000) uses K-band photometry of the Galactic plane to g 2(Ω A) − conclude that the Galaxy contains a two-arm spiral with The epicycle energy is definedas pitchangle18◦.Vallée(2002)reviewsanumberofstudiesof the Galaxy’s spiral structure, mostly based on young stars, Ex ≡ 21 x˙2+κ2(x−xg)2 gasanddust;theseyieldarangeofpitchanglesfrom6◦–20◦, (cid:2) AH2 (cid:3) but Vallée concludesthat thebest overall fit is providedby = E+ 2(Ω A) an m=4 spiral with pitch angle of 12◦. The arm structure − = 1κ2a2 in the Galaxy appears to be intermediate between grand- 2 design and flocculent (Elmegreen 1998). 2Ω2 = 1x˙2+ (y˙+2Ax)2. (10) 2 κ2 Theepicycleenergyiscloselyrelatedtotheradialaction I = 1κa2 =E /κ. (11) 2 x 3 A NUMERICAL MODEL OF DISK HEATING Fora particle in a circular orbit, We model disk heating in the sheared sheet, which approx- imates the local dynamics of a differentially rotating disk (x,y,x˙,y˙)=(xg,yg0 2Axgt,0, 2Axg), (12) − − (Spitzer& Schwarzschild 1953; Goldreich & Lynden-Bell thusE =0. x 1965;Julian & Toomre1966).TheLSRisassumedtotravel The approximations used in deriving the linearized around the Galactic centre in a circular orbit of radius R 0 equationsofmotion(3)areonlymarginallyvalidinthesolar and angular frequency Ω > 0. We use the same Cartesian neighbourhood,sincetheepicycleamplitudesacanbeasig- coordinate system (x,y,z) introduced in the last Section, nificantfractionofR (forapopulationwithradialvelocity 0 restrict ourselves to the z = 0 plane, and denote position dispersion σ in a galaxy with a flat rotation curve, xx and velocity by x=(x,y) and v=(x˙,y˙). For x, y R 0 theequations of motion of a test particle are | | | |≪ a2=σx2x/Ω2. (13) ∂Φ Thusold disk stars in the solar neighbourhood, with σxx x¨−2Ωy˙−4ΩAx = −∂x, 40kms−1, have(a2)1/2 0.18R0). Nevertheless, we believ≃e ≃ ∂Φ that thesheared sheet accurately capturesthemost impor- y¨+2Ωx˙ = , (3) −∂y tant features of theevolution of disk-starkinematics. The kinematicsof a population of stars is described by where Φ(x,t) is the gravitational potential due to sources thedff(x,v,t),wherefdxdvisthenumberofstarsinthe otherthantheaxisymmetricGalacticdisk,andA>0isthe interval(x,v) (x+dx,v+dv).AccordingtoJeans’sthe- Oort constant → df orem (Binney & Tremaine 1987), a stationary can only 1 dΩ depend on the integrals of motion E and H (or Ex, xg, a, A= R , (4) −2 dR etc.). In the solar neighbourhood (x=y =0), the integrals (cid:16) (cid:17)R0 areH =y˙ andE = 1x˙2+2Ω2y˙2/κ2.Thusinasteadystate whereRΩ(R)isthecircularspeedatradiusRandΩ(R0)= thevelocity distrxibuti2on must bean even function of x˙. Ω. Itisstraightforwardtoshowthatthemeanvelocityand IfΦ=0thetrajectoriesgovernedbyequations(3)have velocity-dispersiontensorofanystationary,spatiallyhomo- two integrals of motion related to energy and angular mo- geneous df in thesheared sheet must satisfy therelations mentum, σ2 Ω E ≡ 21(x˙2+y˙2−4ΩAx2), H ≡y˙+2Ωx. (5) vx =0, vy =−2Ax, σx2y =σy2x =0, σxy2yx = Ω−A. (14) 4 R. De Simone, X. Wu & S. Tremaine df A useful model for the sheared sheet is ω=2Ak x κ. (20) y g ± f(x,v) exp E /σ2 , (15) We focus instead on transient spiral structure, which can ∝ − x 0 (cid:0) (cid:1) heat stars over a range of radii. We model each transient where E is defined in equation (10). At a given location, x spiral using a Gaussian amplitude dependence centred at df the (15)leadstoatwo-dimensionalversionofthetriaxial time t with standard deviation σ : s s Gaussianvelocitydistribution(1),inwhichthemeanveloc- 2πǫGΣ (t t )2 iTtyhuasntdhveedlofci(t1y5-d)iisspaelrssoiosnomteentsiomressa(tcisofnyfuthsiengrelyla)tcioanllsed(1t4h)e. Φs(t)= k d exp(cid:20)− −2σs2s +i(θ+2Akyxct)(cid:21). (21) df Schwarzschild . Here θ is a phase constant, and x is the corotation radius c In our simulations the relations (14) are not satisfied relativetotheLSR(thephaseofthespiralisconstantforan exactly,soitisusefultointroducetheprincipal-axissystem observeronacircularorbitatx ).Theparameterǫmeasures c (x1,x2) in which the velocity-dispersion tensor is diagonal; the amplitude of the spiral, and the normalizing constants the x -axis is chosen to lie within 45◦ of the x-axis. The 1 arechosensothatthemaximumsurfacedensityinthespiral corresponding velocity dispersions are σ and σ , which we 1 2 is a fraction ǫ of the surface density of the underlying disk will normally plot instead of σxx and σyy; usually in our (eq.19). simulations andin thedataweshall findthat σ >σ .The 1 2 Ineachsimulationthetrajectories ofthestarswerefol- vertex deviation l is the Galactic longitude of the x -axis, v 1 lowed between time t = 0 and t = t . During this interval 0 and is given by thestarswereperturbedbyN transientspiralwaves.Each s 1 2σ2 wave had phase constant θ chosen randomly from [0,2π], lv ≡−2arctan(cid:18)σx2x xyσy2y(cid:19), (16) corotation radius xc chosen randomly from a Gaussian dis- − tribution centred on the LSR with standard deviation σ , c wbohuerrheo|olvd|h<ave45v◦e.rtSetxeldlaervipaotipounlsatiinontsheinratnhgees0ol<arl n<ei3g0h◦- (atnhdisciesntsrliaglhttilmyelotnsgcehrotsheannrathnedoimntleygrfarotimon[−in2tσersv,at0l,+to2σins-] v (Binney & Merrifield 1998). ∼ ∼ clude the effects of transients whose wings are inside the Throughout the paper, we assume that the rota- integration interval although theircentral times are not). tional curve of the underlying axisymmetric galaxy is flat, The fractional amplitude of the surface-density per- RΩ(R) = constant, so that A = 1Ω and κ = √2Ω. We turbation due to a single transient spiral is ǫ at the wave 2 shall assume that the surface density of the disk in the peak. However, in some ways a better quantity to compare solar neighbourhood is Σd = 50M⊙pc−2; recent observa- withtheobservationaldataonspiralamplitudesistheroot- tional estimates are Σd = 40M⊙pc−2 (Crézé et al. 1998), mean-square (rms) time average of the fractional surface- 48M⊙pc−2 (Holmberg & Flynn 2000), and 42 6M⊙pc−2 density amplitude, ± (Korchagin et al. 2003). The assumed surface density only π1/2N σ 1/2 enters our analysis through the definition of the fractional ǫ =ǫ s s . (22) spiral amplitude ǫ below (eq.21). rms (cid:18) t0 (cid:19) rms Acloselyrelatedquantityisthe potentialperturbation, sinα 3.1 Spiral waves Φrms =2πGǫrmsΣdR0|√2m|; (23) Weapproximatespiralstructureasasuperpositionofwaves (the factor √2 arises because Φ is the rms fluctua- with surface density and potential rms rms tion rather than the amplitude, which is smaller by Σ(x,y,t) = Σs(t) expi(kxx+kyy), √2). Jenkins& Binney (1990) use Fokker-Planck calcula- Φ(x,y,t) = Φ (t) expi(k x+k y), (17) tions of heating by transient spiral structure to estimate s x y that Φ =(9–13kms−1)2. rms where k = (kx,ky) is the wavenumber, and only the real The power spectrum of each transient is a Gaussian partofΣorΦisphysical.Withoutlossofgeneralitywecan with standard deviation (2σ2)−1/2. The central frequencies s assume thatky 0; thenthespiral istrailing if kx >0 and ofthepowerspectraofthetransientsfollowaGaussiandis- ≥ leading if kx < 0. In this paper, we only consider trailing tribution with standard deviation 2Akyσc. Thus the power spiral waves. The number of arms m and the pitch angle α spectrum of the ensemble of N transients is smooth if the s are determined by k through therelations overlap factor ky = Rm0, (cid:12)(cid:12)kkxy(cid:12)(cid:12)=tanα. (18) C ≡Ns(22σAs2k)y−σ1c/2 = 23/2ANksyσcσs (24) Therelationbetwee(cid:12)np(cid:12)otentialandsurfacedensityforspiral is large compared to unity; on the other hand if C 1 the waves is given by(e.g. Binney & Tremaine 1987) ≪ powerspectrumandhencetheheatingislocalizedatnarrow 2πGΣ (t) resonances.AllofoursimulationshaveC 1(seeTable1). Φs(t)=− ks , (19) The dispersion of the power spectrum≫of all the transients is σ , given by where k= k. t | | 1 starsStoenadnyeasrplyiraclisr,cuinlawrhoircbhitΣsso(nt)ly,Φast(tt)h∝e Leixnpd(biωlatd),rheesaot- σt2= 2σs2 +(2Akyσc)2. (25) nances (Lynden-Bell& Kalnajs 1972), which occur at the Note that in this model problem different azimuthal guiding centreradii given by wavenumbers m and m′ = fm are equivalent if we rescale Velocity distribution in solar neighbourhood 5 the other variables by: k′ = fk, α′ = α, t′ = t, (x′,y′) = f(v)= F δ(v v ) (26) i i (x/f,y/f), Φ′s = Φs/f2, Σ′s = Σs/f, ǫ′ = ǫ/f, σi′j = σij/f Xi − σ′ =σ /f, σ′ =σ /f, σ′ =σ . c c 0 0 s s is replaced by 1 (v v )2 f(v)= 2πσ2 Fiexp(cid:20)− 2−σ2i (cid:21), (27) 3.2 Determining the solar neighbourhood velocity ob Xi ob distribution where σ is theassumed observational error. ob df Assume that the disk stars are formed on circular orbits at Wehavedescribedtheprocedureforestimatingthe time0.Wewishtodeterminethevelocitydistributionatthe of stars of age t0, using equation (15) and the epicycle en- df present time t in the solar neighbourhood, (x,y) = (0,0). ergies at t =0. However, we can also determine the for 0 This is atwo-point boundary-valueproblem rather thanan stars having an intermediate age t0 tm in the course of − initial-valueproblem:theboundaryconditionsonthephase- thesameorbitintegration,usingtheepicycleenergiesofthe stp=ace0 (coonordciirncautleasr aorrebit(sx),ya,nxd˙,y(˙)x,=y) =(xg(,0y,g00),0a,t−t2=Axtg)(aint tsthaerdseaptenthdeeninceteorfmtehdeiadtfeotnimsetetllmar.aTghei,saanpdptrhouaschdemteeramsuinreess 0 avr the solar neighbourhood). The solutions to the boundary- the for a given realization of theset of transient spiral value problem will be a set of distinct initial positions arms, using only a single set of orbit integrations. [x (t = 0),y (t = 0)] or final velocities [x˙ (t = t ),y˙ (t = i i i 0 i t )] ( u ,v ),i=1,2,...(thesignofuischosentoagree 0 ≡ − i i 3.3 Units and model parameters with the usual convention that stars moving towards the Galactic centre haveu>0). Because themotion of stars in The time unit is chosen to be Ω−1 and the distance unit the fluctuating gravitational field of the transient spirals is is chosen to be R . The azimuthal orbital period is then 0 complicated, we expect that there will be many solutions. 2π/Ω=2π.Forconversiontophysicalunitsweassumethat In this idealized model the velocity-space df in the solar R =8kpcand ΩR =220kms−1, so that the unit of time 0 0 neighbourhood willconsist ofaset ofdelta-functionsat the is 35.6 Myr. Most of our numerical integrations ran for an velocities (ui,vi),butin practicethesewill besmeared into interval t0 =281, corresponding to 40.8 rotation periods of a continuous distribution by observational errors, the non- theLSR or 10.0 Gyr. zero volume of the “solar neighbourhood”, the small initial We consider only trailing spirals with m=2 or m=4. velocity dispersion of thestars when they are formed, etc. The standard deviation of the distribution of corotation Normally,boundary-valueproblemsaresolvedmostef- radii for the transient spirals was σ = 0.25 correspond- c ficiently by iterative methods. In this case, however, most ing to 2 kpc. With these choices the overlap factor is of the solutions have only a tiny domain of attraction in C =1.41N (m/2)/(Ωσ ). s s the (u,v) plane; thus iterative methods are inefficient. We The radial velocity dispersion of stars at the time of havethereforeadoptedadifferentapproach,basedonMonte their birth is taken to be σ = 3kms−1. The observational 0 Carlo sampling (Dehnen2000). error in the velocities is taken to be σ =3kms−1; for the ob Stars are not born on precisely circular orbits, in part typical Hipparcos parallax error of 1 mas, σ corresponds ob because star-forming clouds are not on circular orbits. to the velocity error arising from the distance uncertainty We may therefore assume that the stars initially have a forastarwithtransversevelocity30kms−1 atadistanceof df Schwarzschild (eq. 15), with a small but non-zero ini- 100 pc. tial velocity dispersion σ0. In each simulation we launch either 1 107 or 7 107 × × Following Dehnen, we integrate orbits backward in stars from the solar neighbourhood and integrate back for time,startingatt=t (whichwecalltheinitialintegration 10 Gyr (the smaller simulations are used to estimate mo- 0 time or “present epoch”) and ending at t = t (the final ments,suchasthoseinFigure1,andthelargersimulations 0 integration time or “formation epoch”). The initial condi- are used to estimate the df, as in Figure 7). In the larger tions are chosen at random from (x,y,x˙,y˙) = (0,0, u,v), simulations, typically 10000–15000 stars contribute signifi- with u, v < v . We choose v = 110kms−1−, large cantlytothesolution,inthesensethattheirepicycleenergy max max | | | | enough to include almost all disk stars in the solar neigh- attheformationepochis<3σ2.Thisnumberislargerthan 0 bourhood. We then integrate equations (3) backwards in the number of stars in the Hipparcos samples to which we time to the formation epoch. The collisionless Boltzmann compare(e.g.4600inthelargestcolor-selectedsampleanal- equation (Binney & Tremaine 1987) states that the phase- ysedbyDehnen1998).Wehaveverifiedthattheconclusions space density around a trajectory is time-invariant. There- belowareunaffectedifwedoublethenumberofstarsinthe fore the phase-space density F at (u ,v ) at the present simulation or change the random-number seed used to se- i i i epoch is given by (15), where the epicycle energy E is lectthesampleofstars(keepingthesameseedtoselect the x measured at theformationepoch. Weassume that thespa- spiral transients). tial volume of our survey of the solar neighbourhood is in- We carried out a number of simulations with different dependent of velocity, so the density in velocity space at valuesofthespiral-waveparametersm(numberofarms),N s thepresentepochisproportionaltoF .Dehnen’sprocedure (number of spirals), ǫ (fractional surface-density perturba- i therefore provides a Monte Carlo sampling of the velocity- tion),σ (duration of thetransient),σ (dispersion in coro- s c df space at thepresent epoch. tation radius), and α (pitch angle). The parameters of the We convert our Monte Carlo realization to a smooth simulationsareshowninTable1.Eachsimulationisdefined df byconvolutionwithobservationalerrors.Wereplacethe by these six parameters and the random number seeds for point solutions v =(u ,v ) by Gaussians, that is, both thestars and thespiral arms. i i i 6 R. De Simone, X. Wu & S. Tremaine Table 1.Simulationparameters Parameter Symbol Run˚5 Run˚10 Run˚20 Run˚40 Runm˚4 Run˚sim Definingparameters Numberofspiralarms m 2 2 2 2 4 2 Numberofspiralwaves Ns 45 18 12 9 12 48 Maximumfractionalsurface densityofspiralwaves ǫ 1.61 0.80 0.80 0.67 0.80 0.31 Timescaleofspiralwaves σs 1 1 1 1 1 1 rms corotationradius(R−R0)/R0 σc 0.25 0.25 0.25 0.25 0.25 0.25 Pitchangle α 5◦ 10◦ 20◦ 40◦ 20◦ 20◦ Derivedparameters Overlapfactor(eq.24) C 64 25 17 13 34 68 rms surface-densityamplitude(eq.22) ǫrms 0.86 0.27 0.22 0.16 0.22 0.17 rmspotential in(kms−1)2 (eq.23) Φrms (16.9)2 (13.4)2 (17.0)2 (19.8)2 (12.0)2 (15.0)2 Our runs have pitch angles of 5◦, 10◦, 20◦, and 40◦, model in which the Galaxy has experienced several short- whichspantherangeofobservedpitchanglesinspiralgalax- lived grand-design spiral phases, possibly caused by minor ies.Wethenchoosetheamplitudeǫandthenumberoftran- mergers. In contrast, run˚sim has a much larger number of sients N to approximately reproduce the observed radial weaker transients (ǫ=0.3, N σ /t =0.17). s s s 0 velocity dispersion of old main-sequence stars in the solar Each run was repeated seven times (realizations a-g), neighbourhood1. using different random-number seeds to generate both the All ofourrunshavetwo-armed spirals (m=2), except stars and the spiral transients, in order to distinguish sys- forrunm˚4,whichhasthesameparametersasrun˚20except tematicfromstochasticvariations.Tostudyextremelylong- that m=4. termbehavior,wealsoranarealizationh,inwhichboththe We also have run˚sim, in which we have used a larger integration time and the number of transients are a factor number of weaker transients than in run˚20, and have also of three larger than in a-g. adjusted some of the central times of the spiral transients tobetterreproducethefeaturesoftheobservedsolarneigh- df bourhood (see 4.4). 4 RESULTS Thermspote§ntialfluctuationΦ =(13–20kms−1)2 rms in our four runs, higher than the velue Φ = 4.1 Age-velocity dispersion relation rms (9–13kms−1)2 derived by Jenkins & Binney (1990) from In 3.2,wehaveshownthatasinglesetoforbitintegrations Fokker-Planck calculations by about a factor of 2–2.4. § df can be used to derivethe for stars of all ages between 0 A minor component of this difference arises because df and t . Thus we may combine s of different ages to pro- gmc 0 Jenkins & Binney (1990) include s, which contribute df duce the that corresponds to a uniform distribution of about 20% of their total heating rate. A more important df ages between 0 and 10 Gyr (hereafter the “combined ”), effect, pointed out to us by A. J. Kalnajs, is that in our aswouldresultfrom auniformstar-formation rate.Theav- model the heating is spatially inhomogeneous: our assump- df eragemomentsofthis forourrunsareshowninTable2, tions place the Sun fairly close to corotation for most of along with the observed values taken from Dehnen (1998). thetransients,sothatheatingattheLindbladresonancesis Ofcourse,theclose agreement oftheradial velocitydisper- more effective inside and outside the solar circle than it is sion σ in each run with the observed value of 37kms−1 xx at thesolar circle. arises because the strength and numberof spiral transients The fractional surface-density fluctuation ǫ required to were chosen to producethecorrect radial dispersion. reproduce the avr in run ˚5 exceeds unity, and the rms The evolution of the principal axes σ and σ of the 1 2 surface-densityfluctuationisclosetounity;thesehighvalues velocity-dispersiontensorinruns˚5,˚10,˚20and˚40isshownin are unrealistic and suggest that a larger number of weaker Figure1.Thethreerowsofpanelsshowtheresultsforthree transients is required if thepich angle is as small as 5◦. different realizations, i.e., three different sets of random- Runs˚10,˚20,˚40,andm˚4havepeakamplitudesǫ=0.7– numberseedsusedtogenerateboththestarsandthespiral 0.8, similar to the strongest grand-design spirals, but these transients. In general, the results in this and other figures large amplitudesareonlypresentfor asmall fraction ofthe areindependentoftherealizationofthestellardistribution, time(N σ /t =0.03–0.06). Therelatively smallnumberof s s 0 andallofthevariationsseenareduetochangesinthereal- high-amplitudetransients in these runsis consistent with a ization of the transient spiral structure—a reflection of the fact that the Monte Carlo simulation of the stellar distribution is a numerical method while the simulation of the 1 More precisely, main-sequence stars redward of the Parenago distribution of spiral transients is a model of a stochastic discontinuity at B − V = 0.6 are all believed to have the physical process. avr same mean age, and have a color-independent radial dispersion The sillustrate the following: ≃38kms−1 (Dehnen&Binney 1998). Wematch this tothe radial dispersion in our simulation by assuming a constant star- Openspirals(largepitchangle)heatthediskstarsmore • formationrate,whichisconsistentwithsolarneighbourhooddata effectivelythantightlywoundspirals.Forexample,boththe withinthelargeuncertainties (Binney,Dehnen&Bertelli2000). numberNs andthesurface-densityamplitudeǫofthespiral Velocity distribution in solar neighbourhood 7 Figure 1.The age dependence ofthe principal axes of the velocity-dispersiontensor inruns˚5,˚10,˚20 and˚40. The solidlineshows σ1 andthedashedlineshowsσ2.Thethreerows(a,b,c)representdifferentrealizations.Thedot-dashedlinesshowthebestfitofequation (29)foreachrealization;seealsoTable3. Table 2.Simulationresultsandobservations Quantity Symbol Run˚5 Run˚10 Run˚20 Run˚40 Runm˚4 Run˚sim Observations Radialvelocitydispersion σxx(kms−1) 36±1 38±2 40±4 39±6 35±3 40±3 36.9 vertexdeviation lv(◦) −3±4 3±3 8±16 −8±12 3±8 −6±11 10.9 axisratioofvelocityellipsoid σ1/σ2 1.3±0.1 1.4±0.1 1.2±0.1 1.3±0.1 1.4±0.1 1.4±0.2 1.48 meanradialvelocity u(kms−1) 0±4 −1±3 0±5 3±1 0±4 0±8 −0.5 meanazimuthalvelocity v(kms−1) 0±2 0±2 −1±3 3±1 −2±3 2±5 −17.7 Resultsgivenarethemeanandstandarddeviationfromrealizationsa–g,forapopulationofstarswithagesuniformlydistributedbetween 0and10Gyr.“Observations”columngivesresultsfromsamplesB4+GIfromDehnen(1998),whichcontainsstarsredwardoftheParenago discontinuity(B−V >0.6,forwhichthelifetimeofmain-sequencestarsexceedstheageoftheGalaxy)andgiantstars.Thelargenegative meanazimuthalvelocityintheobservations isduetotheasymmetricdrift,whichisnotpresentinoursimulations. 8 R. De Simone, X. Wu & S. Tremaine transients are larger in ˚5 than the other runs shown (see run for t = 30 Gyr, three times as long as the other 0 Table 1), even though all runs produce the same velocity runs. The purpose of this long run was to see how well the dispersion by design. parametrization (29) worked as the heating process contin- Theheating occurs in discrete stepsseparated by peri- ued.Thedashedlinesshowthatequation(29)stillprovides • odsofzerogrowth.Thisfeaturearisesbecausewehaveonly areasonable fit;however,thebest-fitexponentp(seeTable a small number of spiral transients over the course of the 3) is systematically smaller than in the shorter realizations simulation, each of which hasthesame amplitudeandlasts a–g,exceptinrun˚5.Therefore,equation(29)shouldbecon- only for a characteristic time σ = 1 out of the total inte- sideredasafittingfunctionvalidoveralimitedtimeinterval, s gration timeof281. Duringtheintervalsbetweenthespiral ratherthan a formula with predictive power. transients, the df is stationary so the velocity dispersions All of our runs have the same value of the parameter are constant. σ (Table 1), which represents the duration of the spiral s Inrun˚5thegrowthofthedispersionslowssharplyafter transients. We have experimented with a range of values of the•radial dispersion reaches 30kms−1. This effect prob- σ . As shown in Figure 3, the dispersion for the combined s ∼ df ably arises because theepicycle amplitudebecomes compa- with a uniform distribution of ages between 0 and 10 rable to theradial spacing between thespiral arms, so that Gyrisalmostindependentofσ solongasσ >0.5.Aprob- s s the forces from the spirals tend to average to zero over an able explanation is that the width of the po∼wer spectrum epicycle oscillation. More specifically, the radial spacing of of the transient perturbations, given by equation (25), is the spiral arms is ∆x = 2πR tanα/m; the rms epicycle dominated by the dispersion in pattern speeds rather than 0 amplitude (a2)1/2 is given by equation (13), and we expect the duration of a single transient once σ >1/(2√2Ak σ ), s y c thespiralheatingtobecomeineffectivewhenthermsepicy- which occurs when σ >1.4 for our param∼eters. Of course, s cle amplitude exceeds the peak-to-trough distance between when σ becomes very∼large, the overlap factor of equation s arms, that is, when (a2)1/2/∆x>0.5. For small pitch angle (24), will decrease below unity and the heating will be lo- α, ∼ calized at discrete resonances. (a2)1/2 mσ m σ 5◦ = xx 0.50 1 . (28) ∆x 2πΩR0tanα ≃ 2 30kms−1 α 4.2 The velocity ellipsoid Therearesubstantialdifferencesamongtherealizations The time evolution of thevertex deviation is shown in Fig- • ofagivenrun;inotherwordstherearesignificantstochastic ure4. Not surprisingly, thestochastic variations in thever- effects in spiral-wave heating (see the discussion of the ex- tex deviation are larger for runs˚10,˚20 and˚40, which have ponentpbelowforfurtherdetail).Thestochasticvariations largerpitchangleandfewertransients,andthefluctuations intheheatingratebecomesmallerasthevelocitydispersion inthevertexdeviationdeclinewithtimeasthegrowthrate avr increases (seealso Figure2,whichfollows the soveran of the velocity dispersion declines. The small but non-zero even longer time interval). vertex deviations seen near the end of the runsare roughly consistent with the vertex deviation 10◦ seen in the old For comparisons with theoretical models, it is useful ∼ stellar population in the solar neighbourhood. There is no avr to fit the simulated s with a parametrized function. A evidence of a systematic preference for positive or negative common parametrization (e.g. Lacey 1991) is vertexdeviations for old stars. σ (t)=(σ1/p+Ct)p, (29) Figure 5 shows the evolution of the ratio of the ma- 1 0 jor and minor principal axes of the velocity ellipsoid. As in where σ0 is the initial dispersion along the major principal the case of the vertex deviation, the fluctuations are much axis, fixed at 3kms−1 as in 3.3. The motivation for this larger in the runs with larger pitch angle. In a steady-state § form is the assumption that the evolution of the dispersion axisymmetric galaxy with aflat rotation curve,theaxis ra- isgoverned byadiffusion equation with velocity-dependent tio should be σ /σ = [Ω/(Ω A)]1/2 = 21/2 = 1.414 (eq. 1 2 − diffusion coefficient, 14), and the observed axis ratio gradually approaches this valueas the runsprogress. dσ2(t) 1 =C σ−q(t), (30) Figure 6 shows the evolution of the mean radial and dt 1 1 azimuthal velocity in the solar neighbourhood. There are whereC =2pC andq=p−1 2,withC,pandqconstants. significant fluctuations, up to 10kms−1 for the larger pitch 1 − For σ σ , equation (29) implies σ (t) tp. angles, even in the oldest stellar populations. These fluc- 1 0 1 ≫ ∝ The dashed lines in Figure 1 show the results of least- tuations limit the accuracy with which we can construct squares fits for the parameters C and p in equation (29). axisymmetric models of local Galactic kinematics. In each case, the fitting function successfully preserves the Wehavealsoplottedtheaxisratioofthevelocityellip- avr overall shape of the but smooths overfluctuationsdue soidagainstthevertexdeviationbuthavefoundnosystem- to the finite number of transients. We fit each of the seven atic correlation between these two quantities. realizations of each runand obtained thevaluesof theavr Theparametersofthedfforapopulationofstarswith exponentpgiveninTable3.Thestandarddeviationinpcan a uniform distribution of ages between 0 and 10 Gyr are beaslargeas35%(inrun˚20)amongrealizationsofasingle shown in Table 2. They are all consistent with theobserva- parameterset,andthevariationduetochangesinthespiral- tions,exceptforthemeanazimuthalvelocityv.Thediscrep- structureparameters isjust aslarge. Thus,theexponentof ancy in v reflects thefact that oursimulations donot show avr the isnotlikelytodiscriminatewellbetweenheatingby asymmetricdrift becauseof theapproximations inherentin transientspiralstructureandheatingbyothermechanisms. the sheared sheet [eqs. (3) with the perturbing potential avr Figure 2 shows the s of realization h, which was Φ=0 are symmetric underreflection through theorigin]. Velocity distribution in solar neighbourhood 9 avr Figure 2.Theevolutionofthe over30Gyr.Thedashedlinesshowthebest-fitmodelsoftheform(29);seealsoTable3. df Figure 3.Thedependence ofthevelocitydispersiononthedurationofthetransients,derivedfromthecombined forapopulation of stars withauniform distributionof ages upto 10Gyr. Allof the runs plotted inagiven panel use the samedefining parameters as inTable 1, except for the transient duration σs. The solidand dashed lines show σxx and σyy, respectively. The velocity dispersionis almostindependent ofthetransientdurationσs,exceptforσs<∼0.5. 4.3 Shape of the velocity distribution 100pc(seeeq.27),fortwodifferentrealizationsofthestellar distributionandthesamerealizationofthespiraltransients. df The difference between these two realizations (the second Figure 7 shows the for stars with age 5 Gyr, from runs ˚5a,˚10a,˚20a and˚40a. We present the df in two forms: (i) and third rows) is an indication of the numerical noise in oursimulations, which is small. we display all of the points from our sampling of velocity df spaceatthepresentepochthatyieldnearlycircularorbitsat Thesimulated sin thesecond andthirdrow ofpan- the formation epoch (i.e. those with epicycle energy 3σ2, els do not resemble the Gaussian or Schwarzschild df (eq. σ = 3kms−1). (ii) We show contour plots of the s≤moot0h 1). They are much lumpier, and the lumps are distributed 0 df thatisexpectedwithHipparcosobservationalerrorsat more-or-less uniformly over a large area in velocity space, ∼ 10 R. De Simone, X. Wu & S. Tremaine avr Table 3.Fittedexponent pof s(eq.29). a b c d e f g a–ga h ˚5 0.31 0.20 0.22 0.27 0.27 0.22 0.31 0.25±0.05 0.24 ˚10 0.30 0.34 0.37 0.32 0.24 0.24 0.25 0.29±0.05 0.25 ˚20 0.33 0.46 0.35 0.59 0.29 0.76 0.57 0.48±0.17 0.29 ˚40 0.37 0.50 0.38 0.53 0.36 0.51 0.40 0.44±0.08 0.23 m˚4 0.48 0.28 0.26 0.44 0.24 0.53 0.26 0.36±0.12 - ˚sim 0.40 0.43 0.55 0.40 0.39 0.41 0.66 0.46±0.10 - a Themeanandstandarddeviationofthevaluesofpfromrunsa–g. Figure 4.Theevolutionofthevertexdeviation lv (eq.16).Thethreerowsrepresentdifferentrealizations. df rather than being concentrated at zero velocity. Qualita- thatobservationssuggestforourGalaxythe sarehighly df tively, the Schwarzschild s resemble a contour map of a irregular. single large mountain, whereas the simulated dfs resemble The dfs in Figure 7 exhibit prominent streaks in the anentiremountainrange.The“lumpiness”ordegreeofsub- direction v =v =constant, particularly at the larger pitch df y structure in the is a strong function of the pitch angle angles. The explanation of these streaks is straightforward, of thespiral arms, and in therange of pitch angles 10◦–20◦ and dates back to Woolley (1961). It is straightforward to