Table Of Content

Mon.Not.R.Astron.Soc.000,1–13(2012) Printed23October2012 (MNLATEXstylefilev2.2) Disc galaxy modelling with a particle-by-particle M2M method 2 1 Jason A. S. Hunt,1⋆ Daisuke Kawata,1† 0 1 Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey, RH5 6NT, UK 2 t c O SubmittedtoMNRAS,27July,2012 9 1 ABSTRACT ] We have developed the initial version of a new particle-by-particle adaptation of the A made-to-measure (M2M) method, aiming to model the Galactic disc from upcoming G Galactic stellar survey data. In our new particle-by-particle M2M, the observables . of the target system are compared with those of the model galaxy at the position h of the target stars (i.e. particles). The weights of the model particles are changed to p reproducetheobservablesofthetargetsystem,andthegravitationalpotentialisauto- - o maticallyadjustedby the changingweightsofthe particles.This paperdemonstrates, r astheinitialwork,thattheparticle-by-particleM2Mcanrecreateatargetdiscsystem t s created by an N-body simulation in a known dark matter potential, with no error in a the observables. The radial profiles of the surface density, velocity dispersion in the [ radial and perpendicular directions, and the rotational velocity of the target disc are 1 all well reproduced from the initial disc model, whose scale length is different from v that of the target disc. We also demonstrate that our M2M can be applied to an in- 1 complete data set and recreate the target disc reasonably well when the observables 2 arerestrictedtoapartofthedisc.Wediscussourcalibrationofthemodelparameters 5 and the importance of regularization. 5 . Key words: methods: N-body simulations — galaxies: structure — galaxies: kine- 0 1 matics and dynamics — The Galaxy: structure 2 1 : v i 1 INTRODUCTION mission, Gaia. Nano-JASMINEisa demonstration mission, X butitwilllikelyimproveonHipparcos’propermotions.Gaia r There is still a gulf between our theoretical galaxy models isexpectedtolaunchin2013andwilloperateforfiveyears, a andtheobservational data,thatmust bebridgedbefore we withapossibleoneortwoyearextension.Gaiawill provide can have a fully dynamical model of the Milky Way which anunprecedentedlylargeamountofinformationwithwhich is consistent with its observed properties. The major devel- to build a more accurate model of theMilky Way. opmentshavebeen localised tocertain regions of theMilky Constructing accurate models of the Milky Way is im- Way and the structure of many other regions of the Milky portant for allowing us to understand and compensate for Way remains largely uncertain. observational bias, which are present in all existing Galac- For the last two decades Galactic astronomy has tic surveys due to dust, gas and our location within the been relying on Hipparcos data (e.g. Perryman & ESA disc. They also allow us to tie together data from differ- 1997). However, new space-based astrometry missions ent surveys, assembling them into a single model. There are going ahead in the near future, and ground- arethreedifferenttypesofgalaxymodel.Massmodelsonly based surveys, e.g. PanStaars, (e.g. Kaiser et al. 2010), describe the density distribution and the galactic potential VISTA (e.g. Minniti et al. 2009), LSST (e.g. Ivezić et al. (e.g.Klypinet al.2002).Kinematicmodelsspecifytheden- 2008), SEGUE (e.g. Yannyet al. 2009), APOGEE (e.g. sity and velocity distributions, but lack the constraint that AllendePrieto et al.2008)andRAVE(e.g.Steinmetz et al. the model must be in a steady state in the galactic poten- 2006),willaddsignificantvaluetothesemissions,whichwill tial(e.g.Robin et al.2004).Amodelwhichalsosatisfiesthis expand our knowledge of the Milky Way. The next space criterionisknownasadynamicalmodel(e.g.Widrow et al. based surveys are Nano-JASMINE and ESA’s cornerstone 2008). There are arguably five different types of dynamical galaxymodel,althoughsometimeswherethelineofdistinc- tion is drawn can be ambiguous. ⋆ E-mail:[email protected] † E-mail:[email protected] Moment based methods find solutions of the Jeans 2 J. A. S. Hunt and D. Kawata equation that best fit the observed moments and ased towards any inaccuracies from this previous model. minimise χ2 (e.g. Young 1980; Binney et al. 1990; The next generation of Milky Way model should be built Magorrian & Binney 1994; Magorrian 1995; Cappellari directly from observational data of the Milky Way, and 2008; Cappellari et al. 2009). The main drawback of this flexibleenough for fitting heterogeneous data. method is that there is no guarantee that there will NMAGIC,developedbydeLorenzi et al.(2007),isthe be a positive distribution function with the required ve- first algorithm to improve upon the initial M2M algorithm locity moments. It is also usually restricted to spheri- by adding the ability to include observational errors in the cally symmetric models as the symmetry allows simpli- constraints. This is an important step forward as it allows fied assumptions to be made. Distribution function based realobservational datatobeusedas constraints.NMAGIC methods fit the distribution function f(r,v) to the data wasalsothefirstM2Malgorithmtousevelocityconstraints, directly. The methods have been applied to spherical intheform of lineofsight spectra. NMAGIChasnowbeen or integrable systems (e.g. Dejonghe 1984; Bishop 1987; applied to several observed galaxies (e.g. deLorenzi et al. Gerhard 1991; Hunter& Qian 1993; Merritt & Tremblay 2008, 2009; Das et al. 2011). Morganti & Gerhard (2012) 1994; Kuijken 1995; Merritt 1996; DeBruyneet al. 2000). also made a recent improvement to the field of M2M mod- Perturbation theory can be used to extend the method to elling, by developing Moving Prior Regularization (MPR) near integrable potentials (e.g. Matthias & Gerhard 1999; which can replace the Global Weight entropy Regulariza- Binney 2010). Schwarzschild’s method works by comput- tion (GWR). Morganti & Gerhard (2012) show that MPR ing a large number of orbits evolved over many orbital is beneficial to accuracy and smoothness in phase space periods in a fixed potential. Information is collected in distributions, and in some circumstances can converge to an orbit library, and they are weighted to produce the a unique solution, independent of the choice of the initial best fit to the target model (e.g. Schwarzschild 1979, 1993; model. To examine M2M’s performance against previous Cretton et al. 1999; Gebhardt et al. 2003; Krajnović et al. methods,Long & Mao(2012)haveperformedadirectcom- 2005; Cappellari et al. 2006; Thomas et al. 2009). This parison between M2M and the better known Schwarzchild method has the advantage of not requiring the distribution method with regard to calculating the mass to light ratios function or the other integrals of motion, and rarely, the of several elliptical and lenticular galaxies. distribution function may even be recovered (Häfner et al. ThesepreviousM2Malgorithmsuseadistributionfunc- 2000). This method is not restricted by symmetry, but due tion or binned density distribution. However, the data that tocomplexityitisusuallyonlyusedforaxisymmetricmod- Gaiaandtherelatedsurveysreturnwillbeintheformofin- els. Recently van den Bosch et al. (2008) have developed a dividualstellardata.Thereforewehavedesignedaparticle- triaxial Schwarzschild method and applied it toNGC 4365. by-particle M2M algorithm that compares the observables Torus methods are very similar to orbit based methods, at the location of each star (or the target particle) with and are often labelled within the same category. The key the model observables at the same locations, and adjusts difference between torus modelling and orbit based mod- the weights in the same fashion as the original algorithm elling is that while in orbit based modelling, the orbits are from Syer& Tremaine (1996). In this paper, we present time series of phase-space points, in torus modelling, these proof of concept of the particle-by-particle M2M by recre- are replaced by orbital tori (e.g. McMillan & Binney 2012; ating disc galaxies, generated with a Tree N-Body code, Binney 2012a,b). For a more detailed explanation, and a GCD+(Kawata & Gibson2003).Ouralgorithmusesaself- list of advantages of torus methods over orbit methods, see consistentpotential,whichevolvesovertimealongwiththe Binney & McMillan (2011). Finally N-body models, which particle weights. We also show a model constructed from a arethesimplest toconstruct,arebasedongravitational at- partial target data set, demonstrating that the observables traction between ‘N’ bodies and can be collisional, or col- of the target galaxy do not have to cover the whole galaxy lisionless. The key assumption of stellar dynamics in the for M2M to work. This is the first step towards the real Galaxy is that these stellar systems are collisionless, hence observational data from Galactic surveys, where the infor- collisionless N-body models are good approximations for mationwillbeprovidedforalimitedregion ofthesky,with galactic dynamics(e.g. Dehnen& Read 2011). a more complicated selection function due to the dust ex- WewilldemonstrateanN-bodymethodthatallows us tinction,crowding and stellar populations. Thepaperis or- torecoveramodelofthedesiredgalaxy,withsomeflexibility ganised as follows. Section 2 describes thetraditional M2M ontheinitialconditions.Ourmethodisbasedupontheorig- methodandSection3describesthemethodsbehindourpar- inal made-to-measure (M2M) method by Syer& Tremaine ticle based adaptation. Section 4 shows the performance of (1996)whichiscapableofconstructingN-bodyequilibrium the particle-by-particle M2M for recreating the target disc systemsbymaximising alinearcombination of theentropy, system.InSection5wediscusstheaccomplishmentsofthis and minimising χ2, the the mean-square deviation between paper, and describe thenext stages of ourwork. the observables and the model. The M2M algorithm has been improved upon by de Lorenzi et al. (2007), Dehnen (2009),Long & Mao(2010)andMorganti & Gerhard(2012) 2 THE M2M ALGORITHM and has been used for a variety of tasks, as detailed below. Bissantz et al. (2004) apply the M2M algorithm from In this section, we will give a brief description of the Syer& Tremaine (1996) to the Milky Way for the first M2M algorithm as detailed in Syer& Tremaine (1996), time, and create a stellar dynamical model of the Milky deLorenzi et al. (2007) and Long & Mao (2010), which Way’s barred bulge. The model is constrained however by forms the base for our work. The M2M algorithm works a previously constructed model of the Milky Way from by calculating observable properties (observables hereafter) Bissantz & Gerhard (2002), so this new model will be bi- from the model and the target, and then adapting parti- A Particle-By-Particle M2M algorithm 3 cleweights suchthatthepropertiesof themodelreproduce This is useful if the total mass of the target system is one thoseofthetarget.Thetargetcanbeintheformofadistri- of the constraints. We do not impose this restriction as we butionfunction,anexistingsimulation,orrealobservational wishtobeabletocreateasystemwithadifferenttotalmass data. The model is always an N-body system. from theinitial model. The observables of thetarget system are described by Once the new entropy term is introduced to the force of change, equation (4) is replaced by Y = K (z)f(z)d6z, (1) j j Z dw (t)=−ǫw (t) Kj[zi(t)]∆ (t)−µδS (t) , (9) wherej representseachindividualobservable,z=(r,v)are dt i i " Yj j δwi # j thephasespacecoordinates,f(z)isthedistributionfunction X or of thetarget galaxy and K is a known kernel.Observables j can come in many forms, including surface or volume den- d K [z (t)] w (t)= − ǫw (t) j i ∆ (t) sities, surface brightness and line of sight kinematics. The dt i i " Yj j corresponding observablefor themodel takesthe form Xj w (t) N + µ ln i +1 , (10) yj = wiKj[zi(t)], (2) (cid:18) (cid:18) wî (cid:19) (cid:19)# i=1 X forthemostcompleteform.NotethatZ hasbeenreplaced j wherewi aretheparticleweightsandzi arethephasespace byYj dueto themaximisation of equation (5). coordinates of the model’s i-th particle. We then calculate It is shown in Syer& Tremaine (1996) and thedifferenceintheobservablesofthetargetandthemodel, deLorenzi et al. (2007) that fluctuations in equation y (t)−Y (3) may be reduced by employing temporal smoothing, ∆ = j j. (3) j Y effectively boosting N without drastically increasing j computation time. This is achieved by replacing ∆ (t) in j We then use this ∆j to determine the so called force of equation (4) with ∆˜ (t), where j change with theequation ∞ dw (t)=−ǫw (t) Kj[zi(t)]∆ (t), (4) ∆˜j(t)=α ∆j(t−τ)e−ατdτ, (11) dt i i Z j Z0 j Xj with α being small and positive. This ∆˜ (t) can be calcu- j where Zj so far is an arbitrary constant, and the factor lated from the differential equation K /Z can be thought of as the degree to which the i-th i j d∆˜(t) particle contributestothej-th observable.ǫ isa parameter =α(∆−∆˜). (12) enablingustocontroltherateofchange.Thelineardepen- dt dence of equation (4) upon w , coupled with the provision Thistemporalsmoothingeffectivelyincreasesthenumberof i that a small enough ǫ is used, ensures that the weights do particles from N to notbecome negative.Syer& Tremaine(1996)showaproof t1 of convergence for equation (4) providing that the system N =N 2 , (13) eff ∆t starts close to thetarget. If N>J, i.e. the number of the model particles, N, where∆tisthelengthofthetimestepandt1 =(ln 2)/αis 2 greatly exceeds thequantityof available constraints, J, the the half life of the ghost particles. Syer& Tremaine (1996) differential equation (4) is ill-conditioned. Syer& Tremaine showthatexcessivetemporalsmoothingisundesirable,and (1996)suggestremovingthisillconditioningbyintroducing should be limited to α>2ǫ. entropy,bymaximising thefunction Theparametersǫ,µandαmustbedeterminedviapa- rametersearch. Wewill discuss ourchoice of theseparame- 1 F =µS− χ2, (5) ters in Section 3.4. 2 where χ2 = ∆2, (6) 3 PARTICLE-BY-PARTICLE M2M j j X This section describes our original adaptation to the M2M and µ is a parameter to control the regularization. The en- algorithm. The majority of the methodology remains the tropy is given by same as described in Section 2, with the most substantial w difference involving the Smoothed Particle Hydrodynamics S =− wiln wî , (7) (SPH) kernel (e.g. Lucy 1977; Gingold & Monaghan 1977), Xi (cid:18) i(cid:19) which will be described in Section 3.1. Syer& Tremaine where wˆ are the priors, a predetermined set of weights, (1996)usedakernelwheretheydividethecoordinatespace i normallyidenticaltoeachothersuchthatwˆ =M/N,where into bins. For example, for the density at the j-th bin with i M is the total mass of the system and N is the number of volumeV ,thekernel,K (r ),issettobe1/V ifr iswithin j j i j i particles.Thesystemcanbeennormalised(deLorenzi et al. the j-th bin. If r is outside the j-th bin, K (r ) = 0. Be- i j i 2007) such that cause K (r ) and K (r ) are the same if r and r are in j 1 j 2 1 2 thesamebin,thislimitstheresolutiontothebinsize.How- N ever as mentioned in Section 1, our ultimate target is the w =1. (8) i Milky Way, and the observables are not binned data, but i=1 X 4 J. A. S. Hunt and D. Kawata the position and velocity of the individual stars which are of the observables, using the same kernel. For example for distributedratherrandomly.Tomaximisetheavailablecon- radial velocity; straints,we evaluatetheobservablesat theposition ofeach N star and compare them with the N-body model, i.e. in a δv = (v −v )m W(r ,h ), (19) t,r,j t,r,k t,r,j t,k kj j particle-by-particlefashion. Tothisendweintroduceaker- nel often used in SPH, Xk=1 where v is the radial velocity of the k-th target parti- t,r,k Kj(ri,vi)=W(|ri−rj |,hj), (14) cle and vt,r,j = (vt,x,jxt,j +vt,y,jyt,j)/(x2t,j +yt2,j)21 is the radial velocity of thetarget system. Equation 19 represents whereourW(r,h)isasphericallysymmetricsplinefunction the weighted mean of the relative velocities of the target given by particles within h of the target particle j. j 1−6(r/h)2+6(r/h)3 if 06r/h61/2, N W(r,h)= πh83 × 2[1−(r/h)]3 if 1/26r/h61, δvr,j = (vr,i−vt,r,j)miW(rij,hj) (20) 0 otherwise.  i=1 X (15) is similarly calculated from the model particles. The same asshowninMonaghan & Lattanzio(1985),wherer=|r − i format is applied for thevertical and rotational velocities. r |. j We then describe the difference in the observables i.e. Below, we describe our particle-by-particle M2M, con- equation (3).For density; sidering that thetarget system is an N-body system whose particle position and velocity are known without any error. ∆ = ρj(t)−ρt,j. (21) Of course in the real data of the Galaxy, there are com- ρj ρ t,j plicated observational errors and selection functions, which For velocity, we normalised them by the target density be- often depend on stellar population and dust extinction. In causeofthedensitydependenceintroducedinequations(19) thispaper,weignoretheseandconsideranidealisedsystem and (20),and therefore for theradial case; foratarget.AsdescribedinSection1,theaimofthispaper istodemonstratehowournewM2Mworksandthepotential ∆ = δvr,j(t)−δvt,r,j. (22) offutureapplication totheGalacticdisc.Webelowassume v σvr,jρt,j thatthetargetsystemconsistsofasinglepopulation,which Notethatσisnotanobservationalerror,butjustanormal- weshallrefertoasparticles,andwhosepositionandvelocity isation constant which wehavearbitrarily set to 10 km s−1 are known without errors. in our demonstration in Section 4. Because ∆ and ∆ are normalised differently, we ρj vj modified their contribution to the force of change by in- 3.1 Method troducinga new parameter ζ such that for oursimulations, We use the kernel of equation (15) to calculate the density equation (10) becomes; at the target particle locations, r , of both the target and j d W(r ,h ) the M2M model. Hereafter we replace the particle weights, m (t)= − ǫm (t) ij j ∆˜ (t) wi,withtheirmassesmi duetoouradoptionofself-gravity dt i i " j ρt,j j,ρ X intheparticle-by-particleM2M.Forexample,thedensityof W(r ,h ) thetarget at rj is evaluated by, + ζ ξr σvirj,j j ∆˜vr,j(t) N ρt,j = mt,kW(rkj,hj), (16) + ξ W(rij,hj)∆˜ (t) z σ vz,j Xk=1 vz,j wheremt,k isthemassofthetargetparticle,rkj =|rk−rj |, + ξ W(rij,hj)∆˜ (t) and hj is thesmoothing length determined by rot σvrot,j vrot,j ! hj =η(cid:18)mρtt,,jj(cid:19)1/3, (17) + µ(cid:18)ln(cid:18)mmî(it)(cid:19)+1(cid:19)#. (23) whereηisaparameterandwehavesetη=3.InSPHsimu- Weuse the additional individual parameters ξ , ξ , ξ for r z rot lations,avalueofη between2and3areoftenused,andwe the different velocity observables, to allow us to fine tune employtherelatively highervaluetomaximisethesmooth- theircontributionstotheforceofchangeevenfurther.Sim- ness. The solution of equation (17) is calculated iteratively ilartode Lorenzi et al.(2007),wewriteǫasǫ=ǫ′ǫ′′ where until the relative change between two iterations is smaller ǫ′′ is given by than 10−3 (Price & Monaghan 2007). Similarly, ′′ 10 ǫ = , (24) N max W(rij,hj)∆˜ (t) ρj = miW(rij,hj), (18) i j ρj,t ρj Xi=1 for thedensity observab(cid:16)lePonly. (cid:17) from the model particles. The target density ρ is calcu- In the previous works (e.g. Syer& Tremaine 1996; t,j lated only once at the beginning of the M2M simulation, Dehnen 2009; Long & Mao 2010; Morganti & Gerhard and the model density ρ is recalculated at every timestep. 2012), the M2M method is applied to a system in a known j For velocity constraints, we define the following form fixed potential, i.e. using test particles. deLorenzi et al. A Particle-By-Particle M2M algorithm 5 (2007) demonstrate that M2M works with a partially self- consistentpotential,inthatthepotentialiscalculatedevery 25 time steps, setting the particle mass m = w M. How- i i everthisrepeatedsuddenchangeofthepotentialcouldcome with some problems that will bediscussed later. We intend to apply our algorithm to the Milky Way, whose mass distribution is poorly known (e.g. McMillan 2011),andoneoftheaimsofapplyingthedynamicalmodel is to reconstruct the mass distribution. Therefore we use a self-consistent potential, setting the particle weight, w , to i the mass, m , allowing the potential to change along with i themodelobservablesandallowingustorecoversimultane- Figure 1.Theendresult(t=2Gyr)ofanN-bodydiscgalaxy ouslythepotentialalongwiththemassandvelocityprofiles. simulation. It had a scale length of 3 kpc initially. This will be In this paper however, we focus on the disc. We ignore the usedasthetargetsystemasshownintheSection4.Theleftand bulge or halo stars, and assume that the dark matter po- rightpanelsshowtheface-onandend-onviewsrespectively. tential is known for this initial demonstration. Note that the previous studies are mainly focused on elliptical galax- halo density profile is taken from Navarroet al. (1997); ies, i.e. systems dominated by velocity dispersion, but not stronglyrotationsupported.Recreatingadiscgalaxywitha ρ = 3H02(1+z )3 Ω0 δc , (26) self-consistent potential has been attempted once before by dm 8πG 0 Ω(z0)cx(1+cx)2 Deg (2010), who highlights some difficulties with an M2M where δ is the characteristic density described by c method that employs self-gravity. He uses a grid to calcu- Navarroet al. (1997). The concentration parameter c = latetheobservables,whichmakeshismethoddifferentfrom r /r and x = r/r , where r is the radius inside 200 s 200 200 ours. which the mean density of the dark matter sphere is equal Oneoftheproblemsarisingfromusingaself-consistent to 200ρ and given by; crit potential as mentioned by Deg (2010) is that the tempo- 1 −1 ral smoothing, which worked well in fixed potential M2M r =1.63×10−2 M200 3 Ω0 3 (1+z )−1h−1kpc. methods, is problematic when used with self-gravity. The 200 (cid:18)h−1M⊙(cid:19) (cid:20)Ω(z0)(cid:21) 0 temporalsmoothingreducesshotnoisebyaveragingthe∆j (27) back along their orbits, which is fine with test particles in WeuseM200 =1.75×1012 M⊙,c=20, Ω0 =0.266, z0 =0 a fixed potential because the orbits are fixed. However in a and H =71 km s−1Mpc−1. 0 self-consistentpotential,thepotentialandthereforeparticle Thestellardiscisassumedtofollowanexponentialsur- orbits change with time, and thus the temporal smoothing face density profile: breaks the self-consistency. Therefore we should be aware that self-gravity M2M models areverysensitive toinstabil- ρ = Md sech2 z e−R/Rd, (28) d 4πz R z ities, and we see substantial disruption when the smooth- d d (cid:18) d(cid:19) ing is first turned on. A way to mitigate this damage due where z is the scale height of the disc and R is the scale d d to the temporal smoothing is described in Section 3.2. In length.Ourtargetdischasz =0.35 kpcandR =3.0 kpc. d d light of thiswe investigated thepossibility of runningmod- ThedischasamassofMd =3.0×1010 M⊙ andconsistsof els without temporal smoothing. However all models had 105 particles, with each particle having a mass of 3.0×105 tobesubstantiallyunder-regularizedtorecoverthevelocity M⊙.Weuseafixedsoftening lengthofε=1.05 kpc,which profiles shown in Section 4, which leads to the continuous we also use for the M2M modelling runs. The velocity dis- fluctuation of the weights, similar to the problems of the persion for each three dimensional position of the disc is under-regularization discussed in Section 4.2. computed following Springelet al. (2005) to construct an We use a standard Euler method for the integration of almost equilibrium condition. We use a high value of the the weight change equation and a leapfrog time integrator free parameter f = σ /σ = 3, which controls the ratio R R z foradvancingtheparticles.Wealsouseindividualtimesteps between the radial and vertical velocity dispersions, to de- fortheparticles,andonlyapplytheweightchangealgorithm liberatelysuppressstructureformationandcreateasmooth, to the particles whose position and velocity are updated at almostaxisymmetricdiscforthisinitialtest.Ourtargetsys- eachtimestep.Thetimestepforeachparticleisdetermined tem is a relatively smooth disc galaxy evolved over 2 Gyr, by as shown in Fig. 1, and it is used for all models in Section 4. 1 dt =C 0.5hi 2 , (25) We set up the initial conditions of the model disc with i DYN(cid:18)|dvi/dt|(cid:19) the same parameters and method, but use a different scale length from that of thetarget galaxy. with C =0.2. DYN 3.3 Procedure 3.2 Target System Setup Thesuddenchangeinpotentialcausedbythechangingpar- Our simulated target galaxy consists of a pure stellar disc ticle weights induces instabilities and potentially unwanted withnobulgeandastaticdarkmatterhalo,setupusingthe structureformation. This effect can be reduced by dividing method described in Grand et al. (2012). The dark matter the modelling process into a series of stages, each with a 6 J. A. S. Hunt and D. Kawata Table 1.M2Mmodelparameters Model Rd,ini (kpc) ǫ′ µ α ζ χ2ρ χ2vr χ2vz χ2vrot Notes A 2.0 0.1 5×105 0.2 0.05 0.0846 7.291 0.918 6.502 Fiducial B 2.0 0.1 104 0.2 0 0.0831 9.599 1.074 10.873 NoVelocity C 2.0 0.1 104 0.2 0.05 0.0912 8.275 1.069 7.464 D 2.0 0.1 105 0.2 0.05 0.0875 7.914 1.005 7.087 E 2.0 0.1 106 0.2 0.05 0.0894 7.099 0.893 6.440 F 2.0 0.1 107 0.2 0.05 0.223 9.395 1.130 9.960 G 2.0 0.1 108 0.2 0.05 0.407 17.291 2.107 17.701 H 5.0 0.1 5×105 0.2 0.05 0.100 10.839 1.414 10.394 I 6.0 0.1 5×105 0.2 0.05 0.111 12.381 1.604 12.094 J 1.5 0.1 5×105 0.2 0.05 0.101 7.849 0.972 6.896 K 2.0 0.1 5×105 0.2 0.05 0.0924 7.309 0.911 6.509 PartialData slightly different level of M2M algorithm. This reduces the values. Note that these parameters are calibrated for this magnitude of the change in potential at any one time. We specifictargetsystem.Itislikelythatweneeddifferentcali- also set a limit on the maximum change in mass any parti- brationfordifferenttargets.Howeverwhatwelearnedfrom clecan experienceinonetimestep.Wesetthislimit toten theparametersearchshouldbeusefulforfutureapplications percent of that particles mass. and developmentsof theimproved version. Initially the model is allowed to relax in a pure self- ǫ provides the balance between the speed of conver- gravityenvironmentwithnoM2Mconstraintsfor0.471Gyr gence, and the smoothness of the process. In this case, we (our N-body code time unit). This relaxation period is im- find that when ǫ′ > 0.1, the weights change too rapidly, portant, as applying the M2M algorithm before the model which induces the sudden potential changes and therefore has settled generates the aforementioned instabilities. Al- more instabilities. This leads to a general decrease in the though our M2M algorithm was still capable of recovering final level of accuracy of both density and velocity profiles. thedesiredprofile,thetimescaleneededwasdrastically in- If ǫ′ 6 0.1 convergence can be achieved and the particle creased because the model had to smooth out again before weights experience a much smoother evolution. However if convergence took place if we turned on the M2M without ǫ′ is too small, the oscillations generated by the temporal therelaxation period. smoothing take too long to damp down, which drastically After this period of relaxation, the M2M algorithm is increases the length of the simulation. In the end, we have activatedandrunswithouttemporalsmoothingforafurther chosen ǫ′ = 0.1 as a balance between accuracy and sim- 1.413 Gyr, which allows thedensity and velocity profiles to ulation time. With more computing power available to us converge quickly. During this time, the contribution of the we would consider running a lower value of ǫ. However if velocity constraint is increased linearly from 0 up until our ǫ′ ≪0.1, it is possible the model will not show any signs of desired ζ. This allows the density profile to converge first. convergenceas theweight change is too slow. We found this slow increase in the velocity constraints to The choice of α, which controls the strength of the be important, because if the velocity constraints were in- temporal smoothing, should depend upon the choice of ǫ troduced simultaneously at full strength, we find the large (α > 2ǫ). Note that ǫ = ǫ′ǫ′′ and ǫ′′ is defined by equation weight changes induce the sudden potential change men- (24). We find that our modelling is not sensitive to α and tioned earlier, which is strong enough to disrupt thedisc. we set α=0.2 in thispaper. Then, after 1.884 Gyr, the temporal smoothing is ζ (and individual ξ) controls the level of the velocity turned on. When the M2M modelling was run with tem- constraints. It is important to strike a balance between the poral smoothing from thebeginning, themass profileexpe- densityandvelocity constraints, becauseif thelevelofcon- rienced large oscillations. The modelling then continues in straints are unbalanced one will dominate in the change of thisstateforaslongasisdesired.OurM2Mmodelsarerun weight and theother observables will not converge. Wecan for a period of 10 Gyr. choose a suitable ζ (and/or ξ ,ξ and ξ ) by comparing r z rot the magnitudes of the individual terms of the right-hand- side of equation (23). We set ζ such that the contribution 3.4 Parameter Calibration of the velocity constraint to the force of change equation As discussed in Syer& Tremaine (1996), deLorenzi et al. is the same magnitude as, or slightly less than, the density (2007), Long & Mao (2010) and Morganti & Gerhard constraints. The individual velocity components may then (2012), the choice of parameters are crucial for the success befinetunedwith ξ . Forour simulations, we findthat the j ofM2Mmodelling.Inthissection,wewilldiscussourchoice followingparametersetworkswell:ζ =0.05,ξ =1,ξ =10 r z oftheparameters;ǫ,α,ζ andµ,andhowwecalibratethese and ξrot=1. A Particle-By-Particle M2M algorithm 7 µcontrolsthestrengthoftheregularization.Wediscuss the importance of µ in greater detail in Section 4.2. In our fiducialmodel shown in Section 4 we adopt µ=5×105. 4 PARTICLE-BY-PARTICLE M2M RESULTS In this section we present the results from our modelling of our target disc galaxy. We will first show the results for ourfiducialmodel, andthencompareit withamodelusing onlydensityconstraints.WeranmultipleM2Mmodelswith different parameters, which can be seen in Table 1, where R is the initial scale length of the model disc. We only d,ini usethe observables within theradius of 10 kpc. 4.1 Fiducial Model Inthis section wepresent Model A,ourfiducial modelcon- structed with the parameters described in Section 3.4, and showninTable1.WestartfromanN-bodydiscwithascale lengthof2kpc,recreatingthetargetdisc(R =3 kpc)with d ourparticle-by-particleM2M,evolvingthemodelfor10Gyr. Fig.2showstheradialprofilesofthesurface density,radial and tangential velocity dispersion, and the mean rotational velocity. The final profiles of Model A reproduce the profiles of the target system remarkably well. Note that these radial profiles are not direct constraints of the particle-by- particleM2M. Especially it israthersurprisingthattheve- locitydispersionprofilesarerecovered.Wethinkthatthisis becausetheparticle-by-particleM2Mforcesthemodelparti- clestofollowthevelocitydistributionofthetargetparticles. We also have no constraints on the total mass of the disc. Note also that the assumed velocity constant is 10 km s−1 inequation(22)yetthevelocityprofilesarereproducedata levelmuchless than10 km s−1.This isnot surprising how- Figure 2. Initial (red dotted), final (green dashed) and target ever, because we have different normalisations for density (black solid), density profile (upper), radial velocity dispersion and velocity, and adjust ζ and ξ to balance their contri- (upper middle), vertical velocity dispersion (lower middle) and butions in equation (23) making the choice of σ arbitrary. rotational velocity (lower) for Model A. The initial model has a Therefore, σ is not indicating an error, but is merely a vr,t scalelengthof2kpc,thetargetmodelhasascalelengthof3kpc. constant value for normalisation. In this paper, we do not includeanyerror.Weplantoaddmorerealisticerrorsinfu- tureworks.Fig.3showstheweightevolutionforaselection ofparticlesfrom ModelA.Weight convergenceisadequate, however it is not as smooth as desired. We find that the particle weight evolution is less smooth for the case where velocity observables have been added. Fig. 4 shows the χ2 evolutionforeach oftheobservables.Forallobservableswe use ∆2 χ2 = X, (29) X N P r where ∆ is equivalent to equations (21), i.e. X = ρ, and X (22),i.e.X =v.Thisisaslightlyunusualdefinitionofχ2for thevelocityobservables.Notethatweincludeonlyparticles within 10 kpcand N is the numberof target particles sat- r isfying this criteria. In Model A, χ2 values rapidly decrease until 2 Gyr, from which point there is almost no improvement.The final valuesof χ2 are also shown in Table 1. Figure 3.The weight evolution foraselection of particles from In comparison we show Model B, with the same initial ModelA. conditions and target with the velocity constraints turned off. We find that µ = 5×105 cause over-regularization for thiscase,andhastobereducedincompensationtoµ=104. 8 J. A. S. Hunt and D. Kawata Figure4.Timeevolutionofχ2fordensity(upper),radialveloc- Figure 5.Same as Fig. 2, but forModel B which uses onlythe ity(uppermiddle),verticalvelocity(lowermiddle)androtational densityobservableasaconstraint. velocity(lower)forModelA. The density observableappreciates a slightly lower valueof Fig. 5 shows the density and velocity profiles for Model B. µ. Thefinalmodel-densityprofileresemblesthetarget.Dueto Althoughthereisnotavastdifferencebetweenthefinal the lack of velocity constraints, while the velocity profiles values of χ2 for µ=104,105,106, Model C with µ=104 is do improve, they do not resemble the target. A compari- foundtobeaninappropriatemodelbecauseofitspoorcon- son between Fig. 2 and Fig. 5 demonstrates how the veloc- vergence. Fig 7 shows thetime evolution of χ2 in Model C. ity constraints improve our reproduction of the dynamical Fig.7showsoscillatorybehaviour.Fig.8showsthetimeevo- properties of thetarget. lutionoftheweightfortheparticlesselectedinFig.3.Com- parison betweenFigs. 3and8demonstratesthatµ=104 is toolowtosuppressthelargeamplitudeofthefluctuationsin 4.2 Effect of Regularization theparticleweights.Theweightsoftheparticleskeepchang- Similar to thepreviousstudies (e.g. Syer& Tremaine 1996; ing and do not converge. Therefore we judge that µ = 104 deLorenzi et al.2007;Long & Mao2010),wealsofindthat is unacceptable for recreating thetarget system. careful choice of the value of µ is key to obtain conver- Fig. 9 displays the distribution of particle weights at gencetoagoodmodel,andreproducethegivenobservables. the final time for Models A and C. The histogram shows a Therefore we discuss in this section how µ affects the mod- wider tail, and lower peak for the under-regularized Model elling. We performed multiple models with the same initial C compared to our fiducial Model A. This is expected be- conditions and parameters as Model A, except the value of causeahigherµrestrictsparticlesfrommovingfarfromthe µ (see Models C-G in Table 1). Fig. 6 shows the χ2 for the initial mass used as a prior. As a result, Model A shows a densityandvelocityat thefinaltime(t=10 Gyr).Thefig- narrower distribution and thus a higher peak close to the uredemonstrates aslowimprovementfor thethreevelocity initial value of w . Fig. 9 also demonstrates that µ=104 is i observableswithanincreasingµupuntilavalueofapproxi- less favourable. matelyµ=106,abovewhichgoodnessoffitdropsoffagain. If we examine substantial over-regularization, i.e. a A Particle-By-Particle M2M algorithm 9 Figure 6. Accuracy of our final M2M model dependent on µ Figure 7.SameasFig.4,butforModelC,withµ=104. as determined by χ2 for density (upper), radial velocity (upper middle), vertical velocity (lower middle) and rotational velocity (lower). 4.3 Different Initial Conditions We also tested thealgorithm on the same target, using ini- highervalueofµ,itiseasytoseethedamagingeffectonthe tial discs with a different scale length, but with the same densityandvelocityprofiles.Fig.10showstheprofilesfrom parametersasModelA.Wehavealreadydiscussedtheben- ModelG,withµ=108,whichshowsthesignificantdiscrep- efits of tailoring µ to the model, so we were not expecting ancy in thedensity and rotational velocity profiles between that these models (Models H-I in Table 1) would recreate the final profiles and the target profiles. The discrepancy theirtarget systemstothesamelevelasModelA.However in the other two velocity observables is not as substantial. for demonstration purposes, we show how the parameter Howeverit is clearly worse when compared with Fig. 2. set in Model A works if the initial conditions are different. Insummary,wefoundthatwerequiredregularizationof When we started from a higher initial scale length (Model aroundµ=105−106asacompromisebetweenthegoodness HwithR =5 kpcand ModelIwith R =6 kpc)we d,ini d,ini offit,andthesmoothnessoftheχ2 andparticleweightevo- attained a reasonable reproduction of the target, However, lution. Both µ=107 and µ=108 show over-regularization thefinalχ2issystematicallyhigherthanModelA(seeTable and the density profiles associated with those values con- 1). Fig. 11 shows the profiles from Model H, which slightly verge to an incorrect profile. µ = 104 shows large oscil- disagree with the targets. This seems to be due to over- lations in both χ2 and particle weights, and convergence regularization, and we would need to adjust µ in order to is not reached. Anything in the range of µ = 105 − 106 obtainabettermodel.Ontheotherhand,whenwestarted appears appropriate and hence our fiducial model adopts from a lower initial scale length (Model J with R =1.5 d,ini µ = 5×105. As can be seen from Table 1, we find under- kpc,theprofilesareshowninFig.12),χ2 wasonlyfraction- regularization is preferable to over-regularization. This is allyworsethanthefiducialcaseModelA(seeTable1).This also the case in previous literature (e.g. deLorenzi et al. demonstrates that it is better to set the initial disc with a 2008; Morganti & Gerhard 2012) implying this is a generic smaller scale length. In the application to the real obser- feature of M2M and not intrinsic to any specific algorithm. vational data of the Milky Way, we do not know the right 10 J. A. S. Hunt and D. Kawata Figure 8.SameasFig,3,butforModelC,withµ=104. Figure 10.SameasFig.2butforModelG,withµ=108. within a 10 kpc sphere around a point in the plane, 8 kpc away from theGalactic centre. Fig. 13 shows the final profiles for Model K, which re- Figure 9. Distribution of particle weights for Model A (solid) producesthetargetprofilesreasonablywell.Comparedwith andModelC(dotted)atthefinaltime,t=10Gyr.w0 indicates Fig.2,Fig.13showsonlyaminordiscrepancytothetarget theinitialparticleweights. profiles, mainly in the outer region. Worse performance in the outer region is unsurprising, as the larger the radii, the smaller percentage of the particles orbits are spent within shape of “the target model”. However, we hope that the thesampledarea.Table1showsthefinalvaluesofχ2,which further studies with these target galaxies would help us to displays a better value of χ2 than over-regularized models, understandmoreabouthowtheM2Mmodellingbehavesin andsimilarlevelsofthegoodnessoffittounder-regularized differentcases,andhowweshouldcalibratetheparameters. ones, without the excessive weight oscillations. Model K demonstrates that it is possible to apply our particle-by- particle M2M to a disc galaxy with only a limited selection 4.4 The Partial Data Case of data. BecauseourgoalistoeventuallyuseourmethodwithGaia data, and Gaia will only survey a section of the Galactic 5 SUMMARY AND FURTHER WORK disc, it isimportant totest ouralgorithm on an incomplete data set (Model K in Table 1). In this paper, a simple se- We have developed the initial version of our new particle- lectionfunctionisappliedforthepurposeofdemonstration. by-particle M2M, where the observables are compared at Remember that our models in the previous sections have the position of the target particles, and the gravitational usedonlythedatawithintheradiusof10kpcfromthecen- potential is automatically adjusted bythe weight changeof tre.Inthissectionweadditionallyrestrictedtheobservables theparticles.Thispaperdemonstratesthattheparticle-by-