Perturbed distribution functions with accurate action estimates for the Galactic disc
H. Al Kazwini, Q. Agobert, A. Siebert, B. Famaey, G. Monari, S. Rozier, P. Ramos, R. Ibata, S. Gausland, C. Riviere, D. Spolyar
AAstronomy & Astrophysics manuscript no. f1timedep ©ESO 2020December 15, 2020
Perturbed distribution functions with accurate actionestimates for the Galactic disc
H. Al Kazwini , Q. Agobert , A. Siebert , B. Famaey , G. Monari , S. Rozier , P. Ramos , R. Ibata ,S. Gausland , C. Rivière , and D. Spolyar Université de Strasbourg, CNRS UMR 7550, Observatoire astronomique de Strasbourg, 11 rue de l’Université, 67000Strasbourg, France Université de Strasbourg, UFR de Mathématique et d’Informatique, 7 rue René Descartes, 67084 Strasbourg, France The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, 10691Stockholm, SwedenReceived; accepted
ABSTRACT
In the Gaia era, understanding the effects of perturbations of the Galactic disc is of major importance in the contextof dynamical modelling. In this theoretical paper, we extend previous work in which, making use of the epicyclicapproximation, the linearized Boltzmann equation had been used to explicitly compute, away from resonances, theperturbed distribution function of a Galactic thin disc population in the presence of a non-axisymmetric perturbationof constant amplitude. Here we improve this theoretical framework in two distinct ways in the new code that we present.First, we use better estimates for the action-angle variables away from quasi-circular orbits, computed from the
AGAMA software, and we present an efficient routine to numerically re-express any perturbing potential in these coordinateswith an accuracy well below the percent level. The use of more accurate action estimates allows us to identify resonancessuch as the outer 1:1 bar resonance at larger azimuthal velocities than the outer Lindblad resonance, and to extend ourprevious theoretical results well above the Galactic plane, where we explicitly show how they differ from the epicyclicapproximation. In particular, the displacement of resonances in velocity space as a function of height can in principleconstrain the 3D structure of the Galactic potential. Second, we allow the perturbation to be time-dependent, therebyallowing us to model the effect of transient spiral arms or of a growing bar. The theoretical framework and toolspresented here will be useful for a thorough analytical dynamical modelling of the complex velocity distribution of discstars as measured by past and upcoming Gaia data releases.
Key words.
Galaxy: kinematics and dynamics – Galaxy: disc – Galaxy: solar neighborhood – Galaxy: structure –Galaxy: evolution – galaxies: spiral
1. Introduction
The natural canonical coordinate system of phase-space forGalactic dynamics and perturbation theory is the system ofaction-angle variables (Binney & Tremaine 2008). In an ax-isymmetric system in equilibrium, Jeans’ theorem impliesthat the phase-space distribution function (DF) of any stel-lar (or dark matter) component can be expressed solely as afunction of the actions, that are labelling the actual orbits(e.g., Binney & Piffl 2015; Cole & Binney 2017). However,the effect of various perturbers of the potential (e.g., theGalactic bar and spiral arms) must be included in this pro-cess, together with the response of the distribution function.Within the resonant regions, to fully capture the behaviourof the DF, one needs to construct for each perturber neworbital tori, complete with a new system of action-anglevariables (e.g., Monari et al. 2017a; Binney 2020a,b). Awayfrom resonances, however, one can simply use the linearizedBoltzmann equation. This also allows us to accurately iden-tify the location of resonances.This is particularly important in the context of the inter-pretation of recent data from the Gaia mission (Gaia Col-laboration et al. 2018a, 2020), which revealed in exquisite details the fine structure of stellar action space (e.g., Tricket al. 2019; Monari et al. 2019b,a). While the existence ofmoving groups of dynamical origin had been known fora long time in local velocity space around the Sun (e.g.,Dehnen 1998; Famaey et al. 2005), Gaia revealed theirstructure in exquisite details (Ramos et al. 2018), as well asan estimate of their age distribution (Laporte et al. 2020),together with the shape of the global velocity field awayfrom the Sun within the Galactic disc (Gaia Collaborationet al. 2018b). One additional major finding of Gaia has beenthe existence of a local phase-spiral in vertical height vs.vertical velocity in the Solar neighbourhood (Antoja et al.2018) which might be related to a vertical perturbation ofthe disc by, e.g., the Sagittarius dwarf galaxy (e.g., Laporteet al. 2019; Binney & Schönrich 2018; Bland-Hawthorn &Tepper-Garcia 2020).In previous theoretical work, Monari et al. (2016) –hereafter M16 – explicitly computed the response of an ax-isymmetric DF in action space, representing a Milky Waythin disc stellar population, to a quasi-stationary three-dimensional spiral potential, using the epicyclic approxi-mation to model the actions, which is only a valid approxi-mation for quasi-circular orbits in the thin disc. It was no-
Article number, page 1 of 14 a r X i v : . [ a s t r o - ph . GA ] D ec &A proofs: manuscript no. f1timedep tably shown that the first order moments of the perturbedDF then give rise to non-zero radial and vertical bulk flows(breathing modes). However, to treat perturbations awayfrom the disc, which is particularly important in the Gaiacontext, one cannot make use of the epicyclic approxima-tion to compute action-angle variables. Moreover, it is well-known that spiral modes in simulations can be transient,remaining stationary for only a few rotations (Sellwood &Carlberg 2014), and the response of the disc should be dif-ferent during the period of rising or declining amplitude.The same can be true for the bar, whose amplitude andpattern speed can also grow or vary with time (e.g., Chibaet al. 2020; Hilmi et al. 2020).Here, we improve on this previous modelling of M16in two ways. First, we will use a better estimate than theepicyclic approximation for the action-angle variables, re-lying on the “Torus Mapping" method of McGill & Binney(1990) to convert from actions and angles to positions andvelocities, as well as on the “Stäckel fudge" (Binney 2012;Sanders & Binney 2016) for the reverse transformation.This will allow us to extend previous results to eccentricorbits and orbits wandering well above the Galactic plane.The routines developed and presented in this paper will alsobe of fundamental importance to study the vertical pertur-bations of the Galactic disc in further works. Second, wewill let the perturbation evolve with time, thereby allowingus to model the effect of a growing bar.The paper is organized as follows. In Sect. 2, after ashort reminder of the theoretical framework of perturbedDF – already exposed in detail in M16 –, we present themethod used to expand in Fourier series the perturbingpotential within the new action-angle coordinate system.Then, a comparison with the results in the epicyclic case ismade in Sect. 3. In Sect. 4, we explore the temporal treat-ment of the DF for an analytic growth of the amplitude ofthe perturber. Finally, we conclude and discuss the possiblefuture applications of these theoretical tools in Sect. 5. Anappendix is dedicated to the presentation of the code.
2. The perturbing potential and perturbed DF
The canonical phase-space action-angle variables ( J , θ ) ofan integrable system are obtained from a canonical trans-formation implicitly using Hamilton’s characteristic func-tion as a type-2 generating function. The actions J aredefined as new generalized momenta corresponding to aclosed path integral of the velocities along their corre-sponding canonically conjugate position variable, namely J i = (cid:72) v i d x i / (2 π ) . Since this does not depend on time,these actions are integrals of motion, and the Hamiltoniancan be expressed purely as a function of these.It turns out that Galactic potentials are close enough tointegrable systems that actions can be estimated for them.For quasi-circular orbits close to the Galactic plane, withseparable motion in the vertical and horizontal directions,one can locally approximate the radial and vertical mo-tions of an orbit with harmonic motions, which is known asthe epicyclic approximation. The radial and vertical actionsthen simply correspond to the radial and vertical energiesdivided by their respective (radial and vertical) epicyclicfrequency. However, the epicyclic approximation is not validanymore when considering orbits with higher eccentricity, or with a large vertical amplitude. More precise ways ofdetermining the action and angle coordinates have beendevised. They typically differ depending on whether onewishes to (i) transform angles and actions to positions andvelocities, or to (ii) make the reverse transformation frompositions and velocities to actions and angles. In the formercase, a very efficient method is the Torus Mapping methodfirst introduced by McGill & Binney (1990) – see also Bin-ney & McMillan (2016) for a recent overview –, while inthe latter case, a “Stäckel fudge” is generally used (Binney2012; Sanders & Binney 2016).The general idea of the Torus Mapping is to first expressthe Hamiltonian in the action-angle coordinates ( J T , θ T ) ofa toy potential, for which the transformation to positionsand velocities is fully known analytically. The algorithmthen searches for a type-2 generating function G ( θ T , J ) ex-pressed as a Fourier series expansion on the toy angles θ T ,for which the Fourier coefficients are such that the Hamilto-nian remains constant at constant J . This generating func-tion fully defines the canonical transformation from actionsand angles to positions and velocities. For the inverse trans-formation, an estimate based on separable potentials can beused. Such potentials are known as Stäckel potentials (e.g.,Famaey & Dejonghe 2003), for which each generalized mo-mentum depends on its conjugated position through threeisolating integrals of the motion. These potentials are ex-pressed in spheroidal coordinates defined by a focal dis-tance. For a Stäckel potential, this focal distance of thecoordinate system is related to the first and second deriva-tives of the potential. One can thus use the true potential atany configuration space point to compute the equivalent fo-cal distance as if the potential was a Stäckel one, and fromthere compute the corresponding isolating integrals of themotion, and the actions. In the following, we will make useof both types of transformations, namely the Torus Map-ping to express the potential in action-angle coordinates,and the Stäckel fudge to represent distribution functions invelocity space at a given configuration space point.Now, let H ( J ) be the Hamiltonian of the axisymmetricand time-independent zeroth order gravitational potential Φ of the Galaxy. The equations of motion connecting ac-tions J and the canonically conjugate angles θ are thensimply ˙ θ = ∂H ∂ J = ω ( J ) , ˙ J = − ∂H ∂ θ = 0 , (1)with ω being the fundamental orbital frequencies. Thus, fora given orbit, the actions J are constant in time, definingan orbital torus on which the angles θ evolve linearly withtime, according to θ ( t ) = θ + ω t . Jeans’ theorem thentells us that the phase-space distribution function (DF) ofan axisymmetric potential f = f ( J ) is in equilibrium. Inother words, f is a solution of the collisionless Boltzmannequation: d f d t = 0 . (2) Let now Φ be the potential of a small perturbation tothe axisymmetric background potential Φ of the Galaxy,with an amplitude relative to this axisymmetric background (cid:15) (cid:28) . Then the total potential is Φ = Φ + Φ and the Article number, page 2 of 14l Kazwini et al.: Perturbed distribution functions with accurate action estimates J ( kpc kms ) ( k m s ) J R = 50 kpc kms J ( kpc kms ) ( k m s ) J R = 50 kpc kms J ( kpc kms ) ( k m s ) J R = 50 kpc kms J R ( kpc kms ) ( k m s ) J = 1760 kpc kms J R ( kpc kms ) ( k m s ) J = 1760 kpc kms J R ( kpc kms ) ( k m s ) J = 1760 kpc kms Fig. 1.
Variations of a few Fourier coefficients φ jml ( J ) of the bar potential of Sect. 2.4, as J R or J ϕ increase separately at J z = 0 .The actions on the abscissa axis are in kpc km s − . The curves are very smooth, which justifies our use of cubic splines method tointerpolate. DF becomes, to first order in (cid:15) , f = f + f , which is still asolution of the collisionless Boltzmann equation. Inserting f = f + f in Eq. (2), and keeping only the terms of order (cid:15) , leads to the linearized collisionless Boltzmann equation d f d t + [ f , Φ ] = 0 . (3)where [,] is the Poisson bracket. Integrating Eq. (3) overtime, from −∞ to the time t , then leads to f ( J , θ , t ) = (cid:90) t −∞ d t (cid:48) ∂f ∂ J (cid:48) ( J (cid:48) ) · ∂ Φ ∂ θ (cid:48) ( J (cid:48) , θ (cid:48) , t (cid:48) ) , (4)where J (cid:48) and θ (cid:48) correspond to the action and angles in theunperturbed case as a function of time (i.e. constant actions J (cid:48) and angles evolving linearly).Since any function of the angles is π -periodic in theangles, the perturbing potential Φ can be expanded in aFourier series as Φ ( J , θ , t ) = Re (cid:26)(cid:88) n φ n ( J , t ) e i n · θ (cid:27) . (5)Hereafter, we will consider in-plane perturbations such asspiral arms, meaning that we can write the time-varyingFourier coefficients in a non-rotating frame as φ n ( J , t ) = g ( t ) h ( t ) φ n ( J ) , where g ( t ) controls the amplitude of theperturbation and h ( t ) controls its pattern speed, with h ( t ) = e − i m Ω p t , where Ω p is the pattern speed of the per-turbation and m its azimuthal wavenumber, i.e. its mul-tiplicity. We will mostly consider hereafter m = 2 pertur-bations. The vector index n is a triplet of scalar integerindices ( j, k, l ) running in principle from −∞ to ∞ , but inpractice limited to a given range sufficient to approximate the perturbing potential. In the case of an m -fold in-planeperturbation, it is sufficient to take k = m . The main goalof this section will be to express typical non-axisymmetricperturbing potentials originally expressed in configurationspace as such a Fourier series in action-angle space. The al-gorithm that we present in Sect. 2.3 can however be appliedto any perturbing potential, including non-plane symmetricvertical perturbations.Once the potential is expressed as a function of anglesand actions as in Eq. (5), then Eq. (4) becomes f ( J , θ , t ) = Re (cid:26) i ∂f ∂ J ( J ) · (cid:88) n n (cid:90) t −∞ d t (cid:48) φ n ( J (cid:48) , t (cid:48) ) e i n · θ (cid:48) ( t (cid:48) ) (cid:27) . (6)In M16, assuming φ n ( J (cid:48) , t (cid:48) ) = g ( t (cid:48) ) h ( t (cid:48) ) φ n ( J ) , with h ( t (cid:48) ) = e − i m Ω p t (cid:48) , and the amplitude of the perturbing po-tential constant in time at present time ( g ( t ) = 1 ), and zeroand constant in time at −∞ , led to the following explicitsolution for f ( J , θ , t ) , f ( J , θ , t ) = Re (cid:26) ∂f ∂ J ( J ) · (cid:88) n n φ n ( J ) e i θ s , n ω s , n (cid:27) , (7)where we defined θ s , n = n · θ − m Ω p t, (8) ω s , n = n · ω − m Ω p . (9)The subscript "s" stands for "slow", because in proximityof a resonance of the type ω s , n = 0 , the angle θ s , n evolvesvery slowly. One can also immediately see that the above Article number, page 3 of 14 &A proofs: manuscript no. f1timedep
Fig. 2.
Top panel: estimate of the central bar potential ofSect. 2.4 at the Sun’s position in the plane with 41 complexFourier coefficients in Eq. 10 and the cubic splines reconstruc-tion. The vertical line denotes the true value. Bottom panel: rel-ative accuracy. The potential is estimated at the same configura-tion space location (the Sun’s position) for different velocities,with a global accuracy well below the percent level, althoughwith a slight bias towards higher amplitudes in the case of thespiral potential. linearized solution is valid only away from such resonances,since it diverges for these orbits. Orbits near these reso-nances are actually trapped, and for them the determina-tion of the linearly perturbed DF becomes inappropriate.A specific treatment for these resonant regions is required,which was addressed in, e.g., Monari et al. (2017a) withinthe epicyclic approximation, and by Binney (2020a) in amore general context.Using the epicyclic approximation, the Fourier coeffi-cients of a spiral potential have been computed analyticallyin M16 with indices n = ( j, k, l ) running over the values j = {− , , } , k = m = 2 , and l = {− , , } , and theperturbed distribution function away from resonances wasthen computed.In the following, we are going to extend the results ofM16 to a more general estimate of the action-angle vari-ables, through the Torus Mapping method. The resultingDF will be plotted in velocity space by making use of theStäckel fudge. For both transformations, we will use the AGAMA code (Action-based Galaxy Modelling Architecture,Vasiliev 2019, 2018).
Fig. 3.
Top panel: estimate of the spiral arms potential ofSect. 2.5 at the Sun’s position in the plane with 41 complexFourier coefficients in Eq. 10 and the cubic splines reconstruc-tion. The vertical line denotes the true value. Bottom panel:relative accuracy. The accuracy is again well below the percentlevel, although with a slight bias towards higher amplitudes inthe reconstruction w.r.t. to the input spiral potential.
M16 worked in the epicyclic approximation as it providesan analytical expression for evaluating actions and anglesfrom cylindrical coordinates. The Fourier coefficients of theFourier series expansion of the perturbing potentials werethen also determined analytically within this approxima-tion. Approximating the vertical component of the perturb-ing potential by an harmonic oscillator, the 9 Fourier coef-ficients φ jml were then limited to the range j = {− , , } ,corresponding to the θ R oscillations of the potential, and l = {− , , } , corresponding to the θ z oscillations of thepotential close to the Galactic plane.However, the epicyclic approximation is not valid any-more when considering eccentric orbits. Thus, apart fromnearly circular orbits, this approximation is not really ef-ficient to expand the potential in Fourier series. Hereafter,the transformations from angles and actions to positionsand velocities (and reciprocally) are all evaluated numer-ically with AGAMA (Vasiliev 2018). The code makes use ofthe Torus Mapping to go from actions-angles to positions-velocities, and uses the Stäckel fudge for the inverse trans-formation. Our goal now is to obtain the Fourier coefficients
Article number, page 4 of 14l Kazwini et al.: Perturbed distribution functions with accurate action estimates
Fig. 4.
Local stellar velocity distribution at axisymmetric equilibrium from Eq. 16 in the ( u, v ) plane at ( R , z , ϕ ) = ( R , , ).Left panel: epicyclic approximation. Right panel: Stäckel fudge with AGAMA . of a known perturbing potential using these numericallycomputed actions (instead of epicyclic).We proceed in the following way for evaluating Fouriercoefficients of the perturbing potential in Eq. (5). The firststep is to choose a set of actions within a range represent-ing all the orbits of interest in the axisymmetric backgroundconfiguration, each triplet of actions representing one of theorbits. For each orbit, we then define an array of angles ( θ R , θ ϕ , θ z ) . These actions and angles can then all be con-verted to positions thanks to the Torus machinery in AGAMA .For each triplet of actions, a range of positions ( R, ϕ, z ) iscovered by the angles, and we look for the best-fitting co-efficients φ jml ( J R , J z , J ϕ ) , satisfying the following equation(setting t = 0 for the time being): Φ ( R, ϕ, z ) = Re (cid:26)(cid:88) j,l φ jml ( J R , J z , J ϕ )e i( jθ R + mθ ϕ + lθ z ) (cid:27) . (10)This is performed with the method of linear least squares byuse of singular value decomposition as proposed in chapter15.4 of Press et al. (1992). We then interpolate the valueof the coefficient φ jml with cubic splines, also as proposedin Press et al. (1992), chapter 3.3. The number of Fouriercoefficients is chosen to be high enough to ensure that allorbits passing through a given configuration space pointyield the same value of the potential at this point within arelative accuracy of less than 1%.Concretely, we apply this hereafter to the potential of acentral bar and of a 2-armed spiral pattern. The backgroundaxisymmetric potential is chosen to be Model I from Binney& Tremaine (2008). This potential has a bulge describedby a truncated oblate spheroidal power-law, a gaseous discwith a hole at the center, a stellar thin disc and a stellarthick disc, both with a scale-length of 2 kpc, and a dark halowith an oblate two-power-law profile. The Galactocentricradius of the Sun is set at R = 8 kpc , and the local circularvelocity is v = 220 km s − . The potential we choose for the bar is a simple quadrupolepotential (Weinberg 1994; Dehnen 2000) with Φ , b ( R, z, ϕ, t ) = Re (cid:26) Φ a , b ( R, z )e i m ( ϕ − ϕ b − Ω b t ) (cid:27) , (11)where m = 2 , Ω b is the pattern speed of the bar (expressedhereafter in multiples of the angular frequency at the Sun Ω = v /R , where v is the local circular velocity at theGalactocentric radius of the Sun R ), and the azimuth isdefined with respect to a line which would correspond tothe Galactic Center-Sun direction in the Milky Way, ϕ b thus being the angle between the Sun and the long axis ofthe bar. We also choose Φ a , b ( R, z ) = − α b v (cid:18) R R b (cid:19) (cid:18) Rr (cid:19) (cid:18) rR b (cid:19) − R ≥ R b , − (cid:18) rR b (cid:19) R < R b , (12)where r = R + z is the spherical radius, R b is the lengthof the bar and α b represents the maximum ratio betweenthe bar and axisymmetric background radial forces at theSun’s Galactocentric radius R = R . We use hereafter,as a representative example, R b = 0 . R , ϕ b = 25 o and α b = 0 . . We will also consider two typical patternspeeds: Ω b = 1 .
89 Ω and Ω b = 1 .
16 Ω .The bar potential is quite easy to reproduce usingFourier coefficients since it varies smoothly on orbits. Thus,for a study in the Galactic plane, only 41 complex Fouriercoefficients for each triplet of actions are sufficient to ap-proximate the value of the potential with a (much) betterthan 1% accuracy. The Fourier coefficients themselves varysmoothly, as illustrated on Fig. 1, which shows the vari-ations of a few Fourier coefficients as J R and J ϕ increaseseparately, justifying the use of cubic-splines interpolationto get the value of the potential at a specific position. Fig. 2demonstrates the accuracy of our method in reproducingthe bar potential in the Solar neighborhood for differentvalues of the local velocities. The potential is estimated at Article number, page 5 of 14 &A proofs: manuscript no. f1timedep the same configuration space location (the Sun’s position)for the whole range of relevant velocities, with a global accu-racy well below the percent level. This tool is of course notlimited to any specific form of the potential, the only ad-justable parameter being the number of Fourier coefficientsnecessary to recover a given perturbing potential with apercent-level accuracy.
The potential we use for the spiral arms is the following(Cox & Gómez 2002; Monari et al. 2016) Φ , sp ( R, z, ϕ, t ) = Re (cid:110) Φ a , sp ( R, z ) e i m ( ϕ − ϕ sp − Ω sp t ) (cid:111) , (13)where m = 2 , Ω sp is the pattern speed of the spiral arms,and Φ a , sp ( R, z ) = − AR sp KD e i m ln( R/R sp)tan( p ) (cid:20) sech (cid:18) Kzβ (cid:19)(cid:21) β , (14)where K ( R ) = 2 R sin( p ) , β ( R ) = K ( R ) h sp [1 + 0 . K ( R ) h sp ] ,D ( R ) = 1 + K ( R ) h sp + 0 . K ( R ) h sp ] . K ( R ) h sp . (15)Here, R sp = 1 kpc is the length parameter of the logarith-mic spiral potential, h sp = 0 . the height parameter, p = − . o the pitch angle, ϕ sp = − o the phase and A = 683 . s − the amplitude. Hereafter, we will adopta pattern speed Ω sp = 0 .
84 Ω , placing the main resonancesaway from (or at high azimuthal velocities in) the solarneighbourhood.We note that the spiral potential tends to vary morethan the bar potential along a specific orbit, especially whenthe orbit is not restricted to the plane, which requires us toincrease the number of Fourier coefficients to approximate itcorrectly in 3D. Within the Galactic plane, Fig. 3 shows ourreconstruction of the spiral potential at the Sun’s position,again with 41 complex Fourier coefficients. The accuracy isagain below the percent level as in the bar case, althoughin the spiral case there is a (small) bias towards higheramplitudes w.r.t. the input spiral potential. This bias ishowever very small and will not affect our results.
3. Results and comparison with the epicyclicapproximation
From here on, we work with a background axisymmetricDF f as a sum of two quasi-isothermal DFs (Binney &McMillan 2011) for the thin and thick disc: f ( J R , J z , J ϕ ) = f thin + 0 . f thick . (16)The form of each DF is: f ( J R , J z , J ϕ ) = Ω exp( − R g /h R )2 (2 π ) / κ ˜ σ R ˜ σ z z exp (cid:18) − J R κ ˜ σ R − J z ν ˜ σ z (cid:19) (17) where R g , Ω , κ and ν are all functions of J ϕ , and ˜ σ R ( R g ) = ˜ σ R ( R ) exp (cid:18) − R g − R h σ R (cid:19) , ˜ σ z ( R g ) = ˜ σ z ( R ) exp (cid:18) − R g − R h σ z (cid:19) (18)For the thin disc DF f thin , we choose h R = 2 kpc , z =0 . , h σ R = h σ z = 10 kpc , ˜ σ R ( R ) = 35 km s − and ˜ σ z ( R ) = 15 km s − . For the thick disc DF f thick , we choose h R = 2 kpc , z = 1 kpc , h σ R = 10 kpc , h σ z = 5 kpc , ˜ σ R ( R ) = 50 km s − and ˜ σ z ( R ) = 50 km s − .The background axisymmetric potential is chosen to beModel I from Binney & Tremaine (2008), in which the aboveequilbrium DF f is a good representation of the thin andthick disc components. In this model, one has R = 8 kpc and v = 220 km s − .Fig. 4 displays the ( u, v ) -plane in the Solar neighbour-hood within the z = 0 plane for this f axisymmetricbackground, where u = − v R and v = v ϕ − v , obtainedby converting velocity-space into action-space through theepicyclic approximation and the Stäckel fudge from AGAMA .The velocity distributions are globally similar.
In the case of a perturbation with quasi-static amplitudethat has reached its plateau, once the Fourier coefficientsrepresenting the perturbing potential have been computed– from the epicyclic approximation or from Eq. 10 –, theexpression for the perturbed DF can be simply expressedaway from resonances with Eq. 7 as f ( J , θ , t ) = Re (cid:26) n (cid:88) j,l = − n φ jml F jml e i[ jθ R + m ( θ ϕ − Ω p t )+ lθ z ] (cid:27) , (19)with n the order of the Fourier series (in this paper, m = 2 in both the bar and spiral cases), and F jml = j ∂f ∂J R + m ∂f ∂J ϕ + l ∂f ∂J z jω R + m ( ω ϕ − Ω p ) + lω z , (20)where ω R , ω ϕ and ω z can be approximated as epicyclicfrequencies in the epicyclic case, or can be determined with AGAMA . The denominator of F jml may lead to a divergencein the DF when it approaches zero. Following our notationin Eq. 9, let us name it ω s ,jml ( J R , J ϕ , J z ) = jω R + m ( ω ϕ − Ω p ) + lω z . (21)The amount of such resonances is limited in the epicycliccase because, by construction, indices run only over the val-ues j = {− , , } and l = {− , , } , but they can be muchmore numerous in the more accurate AGAMA case. For thebar potential of Eq. 11 and Eq. 12, and choosing a patternspeed Ω p = 1 . like for our fiducial bar model, we ex-plore in Fig. 5 and Fig. 6, the values of ω s ,jml ( J R , J ϕ , J z ) inaction space, when varying the pair of integer indices ( j , l ).The actions are renormalized by the radial velocity disper-sion of the thin disc, circular velocity, and vertical velocitydispersion of the thin disc at the Sun, respectively, to onlydisplay a relevant range of actions. Exploring indices in the Article number, page 6 of 14l Kazwini et al.: Perturbed distribution functions with accurate action estimates
Fig. 5.
Values of log( ω s ,jml ) in the ( J R , J ϕ ) plane with fixed J z = 10 kpc km s − , for a few examples of combinations of ( j, l ) indices giving rise to resonant zones in action space (remember that m = 2 ). The pattern speed Ω p is here .
89 Ω . The two actionsare renormalized by the radial velocity dispersion of the thin disc and circular velocity at the Sun, respectively. The deep bluelines correspond to resonant zones. For instance, the (1 , case corrsponds to the traditional OLR (for a non-zero J z ). Most otherlow-order combinations of indices did not give rise to any relevant resonant zone in the region of interest. Fig. 6.
Values of log( ω s ,jml ) in the ( J R , J z ) plane with fixed J ϕ = 1759 kpc km s − . The pattern speed Ω p is the one of our fiducialcentral bar fixed at .
89 Ω . The two actions are renormalized by the radial velocity dispersion and vertical velocity dispersion ofthe thin disc at the Sun, respectively. The deep blue lines correspond to resonance zones. Most combinations of indices exploreddid not give rise to any relevant resonant zone in the region of interest. range [ − , +4] , it appears clearly that most combinationsdo not induce a resonance which is relevant to the dynam-ics of the Solar neighbourhood. We only display on Fig. 5and Fig. 6 the combination of indices for which a resonantzone appears in the plotted region of action space. It ap-pears clearly that very few low order resonances are indeedpresent in the range of actions that are really relevant forthe solar neighbourhood. As displayed on Fig. 7 and Fig. 8,the same is true for a lower pattern speed Ω p = 0 .
84 Ω ,corresponding to the pattern speed of our fiducial spiralpotential.While a specific treatment is needed in these resonantzones (e.g., Monari et al. 2017a), their sparse presence in therelevant zones of action space guarantees us to get meaning-ful results. Note also that, while the form of the DF is notwell estimated in the resonant zones, the signature of the resonances (and thus their location in velocity space) canclearly be identified with this linear perturbation method. We are now in a position to compare the linear deforma-tion of local velocity space for different action estimates,namely the epicyclic case used in previous works and themore accurate
AGAMA action estimates.Fig. 9 displays the f + f linearly perturbed distribu-tion function at the position of the Sun in the Galactic planefor the bar potential of Sect. 2.4 and two different patternspeeds, as well as for the spiral potential of Sect. 2.5. Asin Monari et al. (2017b), whenever f > f , we cap f atthe value of f to roughly represent the resonant zone. Themore rigorous approach, which we leave to further works in Article number, page 7 of 14 &A proofs: manuscript no. f1timedep
Fig. 7.
Same as Fig. 5, with the combinations of indices giving rise to resonant zones for Ω p = 0 .
84 Ω . Fig. 8.
Same as Fig. 6, with some combinations of indices giving rise to resonant zones for Ω p = 0 .
84 Ω . Fig. 9.
The DF of Fig. 4 in velocity space at the Solar position within the Galactic plane, now perturbed to linear order by abar (perturbing potential of Sect. 2.4) with pattern speeds Ω b = 1 .
16 Ω (left) and Ω b = 1 . (middle), or by a spiral pattern(perturbing potential of Sect. 2.5) with pattern speed Ω sp = 0 .
84 Ω (right). Top row: Epicyclic approximation. Bottom row:Stäckel fudge. the context of AGAMA actions (Al Kazwini et al., in prep.),is to treat the DF with the method of Monari et al. (2017a)in these regions. However, while the DF within the reso-nant zone is not well modelled by the present method, the location of resonances is accurately reproduced. The lineardeformation outside of the resonant zones is very accuratelydescribed by our method as well. Interestingly enough, the linear deformation due to the bar is generally stronger inthe
AGAMA case, and that due to the spiral is weaker in the
AGAMA case. This means that reproducing the effect of spi-ral arms on the local velocity distribution might require ahigher amplitude when considering an accurate estimate ofthe action-angle variables rather than the epicyclic approx-imation.
Article number, page 8 of 14l Kazwini et al.: Perturbed distribution functions with accurate action estimates
Fig. 10.
Local stellar velocity distribution perturbed to linear order at the Solar Galactocentric radius and azimuth at threedifferent heights (left: z = 0 kpc , middle: z = 0 . , right: z = 1 kpc ), when perturbed by a bar (perturbing potential ofSect. 2.4) with pattern speed Ω b = 1 . . Top row: Epicyclic approximation. Bottom row: Stäckel fudge. Fig. 11.
Same as Fig. 10, in the Stäckel fudge case, but now for joint perturbation by a bar (perturbing potential of Sect. 2.4)with pattern speed Ω b = 1 . and a spiral pattern (perturbing potential of Sect. 2.5) with pattern speed Ω sp = 0 .
84 Ω . The case of the pattern speed of the bar being .
89 Ω would correspond to a configuration where the ‘Hercules’stream at negative u and negative v corresponds to thelinear deformation below the outer Lindblad reso-nance of the bar (e.g., Dehnen 2000; Minchev et al. 2007;Monari et al. 2017b; Fragkoudi et al. 2019). Interestingly,this feature is less squashed in the more realistic AGAMA case. Moreover, a resonance unnoticed within the epicyclicapproximation appears at high azimuthal velocities: we canidentify this resonance as the outer resonance of thebar (Dehnen 2000). In the spiral case, the resonant ridge atlarge azimuthal velocities can be identified as the corotationof the spiral pattern.Fig. 10 displays the linear deformation due to the bar,in the case of the pattern speed being .
89 Ω , at differentheights above the Galactic plane, both in the epicyclic and AGAMA cases. As can be seen on this figure, the epicyclicapproximation quickly becomes very imprecise at largeheights, because it implies a drastic falloff of the densitywith height while not changing the azimuthal velocity dis-tribution due to the hypothesis of complete decoupling ofvertical motions. In the
AGAMA case, the azimuthal velocity distribution is affected by a larger asymmetric drift at largeheights, and the location of the outer Lindblad resonanceof the bar in the uv -plane is also displaced to lower az-imuthal v at larger heights. This is because at fixed J ϕ , theazimuthal and radial frequencies computed with AGAMA arelower at higher z , meaning that one needs to reach lower J ϕ (corresponding to orbits whose guiding radii are in theinner Galaxy) to reach the resonance. This trend is alreadyvisible at z = 0 . kpc where the ‘Hercules’ feature movesto lower azimuthal v in the AGAMA case, and becomes dra-matic at z = 1 kpc, where the epicyclic approximation isnot representing accurately the velocity distribution at all.Interestingly, comparing the displacement with height ofthe OLR in the case of a bar with pattern speed .
89 Ω with that of the corotation in the case of a .
16 Ω pat-tern speed, we noted that the corotation location in the uv -plane is more displaced than the OLR. This is becausethe corotation only depends on the azimuthal frequency,while the OLR depends on a combination of the azimuthaland radial frequencies. However, this different behaviourclearly depends on the location in the velocity plane: athigh-azimuthal velocities, the corotation of the spiral arm Article number, page 9 of 14 &A proofs: manuscript no. f1timedep is less displaced with z than the outer resonance ofthe bar, causing them to almost merge at z = 0 . when linearly adding the effect of the bar and spiral onFig. 11. The two resonances are however separated againat z = 1 kpc . At this height, other resonances involvingthe radial and vertical frequencies also appear in velocityspace. All these subtle variations however depend stronglyon the background Galactic potential. This means that,once identifying the resonances potentially responsible formoving groups in the Solar neighbourhood, studying theirposition in the uv -plane as a function of z can in princi-ple be a powerful way to constrain the 3D structure of theGalactic potential. This cannot be done within the epicyclicapproximation.
4. Adding the temporal evolution
Up to now, we always considered a constant amplitude forthe perturbing potential in order to determine an analyt-ical expression for the perturbed DF. In this section, weinvestigate the time dependence of the DF by choosing atime-varying amplitude for the perturbing potential.
The expression we will use for the time dependence function g controlling the amplitude of the perturbation during itsgrowth is g ( t ) = 1 − cos ( πt/t f )2 , (22)where t f is the time at which the perturbation is completelyformed, expressed in Gyr. We will consider t f = 0 . .The motivation for this choice of growth function is itsanalytic simplicity, having a function starting from exactlyzero at the origin, and smooth over the whole consideredrange. The first derivative, [ π/ (2 t f )] sin( πt/t f ) , assures thecontinuity at and t f with both stages, fixed at for t ≤ and at for t ≥ t f (the first derivative is thus equal to at and t f ). We now take the integral of Eq. (6), restricted to [0 , t ] (be-cause the g function is equal to on ] −∞ , ) and inte-grate by parts. We take φ n ( J (cid:48) , t (cid:48) ) = g ( t (cid:48) ) h ( t (cid:48) ) φ n ( J ) , with h ( t (cid:48) ) = e − i m Ω p t (cid:48) , and we define η ( t ) ≡ e i θ s , n ( t ) i ω s , n → d η = e i θ s , n ( t ) d t, (23)allowing us to rewrite Eq. (6) as f ( J , θ , t ) = Re (cid:26) i ∂f ∂ J ( J ) · (cid:88) n n φ n ( J ) (cid:90) t g ( t (cid:48) ) d η d t (cid:48) ( t (cid:48) ) d t (cid:48) (cid:27) . (24)We can now integrate by parts (cid:90) t g ( t (cid:48) ) d η d t (cid:48) ( t (cid:48) ) d t (cid:48) = [ g ( t (cid:48) ) η ( t (cid:48) )] t − (cid:90) t d g d t (cid:48) ( t (cid:48) ) η ( t (cid:48) ) d t (cid:48) . (25)Since g (0) = 0 , [ g ( t (cid:48) ) η ( t (cid:48) )] t = g ( t ) η ( t ) . (26) To calculate the second part of the integral, since d g ( t ) / d t = π/ (2 t f ) sin( πt/t f ) , we write, (cid:90) t d g d t (cid:48) ( t (cid:48) ) η ( t (cid:48) )d t (cid:48) = π t f ω s , n (cid:90) t sin (cid:18) πt (cid:48) t f (cid:19) e i θ s , n ( t (cid:48) ) d t (cid:48) . (27)We look for a primitive G of sin( πt/t f )e i θ s , n ( t ) of the form G ( t ) = (cid:20) A cos (cid:18) πtt f (cid:19) + B sin (cid:18) πtt f (cid:19)(cid:21) e i θ s , n ( t ) . (28)Deriving G ( t ) with respect to t , and equating it to the in-tegrand in Eq. 27 we get B πt f + A i ω s , n = 0 and B i ω s , n − A πt f = 1 , (29)which leads to A = π/t f ω s , n − ( π/t f ) and B = − i ω s , n ω s , n − ( π/t f ) . (30)Substituting the Eqs. (26) and (28) into Eq. (25) resultsin the following expression for the perturbed DF f ( J , θ , t ) = Re (cid:40) ∂f ∂ J ( J ) · (cid:88) n n φ n ( J ) × (cid:34) (cid:18) − cos (cid:18) πtt f (cid:19)(cid:19) e i θ s , n ω s , n − π t f ω s , n ω s , n − ( π/t f ) × (cid:18)(cid:18) πt f cos (cid:18) πtt f (cid:19) − i ω s , n sin (cid:18) πtt f (cid:19)(cid:19) e i θ s , n − πt f e i( θ s , n − ω s , n t ) (cid:19) (cid:35) (cid:41) . (31)Note that we do not exactly recover the static case at t = t f because all derivatives of g ( t ) are not strictly zero atthe initial and final time, as assumed in M16. If reachinga true plateau after t f in an analytic fashion, one wouldnevertheless converge towards the static case very quickly.Now we can study analytically how the linear responseto a fiducial bar with Ω b = 1 . evolves with time. Asbefore, the method is not strictly valid at resonances, wherethe treatment of Monari et al. (2017a) must be applied. Itis nevertheless interesting to see on Fig. 12 how the lineardeformation of the velocity plane evolves with time nearresonances (in a patch co-moving with the bar, hence atconstant azimuthal angle to the bar) while the amplitudeof the perturbation grows. The effect of the OLR appearsas soon as the perturbation starts to grow. As it progres-sively grows, the two linear modes in the DF separate, andlead to a velocity plane already very much resembling thestationary form of the perturbed DF after .
25 Gyr , thatis when g ( t ) = 0 . and the perturbation is half-formed. Inthe absence of a pattern speed variation, it is therefore notnecessarily obvious to disentangle the effect of a bar whoseamplitude is growing from that of a fully formed bar withlarger and constant amplitude. Article number, page 10 of 14l Kazwini et al.: Perturbed distribution functions with accurate action estimates
Fig. 12.
Local stellar velocity distribution perturbed to linear order in the Galactic plane by the bar of Sect. 2.4 with Ω b = 1 . with the Stäckel fudge, and an amplitude of the bar growing as described in Sect. 4.1. The first eleven panels display the temporalevolution of the perturbation. The last panel displays the stationary case. The amplitude of the bar goes from at t = 0 to itsplateau ( g ( t ) = 1 ) at t = 0 . .
5. Conclusion
Starting from the formalism exposed in M16, we proposed amore accurate way to determine the DF of the Galactic discperturbed to linear order by a non-axisymmetric perturba-tion, using a more accurate action-angle coordinate system.First, we used the Torus Mapping from
AGAMA to numeri-cally compute the perturbing potential in action-angle co-ordinates, as a Fourier series expansion over the angles. Weshowed that we could estimate typical non-axisymmetricperturbing potentials with an accuracy well below the per-cent level. The algorithm can be applied to any perturbingpotential, including non-plane symmetric vertical perturba-tions, which will be particularly important when studyingthe vertical perturbations of the disc with similar methods(Rozier et al., in prep.).We then compute the DF perturbed to linear order bya typical bar or spiral potential (or a linear combinationof both), and compute the local stellar velocity distribu-tion by converting velocities to actions and angles throughthe Stäckel fudge implemented in
AGAMA . The results arecompared to those obtained by using the epicyclic approx-imation. The linear deformation due to the bar is gener-ally stronger in the
AGAMA case, and that due to the spi-ral is weaker in the
AGAMA case.This means that reproduc-ing the effect of spiral arms on the local velocity distri-bution might require a higher amplitude when consideringan accurate estimate of the action-angle variables ratherthan the epicyclic approximation. Most importantly, theepicyclic approximation is inadequate at large heights abovethe plane, because it implies a drastic falloff of the densitywith height while not changing the azimuthal velocity dis-tribution due to the hypothesis of complete decoupling of vertical motions. In the
AGAMA case instead, the locationsof resonances are displaced to lower azimuthal v at largerheights. This is already visible at z = 0 . but becomesa dramatic effect at z = 1 kpc . Thus, the position of movinggroups in the uv -plane as a function of z can be a powerfulway to constrain the 3D structure of the Galactic potential.The key to exploring this will be the DR3 of Gaia (Brown2019) with its ∼ millions radial velocities allowing tobetter probe the z -axis above and below the Milky Wayplane.Finally, the temporal treatment is also an improvementover M16. We applied it to the case of a bar of growingamplitude, with an analytic evolution of the amplitude. Asthe bar progressively grows, the two linear modes in the DFseparate, and lead to a velocity plane already very muchresembling the stationary form of the perturbed DF oncethe perturbation is half-formed. In the absence of a patternspeed variation, it is therefore not necessarily obvious todisentangle the effect of a bar whose amplitude is growingfrom that of a fully formed bar with larger and constantamplitude. Note that we explored here a peculiar form ofthe growth function motivated by its analytic simplicity. Ifthe perturbation grows by a linear instability, an exponen-tial growth would be more realistic. Numerical experimentsare usually well fitted by a logistic function (exponentialgrowth at the beginning, and saturation to the plateau).One problem for our treatment is that the logistic func-tion is never strictly equal to 0. Besides, there is hope thatsimilar analytical simplifications as those for the amplitudegrowth studied here can also be made with this function,which we will investigate in the future.While the form of the DF is not well estimated in theresonant zones with the linear perturbations presented in Article number, page 11 of 14 &A proofs: manuscript no. f1timedep this paper, the signature of the resonances (and thus theirlocation in velocity space) can clearly be identified with thislinear perturbation method. The more rigorous approach,which we leave to further works in the context of
AGAMA ac-tions (Al Kazwini et al., in prep.), is to treat the DF withthe method of Monari et al. (2017a) in these regions, patch-ing these results over the linear deformation computed here.Another caveat is that the Torus Mapping was used for ex-pressing the perturbing potential in actions and angles, butfor the estimate of the local stellar velocity field, we madeuse of the less precise Stäckel fudge method. Therefore, an-other promising way for improvement would be to use thenew
ACTIONFINDER deep-learning algorithm (Ibata et al.2020) to make the reverse transformation.The tools presented in this paper will be useful for athorough analytical dynamical modelling of the complexvelocity distribution of Milky Way disc stars as measured bypast and upcoming Gaia data releases. These tools will alsobe useful for fully self-consistent treatments of the responseof the disc to external perturbations. The ultimate goalwill be to adjust models to the exquisite data from Gaia,which cannot be done properly with N-body simulationsdue to the vast parameter space to explore. The theoreticaltools and the new code presented in this paper thereforerepresent a useful step in this direction.
Acknowledgements.
AS, BF, GM, SR, PR and RI acknowledge fund-ing from the Agence Nationale de la Recherche (ANR project ANR-18-CE31-0006 and ANR-19-CE31-0017) and from the European ResearchCouncil (ERC) under the European Union’s Horizon 2020 researchand innovation programme (grant agreement No. 834148).
References
Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561,360Binney, J. 2012, MNRAS, 426, 1324Binney, J. 2020a, MNRAS, 495, 886Binney, J. 2020b, MNRAS, 495, 895Binney, J. & McMillan, P. 2011, MNRAS, 413, 1889Binney, J. & McMillan, P. J. 2016, MNRAS, 456, 1982Binney, J. & Piffl, T. 2015, MNRAS, 454, 3653Binney, J. & Schönrich, R. 2018, MNRAS, 481, 1501Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition(Princeton University Press)Bland-Hawthorn, J. & Tepper-Garcia, T. 2020, arXiv e-prints,arXiv:2009.02434Brown, A. G. A. 2019, in The Gaia Universe, 18Chiba, R., Friske, J. K. S., & Schönrich, R. 2020, MN-RAS[ arXiv:1912.04304 ]Cole, D. R. & Binney, J. 2017, MNRAS, 465, 798Cox, D. P. & Gómez, G. C. 2002, ApJS, 142, 261Dehnen, W. 1998, AJ, 115, 2384Dehnen, W. 2000, AJ, 119, 800Famaey, B. & Dejonghe, H. 2003, MNRAS, 340, 752Famaey, B., Jorissen, A., Luri, X., et al. 2005, A&A, 430, 165Fragkoudi, F., Katz, D., Trick, W., et al. 2019, MNRAS, 488, 3324Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018a, A&A,616, A1Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2020, arXive-prints, arXiv:2012.01533Gaia Collaboration, Katz, D., Antoja, T., et al. 2018b, A&A, 616,A11Hilmi, T., Minchev, I., Buck, T., et al. 2020, MNRAS, 497, 933Ibata, R., Diakogiannis, F., Famaey, B., & Monari, G. 2020, arXive-prints, arXiv:2012.05250Laporte, C. F. P., Famaey, B., Monari, G., et al. 2020, A&A, 643, L3Laporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, F. A. 2019,MNRAS, 485, 3134McGill, C. & Binney, J. 1990, MNRAS, 244, 634Minchev, I., Nordhaus, J., & Quillen, A. C. 2007, ApJ, 664, L31 Monari, G., Famaey, B., Fouvry, J.-B., & Binney, J. 2017a, MNRAS,471, 4314Monari, G., Famaey, B., & Siebert, A. 2016, MNRAS, 457, 2569Monari, G., Famaey, B., Siebert, A., et al. 2019a, A&A, 632, A107Monari, G., Famaey, B., Siebert, A., et al. 2017b, MNRAS, 465, 1443Monari, G., Famaey, B., Siebert, A., Wegg, C., & Gerhard, O. 2019b,A&A, 626, A41Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P.1992, Numerical Recipes in C: Second Edition (Cambridge Univer-sity Press)Ramos, P., Antoja, T., & Figueras, F. 2018, A&A, 619, A72Sanders, J. L. & Binney, J. 2016, MNRAS, 457, 2107Sellwood, J. A. & Carlberg, R. G. 2014, ApJ, 785, 137Trick, W. H., Coronado, J., & Rix, H.-W. 2019, MNRAS, 484, 3291Vasiliev, E. 2018, AGAMA: Action-based galaxy modeling framework,Astrophysics Source Code LibraryVasiliev, E. 2019, MNRAS, 482, 1525Weinberg, M. D. 1994, ApJ, 420, 597
Article number, page 12 of 14l Kazwini et al.: Perturbed distribution functions with accurate action estimates
Appendix A: PERDIGAL code presentation
All the results shown in this article were obtained with asingle code, written in C++ and Python, that calculatesboth the perturbing potential in action-angle coordinatesand the perturbed DF. The code is named
PERDIGAL (PER-turbed DIstribution functions for the GALactic disc), andwill be made available on request, although it may even-tually be embedded in a larger Galactic dynamics toolkit.Launching the code without argument gives the explana-tions shown in Fig. A.1.The calculation of the perturbed DF consists in five steps:the creation of a file storing the positions of the DF (DF-pos), the creation (via
AGAMA ) of several files containingpositions from a lot of orbits (FCorb) required for the fol-lowing step, the determination of Fourier coefficients forthe perturbing potential to create the grid (FCgrid), thedetermination of the Fourier coefficients for each positionwhich the DF will be calculated at (FCforDF) and finallythe calculation of the perturbed DF (PertDF). The modeALL is processing the three above steps in one time, with-out saving Fourier coefficients from FCgrid and FCforDF infiles. Then, the mode ALL regroups in one the three modes:FCgrid, FCforDF and PertDF. However, this mode shouldnot be used without certainty about the decomposition ofthe potential. Indeed, this process depends a lot on the ini-tial conditions imposed to the code, and a check after eachpart of the calculation is recommended. Finally, PDFdispdisplays the perturbed DF using Matplotlib.All useful parameters are stored in two particular filesand can be modified without compiling the code. They con-tain parameters for the calculation of Fourier coefficients,parameters for the perturbing potentials and others forwhether the DF is determined at first or second order, orwith a time-dependent amplitude for the perturbing poten-tial. Figure A.2 is a logigram reminding the procedure ofthe calculation of the perturbed DF.
Article number, page 13 of 14 &A proofs: manuscript no. f1timedep
Fig. A.1.
Guide for the use of the code. The first three modes deal with the determination of the Fourier coefficients whenexpanding in series the perturbing potential (FCgrid and FCforDF), and the calculation of the perturbed distribution function(PertDF). PDFdisp allows this function to be displayed as a distribution of velocities. ALL regroups the three above modes in one,but is not recommended because of the need to check the Fourier decomposition. EPI calculates directly the perturbed distributionfunction under the epicyclic approximation, as explained in M16. The FUNCTIONs correspond to different operations requiredbefore using the MODEs. DFpos creates the file storing the positions at which the distribution function is going to be determined.FCorb creates several files containing positions from a lot of orbits, needed to calculate the Fourier coefficients in the FCgrid step.