Mon.Not.R.Astron.Soc.000,1–9(2014) PrintedJanuary12,2015 (MNLATEXstylefilev2.2) The stellar kinematics of co-rotating spiral arms in Gaia mock observations 5 Jason A. S. Hunt1⋆, Daisuke Kawata1, Robert J. J. Grand2,3, Ivan Minchev4, 1 0 Stefano Pasetto1 and Mark Cropper1 2 n 1 Mullard Space Science Laboratory, UniversityCollege London, Holmbury St. Mary, Dorking, Surrey, RH5 6NT, UK a 2 HeidelbergerInstitut fur Theoretische Studien, Schloss-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany J 3 Zentrumfu¨rAstronomie derUniversita¨t Heidelberg, Astronomisches Recheninstitut, Mo¨nchhofstr. 12-14, 69120 Heidelberg, Germany 8 4 Leibniz-Institut fu¨rAstrophysik Potsdam (AIP), Ander Sternwarte 16, 14482 Potsdam, Germany ] A G SubmittedtoMNRAS:22nd December2014. . h p ABSTRACT - o We have observed an N-body/Smoothed Particle Hydrodynamics simulation of a r Milky Way like barred spiral galaxy. We present a simple method that samples N- t s bodymodelparticlesintomockGaiastellarobservationsandtakesintoaccountstellar a populations, dust extinction and Gaia′s science performance estimates. We examine [ the kinematics arounda nearbyspiralarmata similarpositionto the Perseusarmat 1 three lines of sight in the disc plane; (l,b)= (90,0),(120,0) and (150,0) degrees. We v find that the structure of the peculiar kinematics around the co-rotating spiral arm, 9 which is found in Kawata et al. (2014b), is still visible in the observational data ex- 6 pectedtobeproducedbyGaiadespitethedustextinctionandexpectedobservational 9 errorsof Gaia. These observable kinematic signatures will enable testing whether the 1 PerseusarmoftheMilkyWayissimilartotheco-rotatingspiralarmscommonlyseen 0 in N-body simulations. . 1 Key words: methods:N-bodysimulations—methods:numerical—galaxies:struc- 0 5 ture — galaxies: kinematics and dynamics — The Galaxy: structure 1 : v i X 1 INTRODUCTION density wave which can rotate rigidly as a feature with a r constant pattern speed and thusbe long lived. a The spiral features visible in many galaxies have long However, no N-body simulations have yet been able been the subject of debate. Although it has been almost to reproduce these long lived stable spiral arms, de- a century since the resolution of the “great debate” of spite the increase in computational power and resolu- Shapley & Curtis (1921), when it was argued over whether tion which has occurred in recent years (e.g. Sellwood these beautiful spiral structures were nebulae within our 2011; Dobbs& Baba 2014). Recent work has shown galaxyorgalaxiesintheirownright,themechanismswhich spiral modes and waves which survive over multi- generatethemarestilluncertain.Oneoftheproblemswith ple rotations (Quillen et al. 2011; Roškar et al. 2013; developing a comprehensive theory of spiral arms is the so Sellwood & Carlberg 2014) while the spiral arm features called “winding dilemma". It is known from observations of in the stellar mass are short-lived but recurrent (e.g. discgalaxiesthatthestarsintheinnerregionhaveahigher Sellwood & Carlberg 1984; Carlberg & Freedman 1985; angular velocity than those in the outer region. Therefore Bottema2003;Fujii et al.2011;Grand et al.2012a,b,2013; the spiral structure should “wind up" relatively quickly if Baba et al.2013;Roca-Fàbrega et al.2013;D’Onghia et al. the spiral arms rotate at the mean rotation velocity of the 2013) including in galaxies with a central bar (e.g. stars (e.g. Wilczynski 1896), contrary to observations of Grand et al.2012b),implyingthatthelargespiralarmsvis- many “grand design” spiral galaxies. A proposed solution ible in external galaxies may only appear to be rigid struc- tothewindingdilemmaisgivenbyspiraldensitywavethe- tures extending over the disc, while in fact being made of ory (Lin & Shu1964) which treats the spiral structureas a transient reforming features. Theinterpretationofthetransientandrecurrentspiral arm features observed in N-body simulations is still in de- ⋆ E-mail:[email protected] bate.Forexample,Minchev et al.(2012)showedforthefirst 2 J. A. S. Hunt et al. time(bystudyingthetimeevolutionofthediscpowerspec- MeasureGalaxy modelling code,primal (Huntet al.2013; trum) that spiral wave modes in N-body simulations can Hunt& Kawata 2013, 2014), and make a good estimation lastforaslongas1Gyr,whichcanjustifytreatingthewave of the patten speed of the bar, using tracer populations of modes as quasi-stationary structure, and the transient and M0III and Red Clump stars with the Gaia selection func- recurrentspiralarmfeaturescanbeexplainedbythesuper- tion, errors and dust extinction. positionofdifferentmodeswithdifferentpatternspeeds(see There exist full mock catalogues of Gaia stars, e.g. the also Roškar et al. 2012;Sellwood & Carlberg 2014). Onthe Gaia Universe Model Snapshot (gums) which provides a other hand, Grand et al. (2012a); D’Onghia et al. (2013); view of the Besançon Galaxy model as seen from Gaia Baba et al. (2013) demonstrated non-linear growth of the (Robin et al. 2012), taking into account dust extinction spiral arm features due to similar but different (in terms while assuming there are no observational errors. This de- ofevolution)mechanismsfromswing-amplification(Toomre tailed prediction of Gaia observations gives an excellent in- 1981),whichcouldbedifficulttoexplainwiththelinearsu- dication of the volume and quality of data which will be- perposition of thewave modes. come available from Gaia, predicting 1.1 billion observable Our position within the Milky Way gives us a unique stars, almost 10,000 times more than from its predecessor view of these spiral structures seen in external galaxies, Hipparcos. gums can be extended through the Gaia Ob- but it comes with its own set of problems which we must ject Generator (gog) (Luriet al. 2014) to simulate inter- overcome when studying them. The location and kinemat- mediateandfinalcatalogue dataincludingtheintroduction ics of the gaseous component of the arms may be deter- of realistic astrometric, photometric and spectroscopic ob- minedfromHIandCOobservations(e.g.Dame et al.2001; servationalerrorstothecatalogue baseduponGaiascience Nakanishi& Sofue 2003; Kalberla & Kerp 2009). However performance estimates. While these mock data provide an to observe the kinematics of the stellar component in and excellent example of the capabilities of Gaia, the Besançon aroundthespiralarmswemustlookthroughthediscplane, galaxy model is an axisymmetric model and a kinematic which carries the heaviest levels of dust and gas, and thus modelnot a dynamicalmodel. AlthoughGaia will not pro- high levels of extinction. vide accelerations, the kinematics it will provide are from DustextinctionhaslongbeenaproblemforMilkyWay a dynamical system, the Milky Way. Thus it is important model construction. Although there are reasonably reliable for our purpose to generate catalogues from fully dynami- extinction maps for extra galactic sources whose extinc- calmodelswithnon-axisymmetricstructures,suchasspiral tion by the interstellar medium of the Milky Way can be armsandabar,whichforexampleN-bodydiscgalaxymod- corrected as a function A (l,b) (e.g. Schlegel et al. 1998), els can provide. λ threedimensionalextinctionmappingforsourceswithinthe Therefore we propose here to create mock Gaia obser- Milky Way i.e. A (l,b,d) is more challenging. There are vationsfromanN-bodymodelusingapopulationsynthesis λ threedimensionalextinctionmapsforindividualsectionsof codesuchasgalaxia (Sharma et al.2011),orthemethod- thesky(e.g.Drimmel & Spergel2001;Marshall et al.2006; ology presented in Pasetto et al. (2012) or Lowing et al. Hanson & Bailer-Jones 2014; Sale & Magorrian 2014) and (2014). galaxia is a flexible population synthesis code for two dimensional maps have been extended to three dimen- generating a synthetic stellar catalogue from an N-body sions (e.g. Drimmel et al. 2003). However a truly Galactic or an analytical galaxy model over wide sections of the 3D extinction map does not yet exist (Rix & Bovy 2013). sky, with a sampling scheme which generates a smoothly The European Space Agency (ESA)’s Gaia mission will distributed sample of stars. Synthetic catalogues generated help us map the stellar structure and kinematics of the from dynamical Galaxy models are essential for preparing Milky Way, and help constrain extinction at the same time toexploittherealGaiacatalogueandcanbeusedtodeter- (Bailer-Jones et al. 2013). minewhethercertainfeatureswithintheMilkyWaywillbe Gaia, which was launched on the 19th December visible to Gaia. 2013 will providedetailed astrometric (e.g. Lindegren et al. In our previous work (Kawata et al. 2014b) we exam- 2012), spectroscopic (e.g. Katz et al. 2011) and photomet- ined the kinematics of both thestellar and gas components ric (e.g. Jordi et al. 2010) information for around one bil- around a transient, co-rotating spiral arm in a simulated lion stars in the Milky Way. Detailed information on Gaia barred spiral galaxy similar in size to the Milky Way. Al- scientific accuracies is available in, for example, de Bruijne though this arm is transient, similar arms recur during the (2012). Synthetic Gaia mock data have already been used evolution of the galaxy. We made predictions of observable to demonstrate different applications of the real Gaia data kinematicsignaturesthatmaybevisibleintheMilkyWay’s set.Forexample,Abediet al.(2014)usethreetracerpopu- Perseusarmifitisalsoatransient,recurrentandco-rotating lations (OB, A and Red Clump stars) with the Gaia se- spiralarm.Wethencomparedoursimulationwithdatafrom lection function, errors and dust extinction, and demon- APOGEE and the maser sources from Reid et al. (2014) stratedthattheGaiamockdatacanrecovertheparameters measuredbytheBarandSpiralStructureLegacy(BeSSeL) of the Galactic warp. Romero-Gomez et al. (2014) examine survey and the Japanese VLBI Exploration of Radio As- the Galactic bar in the Gaia observable space using Red tronomy (VERA),findingtentativeagreement between our Clump tracers with the Gaia selection function, errors and simulation and the observations. Owing to the low number dust extinction combined with selected Red Clump stars ofmasersourcesandthelackofdistanceinformationforthe from theApachePoint ObservatoryGalactic Evolution Ex- APOGEEstarsnofirmconclusionscouldbedrawn;however periment (APOGEE DR10, e.g. Ahn et al. 2014) showing it is encouraging to see similar features in both, including the value of combining data from complimentary surveys. thepossible signatures of a co-rotating spiral arm. In Hunt& Kawata (2014) we show that we can recover the In this paper we build upon the previous work by gen- large scale structureof theGalactic discwith ourMade-to- erating a stellar sample with different populations from the Spiral arm kinematics: Gaia Mock Observations 3 simulation data in Kawata et al. (2014b) and makingmock observationsofthesestarstakingintoaccounttheexpected Gaiascienceperformanceestimates.Theaimisnottomake further predictions about the kinematics of transient, re- current and co-rotating spiral arms but rather to examine whetherthesesignatures,remain visibleintheGaiadataif they exist in theMilky Way. 2 SIMULATION Figure 1. Snapshot of the simulated galaxy from Kawataetal. We use the simulated galaxy which is presented in (2014b) which is also used in this paper. The left (right) panel Kawata et al.(2014b)andGrand et al.(2014b).Thedetails showstheface-onviewofthestar(gas)particledistribution.The of thenumericalsimulation code, and thegalaxy model are solidlineindicates the position of the spiral arm identified. The described in Kawata et al. (2014b). We briefly describe the observer is assumedto be located at (x,y)=(−8,0) kpc. Three galaxymodelinthissection.Thegalaxyissetupinisolated line-of-sight directions (lLOS = 90,120 and 150 deg) are high- conditions,andconsistsofagasandstellardiscbutnobulge lightedwiththedotted lines.Thegalaxyisrotatingclockwise. component. The discs are embedded in a static dark mat- ter halo potential (Rahimi & Kawata 2012; Kawata et al. 2014b). The dark matter halo mass is M = 2.5×1012 dm M⊙, and the dark matter density follows the density pro- filefromNavarroet al.(1997),withaconcentrationparam- tionsaremadeoftheobservationalsignaturesofco-rotating eter of c = 10. The stellar disc is assumed to follow an spiralarmsnotablythedifferenceinkinematicstructurebe- exponential surface density profile with the initial mass of tweenthetrailingnearsideandleadingfarsideofthespiral Md,∗ = 4.0×1010 M⊙, a radial scale length of Rd,∗ = 2.5 arm. In general, in Kawata et al. (2014b) (as also shown in kpc and a scale height of zd,∗ =350 pc. The gas disc is set Grand et al. (2014a)) the stars in the trailing near side ro- up following the method of Springelet al. (2005), and has tate slower because they tend to be at the apo-centre and an exponential surface density profile with the scale length migrate outward, and the stars in the leading far side ro- of Rd,g = 8.0 kpc. The total mass of the gas is 1010 M⊙. tatefasterastheytendtobeattheperi-centreandmigrate Thesimulation comprises 106 gas particlesand 4×106 star inward. There are however some stars which follow the op- particles;thereforeeachparticlehasamassof104 M⊙.The positetrend,leadingtomultiplepopulationsseeninthero- resolution is sufficient to minimise numerical heating from tational velocity in the leading far side; one faster, and one Poisson noise (Fujii et al. 2011;Sellwood 2013).Weemploy slower than the single population in the trailing near side. aminimumsofteninglengthof158pc(equivalenttoaPlum- These features which will be discussed later may be caused mersofteninglengthof53pc)withthesplinesofteningand by the co-rotation resonance of the spiral arm, and are vis- variable softening length for gas particles as suggested by ible at different galactic longitudes because the arm in the Price & Monaghan (2007). simulationco-rotatesatalltheexaminedradialrange.How- The radial profile of the mean metallicity of stars and ever,inKawata et al.(2014b),thespiralarmkinematicsare gas is initially set by[Fe/H](R)=0.2−0.05(R/1 kpc), and examined using the full, error and extinction free N-body themetallicitydistributionfunctionateachradiusiscentred data and thus such trends, when present, are easy to iden- on the mean metallicity value with the dispersion set to a tify. Gaussiandistributionof0.05dexforthegasand0.2dexfor In this Section we describe how we generate a sample the stars. The stellar ages are set randomly between 0 and ofstarsfromtheN-bodymodelofKawata et al.(2014b)to 10 Gyrfor stars present at thebeginning of thesimulation. produceamockGaiacatalogue. Itisworthnotingthatthe The simulation was run for 1 Gyr from the ini- population synthesis code, galaxia (Sharmaet al. 2011) tial conditions with the N-body smoothed particle hy- providesatooltogeneratestellarpopulationsfromN-body drodynamics code, gcd+ (e.g. Kawata & Gibson 2003; simulationdata.However,becauseweplantocombinesuch Rahimi & Kawata 2012; Barnes et al. 2012; Kawata et al. a tool with our Made-to-Measure Galaxy modelling code, 2013, 2014a) without the inclusion of any continuous ex- primal, we have developed our own simplified version of ternal inflow of gas for simplicity. In this paper we use the galaxia,apopulationsynthesiscodecalledsnapdragons, samesnapshotofthegalaxyasusedinKawata et al.(2014b) (Stellar Numbers And Parameters Determined Routinely which is taken at t = 0.925 Gyr, as this snapshot shows a AndGeneratedObservingN-bodySystems).snapdragons spiralarmatasimilarlocation tothePerseusarmfromthe uses the same isochrones and extinction map as galaxia, Milky Way(see Fig. 1). but uses a different and more simplistic process to gener- ate the stellar catalogue which is described in Section 3.2. snapdragons allows us to add the expected Gaia errors more easily, and enables us to track the link between sam- 3 GAIA MOCK CATALOGUE pled stars and their parent N-body particle for our future In Kawata et al. (2014b) the kinematics of the spiral arm studies,e.g.primalmodellingoftheGalacticdiscbyfitting shown in Fig. 1 are examined at three lines of sight l = tracersfrommultiplestellarpopulations,andidentifyingra- LOS 90,120 and 150 deg, with blos = 0 because of the lower ex- dially migrating stars and non-migrating stars trapped by tinction relative to other lines of sight in theplane. Predic- thespiral arm (Grand et al. 2014a). 4 J. A. S. Hunt et al. 3.1 Extinction age from the grid of isochrones which are extracted from We use the extinction map of the Milky Way taken from galaxia.Onceanisochroneisselected,weidentifym⋆,i,max from the isochrone. We then determine how many stars to galaxia (Sharma et al. 2011), which is a 3D polar loga- sample from the N-body particle by integrating the IMF rithmic grid of the dust extinction constructed using the overthedesired mass range; method presented in Bland-Hawthorn et al. (2010) and the dust maps from Schlegel et al. (1998). The same extinc- N =A m⋆,i,max m−(x+1)dm, (4) tion is applied in Hunt & Kawata (2014) and more detail s Z m⋆,i,<Vlim is given there. In an update from Hunt & Kawata (2014) we follow the correction to the Schlegel EB−V presented in wherem⋆,i,<Vlim isminimummassrequiredforthestarpar- Sharma et al. (2014) such that ticle to be brighter than our apparent magnitude selection limit, Vlim, taking into account the extinction value at the EB−V =EB−V(cid:18)0.6+0.2(cid:18)1−tanh(cid:18)EB−V0.−1 0.15(cid:19)(cid:19)(cid:19). positionoftheparentparticle.Starssmallerthanm⋆,i,<Vlim are not used in the subsequent analysis, to save on compu- (1) tational time. This correction is made as it has been suggested (e.g. We then randomly sample stellar masses from the sec- Arce& Goodman 1999; Yasudaet al. 2007) that the red- tion of the isochrone N times. We have weighted the ran- dening is overestimated by the maps from Schlegel et al. s dom selection by theIMF using theequation (1998) by ∼1.3-1.5 in regions with high extinction with AtioVn>by0∼.5 4(E0%B−fVor>low0.1la5t)i.tuTdheishcigohrreexcttiionnctiroendurceegsioenxstibnuct- m⋆ =(Rm−⋆,xi,max+(1−R)m−⋆,xi,<Vlim)−1x, (5) has minimal effect on high latitude low extinction regions. where R is a random numberbetween 0 and 1. The isochrones are comprised of discrete stellar data, and therefore we then interpolate within the nearest 3.2 Population Synthesis: snapdragons isochrone values of M and V −I to determine M and V c V⋆ ThegoalofthispopulationsynthesiscodeistospliteachN- V −Ic⋆ for the generated m⋆. At this stage we assume the bodyparticlefromthegalaxysimulationintoanappropriate generatedstarshavethesamepositionandvelocityastheir numberof stellar particlescreating amockcatalogue of ob- parent particles. servable stars from our N-body model. We must choose an IMFandaset ofisochrones withwhichtowork.Wechoose a Salpeter IMF (Salpeter 1955) where the IMF, Φ(m), is 3.3 Observational Errors definedin each mass interval dm as Having generated the visible stellar catalogue we then add Φ(m)dm=Am−(x+1)dm, (2) observational errors based upon the Gaia Science Perfor- mance estimates1. We use the post launch error estimates where x = 1.35 is the Salpeter index, and A is a constant approximatedfromtheestimatesinpre-launchperformance for normalisation in the desired mass range. We set this byMercè Romero-Gómez (e.g. Romero-Gomez et al. 2014), constant as providedthroughtheGaiaChallengecollaboration2.Weas- A =m m⋆,i,maxm−xdm −1, (3) sume the position and velocity of the Sun is known. We i i (cid:18)Zm⋆,min (cid:19) locate the observer at (−8,0,0) kpc as shown in Fig. 1, and themotionoftheSunisassumedtobe228kms−1.Forthis wheremi istheN-bodyparticlemass,m⋆,i,max isthemax- work, while generating the stellar catalogue we produced imum initial mass of any surviving star and m⋆,min is the minimum stellar mass to be considered. We make use of starsonlybrighterthanVlim≤16mag,whichiswellwithin Gaia′s m ≤ 20 mag magnitude limit for the astrometry. thePadovaisochrones(e.g.Bertelli et al.1994;Marigo et al. G However, because we are interested in the Galactic radial 2008), although thechoice of isochrones (and IMF) may be and rotation velocity for the stars, which requires the full substitutedwithotherswithnochangetothemethodology. 6D phase space information, we chose the lower magnitude ItisworthnotingthatthePadovaisochronesareavail- able only for stellar masses above 0.15 M⊙. galaxia for limit whereGaiaRVScanproducethereasonably accurate line-of-sight velocity. Note that the errors are added to the example uses the isochrones from Chabrier et al. (2000) to parallax, propermotion and line-of-sight velocities. extend the mass limit down to 0.07 M⊙, which is the hy- A full description of the method to add the pre-launch drogenmassburninglimit.Wesetourlowerlimitonstellar Gaia error is available in Hunt& Kawata (2014). However mass as m⋆,min = 0.1 M⊙ to correspond with the simula- the Gaia science performance estimates have been revised tion from Kawata et al. (2014b) and extrapolate from the Padova isochrones for 0.1≤M⊙ ≤0.15. It is relatively safe after launch, and as such a correction must be made. The error in parallax has increased, and although it has little to do this because all such stars lie on the main sequence. effectforstarswithm ≤16magwhichweworkwithinthis Additionallytheseexceedinglyfaint starswillnotbevisible V paper, the coefficients within the equation to describe the atthedistanceofthespiralarmswhicharethefocusofthis pre-launchparallax performance(providedbyKazi, Antoja work. & DeBruijne (Oct. 2014) by fitting to the new estimations As discussed in Section 2 each N-body star particle in thesimulatedgalaxyhasbeenassignedanageandmetallic- itywithinthechemodynamicalcodegcd+,thenit ismade to evolve. When we examine the snapshot, each particle is 1 http://www.cosmos.esa.int/web/Gaia/science-performance matched to its nearest isochrone in both metallicity and 2 http://astrowiki.ph.surrey.ac.uk/dokuwiki/doku.php Spiral arm kinematics: Gaia Mock Observations 5 on theGaia science performance web page) are revised to σ = (−11.5+706.1z+32.6z2)1/2 π ×(0.986+(1−0.986)(V −I )), (6) c where z=max(100.4(12−15),100.4(G−15)), (7) correcting also the typo for equations (6) and (7) in Hunt& Kawata (2014). Additionally, because of the loss of spectroscopic accu- racy by ∼1.5 mag in the RVS post launch performance we also apply a correction to the error function for the end of missionradialvelocity.Wechangethetable3 ofvaluesaand b, again determined by fitting the revised performance es- timates on the Gaia science performance web page, for the equation σ =1+bea(V−14), (8) vr whereaandbareconstantsdependantonthespectraltype of the star. The new table along with the code to add the Gaia error is available online4. Figure 2. Colour magnitude diagram for stars generated by snapdragons from particles within a square region of ±2 deg (upper) and ±5 deg (lower) around (l,b) = (90,0) deg. Stars withapparent magnitudeofmV ≤16onlyareincluded. 4 RESULTS As discussed in Section 3, it was shown in Kawata et al. many particles to cover a broad range of stellar ages and (2014b) that in general thestars in thetrailing near side of metallicities in the CMD. Therefore, care is required with thespiralarmrotateslowerthanaveragebecausetheytend the resolution of the N-body simulation and the selection tobeattheapo-centre,andthestarsintheleadingfarside function if we discuss in detail thestellar population distri- of the spiral arm rotate faster than average as they tend bution in the CMD. However, this is unlikely to affect the to be at the peri-centre. However, there are groups of stars studyin this paper. whichfollowdifferenttrendsleadingtomultiplepopulations which will be discussed later. It is important to determine whether such features will still be visible in the Gaia cata- 4.2 Observable Spiral Arm Kinematics logue, notjust theerrorand extinction-freeN-body model. In this section we examine if the possible kinematic signa- InthisSectionweshowtheresultofsamplingtheseN-body turesofco-rotatingtransientandrecurrentspiralarmsiden- data into stellar data, first looking at the properties of the tified in Kawata et al. (2014b) will be visible in the Gaia resultingmockstellarcatalogue,andthenexaminingthespi- data even given the dust extinction in the disc and Gaia′s ralarmkinematicswiththestellardatatakingintoaccount science performance accuracy. A detailed analysis of the dust extinction and Gaia science performance estimates. kinematicsthemselvesisthefocusofKawata et al.(2014b), while thiswork isconcerned with thevisibility of thiskine- 4.1 Population synthesis maticstructureintheGaiadata.Weexaminetherotational velocitiesofthestarsinthecataloguefordifferentdistances In this section we describe the stellar catalogue produced because in Kawata et al. (2014b) we found the rotation ve- bysnapdragons,andshowtheresultingcolourmagnitude locity is most affected by the transient co-rotating spiral diagram(CMD)varyingtheareaoftheskycoverage.Fig.2 arm. Then we calculated the Probability Density Function shows the CMD for stars generated by snapdragons from (PDF) of the rotation velocity of stars behind and in front particles within a square region of ±2 deg (upper) and ±5 of the spiral arm using Kernel Density Estimation (KDE) deg (lower) around (l,b) = (90,0) deg. The upper panel which we are using as a desirable alternative to histograms of Fig. 2 shows clearly the individual stellar isochrones be- (e.g. Wasserman 2006). cause there are only a small numberof N-body particles in Fig. 3 shows a smoothed contour plot of the galacto- the selected region, and each particle has only one age and centricrotational velocity against distanceforparticles and metallicity. These problems are resolved when smoothing is stars within a square region of ±5 degrees around (l,b) = applied in the phase space distribution and age-metallicity (90,0) (left), (l,b) = (120,0) (middle) and (l,b) = (150,0) distribution(e.g.Sharma et al.2011).However,asdiscussed (right). This compares thekinematics of theunderlyingN- inSection3.2wedeliberatelyavoidthissmoothingtomain- body model (upper) with the stellar catalogue generated tain the clear particle-star relation. The lower panel of Fig. with snapdragons, before (middle) and after (lower) the 2 shows no such discrete structure, as there are sufficiently addition of the errors from the Gaia science performance estimates.Owingtothehighpercentageoflowmassandlu- 3 http://www.cosmos.esa.int/web/Gaia/table-5 minositystellartypeswhichwoulddominatetheselectedre- 4 https://github.com/mromerog/Gaia-errors gionandsaturatetheplotatsmalldistances,wehavemade 6 J. A. S. Hunt et al. cuts to our sample to visualise the underlying kinematic inKawata et al.(2014b).Thisisapositiveoutcomeconsid- structure from the stellar catalogue. We have first cut the ering the loss of data from the dust extinction. When com- sampleofstarsinallthreelinesofsightwithabsolutemag- paringthe‘true’(blacksolid)stellarcataloguedatawiththe nitude, M ≤−1, calculated from the apparent magnitude stellar data taking into account dust extinction and Gaia’s V mV and observed distance dobs, assuming the dust extinc- expectederrors(reddashed)ageneralsmoothingoutofthe tion is known. We then cut with σvlos/(vlos×dobs)≤0.015 structureisevidentinthe‘observed’data.Theupperpanels kpc−1 toselectthestarswithlowererrorinthelineofsight ofFig.4showingthetrailingnearsideofthearmshowvery velocitiesatasmallerdistancetogeneratesimilarquantities similar PDF’s when comparing the true and observed stel- of data at different distance scales. This is purely for illus- lar data, whereas the lower panels showing the leading far tration purposes and we are not suggesting that this is the side show an information loss, especially at (l,b) = (90,0), best possible selection function. The upper panels of Fig. 3 where the three peaks are no longer resolved. This is to be showthedifferentkinematicstructureintheN-bodymodel expected because of the higher distances and therefore ad- at the different lines of sight. These are the same data as ditional extinction; however at (l,b)= (120,0) and (150,0) those shown in the top panels of Fig. 4 from Kawata et al. even on the far side of the spiral arm the structure within (2014b). Note that the density colour scale for the N-body thedistribution is still clearly visible. dataisdifferentfromthestellardatainthemiddleandlower When comparing the ‘observed’ data in Fig. 4 in front panels. and behind the spiral arms, we see a clear difference in the ThemiddlerowofpanelsofFig.3showthevelocitiesof PDF at all three lines of sight. In each case, the PDF in theselectedstars,whichappearslightlydifferentfromthose the trailing near side of the spiral arm forms a single cen- of the N-body data owing to the selection function. While tralpeaksimilartothemeanrotationvelocity,withasmall the general shape of the distribution has been recovered, tail towards faster rotation velocities whereas the leading at (l,b) = (90,0) deg (middle left) the fast rotating stars far side of the spiral arm shows a broader distribution of withinthearmdominatethedensityscaleandwashoutthe velocities with a peak velocity faster than the peak for the rest of the plot slightly. At (l,b) = (120,0) deg (middle), although there is some saturation around 220 km s−1 the trailing nearside. Thedifferenceis particularly apparentat (l,b) = (120,0) deg where the leading far side shows two kinematicstructureisclearlyvisibleandisagoodmatchto clear peaks, one faster and one slower than the single peak the particle data. Similarly at (l,b) = (150,0) deg (middle in thetrailing near side.This bimodal distribution can also right), despite the lower number of counts, the kinematic beseeninthelowermiddlepanelofFig.3between4.39and structureis clearly shown. 5.39kpc(althoughnotethatFig.3usesadifferentselection ThelowerpanelsofFig.3showtheerroraffectedrota- function).Alsoat(l,b)=(150,0)degthesinglebroadpeak tionvelocityanddistancefortheselectedstarstakingGaia in the trailing near side is easily distinguishable from the scienceperformanceestimatesintoaccount.Therotationve- leadingfarsidewhichshowsthreepeaks.Thesethreepeaks locity is calculated from the observed parallax, proper mo- are also partially visible in the lower right panel of Fig. 3 tion and line of sight velocities. At (l,b) = (90,0) (lower between 4.29 and 5.29 kpc. These features all match those left) the shape of the distribution remains relatively un- observed in Kawata et al. (2014b) despite the addition of changed, with the main loss in accuracy occurring around dust extinction and observational errors to thedata. dobs ≈ 7− 10 kpc. The recovery of the kinematic struc- ture around the spiral arm around dobs ≈ 4 kpc remains In general, as shown in Grand et al. (2014a), the stars almostidenticaltothecasewithoutobservationalerrors.At in the leading side rotate faster as they tend to be at peri- (l,b) = (120,0) (lower middle) the visible loss of accuracy centrephaseandmigratinginward,andstarsinthetrailing is again in the outer region of dobs ≈ 7−10 kpc, with the side rotate slower as they tend to be at apo-centre phase region containing the spiral arm remaining very similar to and migrating outward. This explains the single large peak that of the error free case. At (l,b)=(150,0) (lower right), inthetrailingside,andthelargest peakontheleadingside the entire distribution remains very similar to the middle which has a higher rotational velocity than the single peak rightpanel,thecasewithoutGaialikeobservationalerrors. on the trailing side as shown in Fig. 4. However, when the Fig.4showsthePDF’s,withaKDEbandwidthof4,for transient spiral arm starts forming, stars which are close to the rotational velocity of the stars in the catalogue within thearm on thetrailing side andare close to theperi-centre a square region of ±5 degrees around (l,b) = (90,0) (left), phase,areacceleratedtowardsthearm,passingthroughand (l,b) = (120,0) (middle) and (l,b) = (150,0) (right) in the thenslowingdownastheyreachtheapo-centreonthelead- trailing nearside, between 1 and 2 kpccloser than thecen- ing side as discussed in Kawata et al. (2014b). These stars tre of the arm (upper) and leading far side between 1 and correspond to the ‘slower’ peaks visible in the lower panels 2 kpcfurtherthan the centreof the arm (lower). Note that ofFig.4.Similarly,thestarswhichareclosetothearmand these distance bins were chosen as they show the discussed close to the apo-centre phase on the leading side are decel- structure most clearly; the same features are present closer eratedbythearm,andareovertakenbythearm.Thenthey tothearmbutarelessclear.Thecentreofthearmwasde- areaccelerated again bythearm oncetheyareonthetrail- terminedtobeatd=4.0kpcat(l,b)=(90,0),d=3.4kpc ingsideatperi-centrephase,whichcorrespondstothesmall at (l,b)=(120,0) and d=3.3 kpc at (l,b)=(150,0). Note tail present at high velocities in the upper panels of Fig. 4. that Fig. 4 uses all the stars with m ≤ 16, not applying The difference in the rotation velocity distribution between V the selection function used for illustration purposes in Fig. the leading and trailing side of the spiral arm seen in Figs. 3.AtallthreelinesofsightFig.4showsacleardifferencein 3 and 4 is that the latter population is smaller than the thedistributionofvelocitiesforthe‘true’data(blacksolid) former. It appears that it is easier for stars to escape from when comparing the different observed distances, as shown thearmontheleadingsidethanthetrailingside.Fromour Spiral arm kinematics: Gaia Mock Observations 7 Figure 3.Smoothed linearscalecontourplotofgalactocentric rotationvelocityofsimulationparticles(upper), selectedsnapdragons stars (middle)andselected snapdragons stars observed withGaia error (lower) for(l,b)=(90,0)(left), (l,b)=(120,0) (middle)and (l,b)=(150,0)(right).Forthesnapdragonsstars(middleandlowerpanels),alimitedselectionofMV ≤−1calculatedusingmV and dobs andassumingaknownextinction, alongwithσvr/(vr×dobs)≤0.15isshowntoavoidoverlydensepopulations offainterstars at smallerdistances.Thisistovisualisethedataset,andthesefaintstarscontribute tothesubsequentanalysis.Notewhileconsistentfor snapdragonsstarsacrossthedifferentlinesofsight,thedensityscaleisdifferentforthesimulationparticles;howeverthechoiceofthe scaleisarbitrary. analysisofN-bodysimulationsthisappearstobeacommon thekinematicsintheleadingandtrailing sidesofthespiral feature of transient and co-rotating spiral arms. arm, notably the difference in the number and locations of Comparetta & Quillen (2012) propose that the radial the peaks, and the small high velocity tail present in the overlapofmultiplelonger-livedpatternsmovingatdifferent trailing near side. The comparison between the middle and pattern speeds can reproduce the transient spiral features, lower panels of Fig. 3 shows little difference, implying that whichwhenstrongenoughcanleadtoradialmigrationaway theobservationalerrorfromGaiawillhavelimitedeffecton from theco-rotation radius associated with co-rotating spi- our ability to study the kinematics of the spiral arms. Fur- ral arms as seen, for example in Grand et al. (2012a,b). ther examination of galaxy models constructed using the In such a scenario, the spiral arm features are co-rotating, differenttheoriesofspiralarmformationwillbeessentialto which may give rise to the co-existence of many inner and determinethe distinct kinematic signatures of each theory. outerLindbladresonancesinarangeofradiiandleadtothe features visible in Figs. 3 and 4. However, further analysis of the spiral arms in N-body simulations is required before 5 SUMMARY drawing firm conclusions on the mechanism that generates such kinematic signatures, which we will tackle in future We observed our N-body/SPH simulation of a Milky Way studies. like barred spiral galaxy to create a mock Gaia stellar cat- From Figs. 3 and 4 we find thatGaia′s scientific accu- alogue, with particular interest in the stellar kinematics in racy ought to be sufficient to examine the kinematic struc- and around the spiral arms. We focused on the same three ture of the nearby spiral arms in the Milky Way, even on lines of sight in the disc plane as Kawata et al. (2014b), the far side of the arm. Fig. 4 shows clear differences in (l,b) = (90,0),(120,0) and (150,0) deg and analysed the 8 J. A. S. Hunt et al. 0.025 Trailing near side (l,b)=(90,0) (l,b)=(120,0) (l,b)=(150,0) 0.020 0.015 0.010 0.005 F D0.000 P000000......000000012012000555 Leading1 f5a0r side 200 250 300 150 200 250 300 150 200 250 300 Vrot (km s−1) Figure4.Comparisonofthedistributionofgalactocentricrotationalvelocitiesforthestarsgeneratedbysnapdragonswithinasquare region of ±5degrees around (l,b)=(90,0) (left), (l,b)=(120,0) (middle) and (l,b)=(150,0) (right) inthe trailing near side(upper) and leading farside (lower) of the spiral arm which meet the mV ≤16selection limit.The black solidcurve shows the true velocities, and the red dashed curve shows the distribution once the Gaia errors have been applied. The vertical lines show the circular velocity (dotted) andthemeanrotation velocity(dash-dotted) attheradiusofthespiralarm. galactocentric rotational and line of sight velocities of the ACKNOWLEDGEMENTS selected stars as a function of the distance from the ob- We gratefully acknowledge the support of the UK’s server.InagreementwithexistingliteratureonN-bodyspi- Science & Technology Facilities Council (STFC Grant ral galaxy simulations the spiral arm features seen in the ST/H00260X/1 and ST/J500914/1). The calculations for stellar mass in our model are transient, recurrent and co- thispaperwereperformedonCrayXT4atCenterforCom- rotating,i.e.thespiralarmisrotatingatthecircularveloc- putational Astrophysics,CfCA, of theNational Astronomi- ity of thestars at theselected lines of sight. calObservatoryofJapanandtheDiRACfacilities(through the COSMOS consortium) including the COSMOS Shared We show that the structure in the kinematics identi- Memory system at DAMTP, University of Cambridge op- fiedinKawata et al.(2014b)remainsvisibleaftertheinclu- erated on behalf of the STFC DiRAC HPC Facility. This sion of dust extinction and observational errors based upon equipment is funded by BIS National E-infrastructure cap- Gaiascienceperformanceestimates.Althoughtheinclusion ital grant ST/J005673/1 and STFC grants ST/H008586/1, of these observational effects makes the trends less clear, ST/K00333X/1&ST/J001341/1. Theauthorsacknowledge they are still observable in the mock Gaia data in front of, the use of the IRIDIS High Performance Computing Fa- insideandbehindthespiralarm.Thestructureonthetrail- cility, and associated support services at the University of ing nearside is relatively unchanged,whereas thestructure Southampton. We would also like to thank PRACE for the on the leading far side is, unsurprisingly, more affected, al- use of the Cartesius facility. This work was carried out, though the bi-modal (or more) and broader distribution of in part, through the Gaia Research for European Astron- therotation velocities is stillclearly visible. Because webe- omy Training (GREAT-ITN) network. The research lead- lievethatthesekinematicsignaturesareindicationsoftran- ingtotheseresultshasreceivedfundingfrom theEuropean sient and co-rotating spiral arms owing to the co-rotation Union Seventh Framework Programme ([FP7/2007-2013] resonance at all radii, we predict they should be visible in undergrant agreement number264895). Wewould also like the Gaia data at different longitudes if the Milky Way’s to thank Mercè Romero-Gómez and Francesca Figueras for Perseus arm is also a transient and co-rotating spiral arm. providingthesubroutinetocalculatetheGaiaperformance errors, including the update to post launch estimates, San- jib Sharma for providingthegalaxia extinction maps and Encouraged by the success of this study, we intend isochronesandSamiNiemiforthesuggestionofusingKDE’s to repeat the analysis with simulated galaxies which use to visualise the velocity distributions. different theories of spiral structure formation, for exam- ple test particle simulations (e.g. Minchev & Quillen 2008; Minchev et al. 2010; Minchev & Famaey 2010; Faure et al. References 2014; Antoja et al. 2014) and N-body simulations with a fixed spiral arm potential (e.g. Wadaet al. 2011). From Abedi H., Mateu C., Aguilar L. A., Figueras F., Romero- these analyses we expect to make predictions of the kine- Gómez M., 2014, MNRAS,442, 3627 matic signatures of different spiral arm theories, which can Ahn C. P. et al., 2014, ApJS,211, 17 betested by theGaia stellar catalogue. Antoja T. et al., 2014, A&A,563, A60 Spiral arm kinematics: Gaia Mock Observations 9 Arce H.G., Goodman A. A., 1999, ApJ, 512, L135 Minchev I., Boily C., Siebert A., Bienayme O., 2010, MN- Baba J., Saitoh T. R.,Wada K., 2013, ApJ, 763, 46 RAS,407, 2122 Bailer-Jones C. A. L. et al., 2013, A&A,559, A74 Minchev I., Famaey B., 2010, ApJ, 722, 112 BarnesD.J.,KawataD.,WuK.,2012,MNRAS,420,3195 Minchev I., Famaey B., Quillen A. C., Di Matteo P., Bertelli G., Bressan A., Chiosi C., Fagotto F., Nasi E., CombesF.,VlajićM.,ErwinP.,Bland-HawthornJ.,2012, 1994, A&AS,106, 275 A&A,548, A126 Bland-Hawthorn J., Krumholz M. R., Freeman K., 2010, Minchev I., Quillen A. C., 2008, MNRAS,386, 1579 ApJ,713, 166 Nakanishi H.,Sofue Y., 2003, PASJ,55, 191 Bottema R., 2003, MNRAS,344, 358 NavarroJ.F.,FrenkC.S.,WhiteS.D.M.,1997,ApJ,490, Carlberg R. G., Freedman W. L., 1985, ApJ, 298, 486 493 ChabrierG.,BaraffeI.,AllardF.,HauschildtP.,2000,ApJ, Pasetto S., Chiosi C., Kawata D., 2012, ArXiv e-prints: 542, 464 1207.5048 Comparetta J., Quillen A.C., 2012, ArXiv e-prints Price D.J., Monaghan J. J., 2007, MNRAS,374, 1347 Dame T. M., Hartmann D., Thaddeus P., 2001, ApJ, 547, Quillen A. C., Dougherty J., Bagley M. B., Minchev I., 792 Comparetta J., 2011, MNRAS,417, 762 de Bruijne J. H. J., 2012, Ap&SS,341, 31 Rahimi A., Kawata D., 2012, MNRAS,422, 2609 Dobbs C., Baba J., 2014, PASA,31, 35 Reid M. J. et al., 2014, ApJ, 783, 130 D’Onghia E., Vogelsberger M., Hernquist L., 2013, ApJ, Rix H.-W., Bovy J., 2013, A&A Rev.,21, 61 766, 34 Robin A. C. et al., 2012, A&A,543, A100 Drimmel R., Cabrera-Lavers A., López-Corredoira M., Roca-Fàbrega S., Valenzuela O., Figueras F., Romero- 2003, A&A,409, 205 Gómez M., Velázquez H., Antoja T., Pichardo B., 2013, Drimmel R., Spergel D. N., 2001, ApJ, 556, 181 MNRAS,432, 2878 FaureC.,SiebertA.,FamaeyB.,2014,MNRAS,440,2564 Romero-Gomez M., Figueras F., Antoja T., Abedi H., Fujii M. S., Baba J., Saitoh T. R., Makino J., Kokubo E., Aguilar L., 2014, ArXiv e-prints WadaK., 2011, ApJ,730, 109 Roškar R., Debattista V. P., Loebman S. R., 2013, MN- Grand R. J. J., Kawata D., Cropper M., 2012a, MNRAS, RAS,433, 976 421, 1529 RoškarR.,DebattistaV.P.,QuinnT.R.,WadsleyJ.,2012, Grand R. J. J., Kawata D., Cropper M., 2012b, MNRAS, MNRAS,426, 2089 426, 167 Sale S. E., Magorrian J., 2014, MNRAS,445, 256 Grand R. J. J., Kawata D., Cropper M., 2013, A&A, 553, Salpeter E. E., 1955, ApJ, 121, 161 A77 SchlegelD.J.,FinkbeinerD.P.,DavisM.,1998,ApJ,500, Grand R. J. J., Kawata D., Cropper M., 2014a, MNRAS, 525 439, 623 Sellwood J. A.,2011, MNRAS,410, 1637 Grand R. J. J., Kawata D., Cropper M., 2014b, ArXiv e- Sellwood J. A.,2013, ApJ,769, L24 prints Sellwood J. A.,Carlberg R.G., 1984, ApJ,282, 61 Hanson R. J., Bailer-Jones C. A. L., 2014, MNRAS, 438, Sellwood J. A.,Carlberg R.G., 2014, ApJ,785, 137 2938 Shapley H., Curtis H. D., 1921, Bulletin of the National Hunt J. A.S., Kawata D., 2013, MNRAS,430, 1928 Research Council, Volume2, Part 3, 11, 171 Hunt J. A.S., Kawata D., 2014, MNRAS,443, 2112 Sharma S.et al., 2014, ApJ, 793, 51 Hunt J. A.S., Kawata D., Martel H., 2013, MNRAS,432, SharmaS.,Bland-HawthornJ.,JohnstonK.V.,BinneyJ., 3062 2011, ApJ, 730, 3 Jordi C. et al., 2010, A&A,523, A48 Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, Kalberla P. M. W., Kerp J., 2009, ARA&A,47, 27 361, 776 Katz D. et al., 2011, in EAS Publications Series, Vol. 45, Toomre A., 1981, in Structure and Evolution of Normal EAS Publications Series, pp.189–194 Galaxies, Fall S. M., Lynden-BellD.,eds., pp.111–136 Kawata D., Gibson B. K., 2003, MNRAS,340, 908 Wada K., Baba J., Saitoh T. R., 2011, ApJ, 735, 1 Kawata D., Gibson B. K., Barnes D. J., Grand R. J. J., Wasserman L., 2006, All of Nonparametric Statistics Rahimi A.,2014a, MNRAS,438, 1208 (SpringerTextsin Statistics). Springer-VerlagNewYork, KawataD.,HuntJ.A.S.,GrandR.J.J.,PasettoS.,Crop- Inc.,Secaucus, NJ, USA perM., 2014b, MNRAS,443, 2757 Wilczynski E. J., 1896, ApJ, 4, 97 Kawata D., OkamotoT., Gibson B. K., Barnes D. J., Cen Yasuda N., Fukugita M., Schneider D. P., 2007, AJ, 134, R.,2013, MNRAS,428, 1968 698 Lin C. C., ShuF. H., 1964, ApJ,140, 646 LindegrenL.,LammersU.,HobbsD.,O’MullaneW.,Bas- tian U.,HernándezJ., 2012, A&A,538, A78 Lowing B., Wang W., Cooper A., Kennedy R., Helly J., FrenkC., Cole S., 2014, ArXiv e-prints Luri X.et al., 2014, A&A,566, A119 Marigo P.,GirardiL.,Bressan A.,Groenewegen M.A.T., SilvaL., Granato G. L., 2008, A&A,482, 883 Marshall D. J., Robin A. C., Reylé C., Schultheis M., Pi- caud S., 2006, A&A,453, 635