Perturbed distribution functions with accurate action estimates for the Galactic disc

In the Gaia era, understanding the effects of the perturbations of the Galactic disc is of major importance in the context of dynamical modelling. In this theoretical paper we extend previous work in which, making use of the epicyclic approximation, the linearized Boltzmann equation had been used to explicitly compute, away from resonances, the perturbed distribution function of a Galactic thin-disc population in the presence of a non-axisymmetric perturbation of 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 coordinates with a typical accuracy at the per cent level. The use of more accurate action estimates allows us to identify resonances such as the outer 1:1 bar resonance at higher azimuthal velocities than the outer Lindblad resonance (OLR), and to extend our previous theoretical results well above the Galactic plane, where we explicitly show how they differ from the epicyclic approximation. In particular, the displacement of resonances in velocity space as a function of height can in principle constrain the 3D structure of the Galactic potential. Second, we allow the perturbation to be time dependent, thereby allowing us to model the effect of transient spiral arms or a growing bar. The theoretical framework and tools presented here will be useful for a thorough analytical dynamical modelling of the complex velocity distribution of disc stars as measured by past and upcoming Gaia data releases.


Introduction
The natural canonical coordinate system of phase-space for Galactic dynamics and perturbation theory is the system of action-angle variables (Binney & Tremaine 2008). In an axisymmetric system in equilibrium, the Jeans theorem implies that the phase-space distribution function (DF) of any stellar (or dark matter) component can be expressed solely as a function 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., the Galactic bar and spiral arms) must be included in this process, together with the response of the distribution function. Within the resonant regions, to fully capture the behaviour of the DF, one needs to construct for each perturber new orbital tori, complete with a new system of action-angle variables (e.g., Monari et al. 2017a;Binney 2020a,b). Away from resonances, however, one can simply use the linearized Boltzmann equation. This also allows us to accurately identify the location of resonances. This is particularly important in the context of the interpretation of recent data from the Gaia mission (Gaia Collaboration 2018a, 2021, which revealed in exquisite detail the fine structure of stellar action space (e.g., Trick et al. 2019;Monari et al. 2019a,b). While the existence of moving groups of dynami-cal origin had been known for a long time in local velocity space around the Sun (e.g., Dehnen 1998;Famaey et al. 2005), Gaia revealed their structure in exquisite detail (Ramos et al. 2018) and also provided an estimate of their age distribution (Laporte et al. 2020), together with the shape of the global velocity field away from the Sun within the Galactic disc (Gaia Collaboration 2018b). One additional major finding of Gaia is the existence of a local phase-spiral in vertical height versus vertical velocity in the solar neighbourhood , which might be related to a vertical perturbation of the disc, for example by the Sagittarius dwarf galaxy (e.g., Laporte et al. 2019;Binney & Schönrich 2018;Bland-Hawthorn & Tepper-García 2021).
In previous theoretical work, Monari et al. (2016) (hereafter M16) explicitly computed the response of an axisymmetric DF in action space, representing a Milky Way thin-disc stellar population, to a quasi-stationary three-dimensional spiral potential, using the epicyclic approximation to model the actions, which is only a valid approximation for quasi-circular orbits in the thin disc. It was notably shown that the first-order moments of the perturbed DF then give rise to non-zero radial and vertical bulk flows (breathing modes). However, to treat perturbations away from the mid-plane, which is particularly important in the Gaia context, one cannot make use of the epicyclic approximation 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 different during the period of rising or declining amplitude. The same can be true for the bar, whose amplitude and pattern speed can also grow or vary with time (e.g., Chiba et al. 2021;Hilmi et al. 2020).
Here we improve on this previous modelling of M16 in two ways. First, we use a better estimate than the epicyclic approximation for the action-angle variables, relying on the torus mapping method of McGill & Binney (1990) to convert from actions and angles to positions and velocities, and on the Stäckel fudge (Binney 2012;Sanders & Binney 2016) for the reverse transformation. This will allow us to extend previous results to eccentric orbits and orbits wandering well above the Galactic plane. The routines developed and presented in this paper will also be of fundamental importance to study the vertical perturbations of the Galactic disc in further works. Second, we let the perturbation evolve with time, thereby allowing us to model the effect of a growing bar.
The paper is organized as follows. In Sect. 2, after a short reminder of the theoretical framework of perturbed DF, already given in detail in M16, we present the method used to expand in Fourier series the perturbing potential within the new actionangle coordinate system. Then a comparison with the results in the epicyclic case is made in Sect. 3. In Sect. 4 we explore the temporal treatment of the DF for an analytic growth of the amplitude of the perturber. Finally, we conclude and discuss the possible future applications of these theoretical tools in Sect. 5. The Appendix is dedicated to the presentation of the code.

Action-angle variables
The canonical phase-space action-angle variables (J,θ) of an integrable system are obtained from a canonical transformation implicitly using Hamilton's characteristic function as a type 2 generating function. The actions J are defined as new generalized momenta corresponding to a closed path integral of the velocities along their corresponding canonically conjugate position variable, namely J i = v i dx i /(2π). Since this does not depend on time, these actions are integrals of motion, and the Hamiltonian can be expressed purely as a function of these actions.
It turns out that Galactic potentials are close enough to integrable systems that actions can be estimated for them. For quasicircular orbits close to the Galactic plane, with separable motion in the vertical and horizontal directions, one can locally approximate the radial and vertical motions of an orbit with harmonic motions, which is known as the epicyclic approximation. The radial and vertical actions then simply correspond to the radial and vertical energies divided by their respective (radial and vertical) epicyclic frequency. However, the epicyclic approximation is no longer valid when considering orbits with higher eccentricity, or with a large vertical amplitude. More precise ways of determining the action and angle coordinates have been devised. They typically differ depending on whether one wishes to transform angles and actions to positions and velocities or instead to make the reverse transformation from positions and velocities to actions and angles. In the first case a very efficient method is the torus mapping method first introduced by McGill & Binney (1990) (see also Binney & McMillan 2016 for a recent overview), while in the second case a Stäckel fudge is generally used (Binney 2012;Sanders & Binney 2016).
The general idea of torus mapping is to first express the Hamiltonian in the action-angle coordinates (J T , θ T ) of a toy potential, for which the transformation to positions and velocities is fully known analytically. The algorithm then searches for a type 2 generating function G(θ T , J) expressed as a Fourier series expansion on the toy angles θ T , for which the Fourier coefficients are such that the Hamiltonian remains constant at constant J. This generating function fully defines the canonical transformation from actions and angles to positions and velocities. For the inverse transformation, an estimate based on separable potentials can be used. These potentials are known as Stäckel potentials (e.g., Famaey & Dejonghe 2003), for which each generalized momentum depends on its conjugated position through three isolating integrals of the motion. These potentials are expressed in spheroidal coordinates defined by a focal distance. For a Stäckel potential, this focal distance of the coordinate system is related to the first and second derivatives of the potential. One can thus use the true potential at any configuration space point to compute the equivalent focal distance as if the potential were of Stäckel form, and from there compute the corresponding isolating integrals of the motion and the actions. In the following we make use of both types of transformations, namely the torus mapping to express the potential in action-angle coordinates and the Stäckel fudge to represent distribution functions in velocity space at a given configuration space point.
We now let H 0 (J) be the Hamiltonian of the axisymmetric and time-independent zeroth order gravitational potential Φ 0 of the Galaxy. The equations of motion connecting actions J and the canonically conjugate angles θ are then simplẏ with ω being the fundamental orbital frequencies. Thus, for a given orbit, the actions J are constant in time, defining an orbital torus on which the angles θ evolve linearly with time, according to θ(t) = θ 0 + ωt. The Jeans theorem then tells us that the phasespace distribution function (DF) of an axisymmetric potential f = f 0 (J) is in equilibrium. In other words, f 0 is a solution of the collisionless Boltzmann equation:

Perturbed distribution functions
We now let Φ 1 be the potential of a small perturbation to the axisymmetric background potential Φ 0 of the Galaxy, with an amplitude relative to this axisymmetric background 1. Then the total potential is Φ = Φ 0 + Φ 1 and the DF becomes, to first order in , f = f 0 + f 1 , which is still a solution of the collisionless Boltzmann equation. Inserting f = f 0 + f 1 in Eq. (2), and keeping only the terms of order , leads to the linearized collisionless Boltzmann equation where [,] is the Poisson bracket. Integrating Eq.
(3) over time, from −∞ to the time t, leads to where J and θ correspond to the actions and angles in the unperturbed case as a function of time (i.e. constant actions J and angles evolving linearly). Since any function of the angles is 2π-periodic in the angles, the perturbing potential Φ 1 can be expanded in a Fourier series as Hereafter we consider in-plane perturbations such as spiral arms, meaning that we can write the time-varying Fourier coefficients in a non-rotating frame as φ n (J, t) = g(t) h(t) φ n (J), where g(t) controls the amplitude of the perturbation and h(t) controls its pattern speed, with h(t) = e −imΩ p t , where Ω p is the pattern speed of the perturbation and m its azimuthal wave number (i.e. its multiplicity). Hereafter we mainly consider m = 2 perturbations. The vector index n is a triplet of scalar integer indices ( j, k, l) running in principle from −∞ to ∞, but in practice limited to a given range sufficient to approximate the perturbing potential.
In the case of an m-fold in-plane perturbation, it is sufficient to take k = m. The main goal of this section is to express typical non-axisymmetric perturbing potentials originally expressed in configuration space as such a Fourier series in action-angle space. The algorithm that we present in Sect. 2.3 can be applied to any perturbing potential, however, including non-plane symmetric vertical perturbations. Once the potential is expressed as a function of angles and actions as in Eq. (5), then Eq. (4) becomes e −imΩ p t , and the amplitude of the perturbing potential constant in time at present time (g(t) = 1), and zero and constant in time at −∞, led to the following explicit solution for f 1 (J, θ, t), where we defined The subscript 's' stands for slow, because in the proximity of a resonance of the type ω s,n = 0, the angle θ s,n evolves very slowly. One can also immediately see that the linearized solution above is valid only away from such resonances since it diverges for these orbits. Orbits near these resonances are actually trapped, and for them the determination of the linearly perturbed DF becomes inappropriate. Specific treatment for these resonant regions is required, which was addressed in Monari et al. (2017a) within the epicyclic approximation, and by Binney (2020a) in a more general context. Using the epicyclic approximation the Fourier coefficients of a spiral potential have been computed analytically in M16 with indices n = ( j, k, l) running over the values j = {−1, 0, 1}, k = m = 2, and l = {−2, 0, 2}, and the perturbed distribution function away from resonances was then computed. In the following we extend the results of M16 to a more general estimate of the action-angle variables through the torus mapping method. The resulting DF is plotted in velocity space by making use of the Stäckel fudge. For both transformations we use the Actionbased Galaxy Modelling Architecture (AGAMA; Vasiliev 2019, 2018).

Perturbing potential in actions and angles
In previous work, M16 worked in the epicyclic approximation as it provides an analytical expression for evaluating actions and angles from cylindrical coordinates. The Fourier coefficients of the Fourier series expansion of the perturbing potentials were then also determined analytically within this approximation. Approximating the vertical component of the perturbing potential by a harmonic oscillator, the nine Fourier coefficients φ jml were then limited to the range j = {−1, 0, 1}, corresponding to the θ R oscillations of the potential, and l = {−2, 0, 2}, corresponding to the θ z oscillations of the potential close to the Galactic plane.
However, the epicyclic approximation is only valid for nearly circular orbits, but not when considering eccentric orbits. Hereafter the transformations from angles and actions to positions and velocities (and reciprocally) are instead evaluated numerically with AGAMA (Vasiliev 2018). The code makes use of torus mapping to go from action-angle to position-velocity, and uses the Stäckel fudge for the inverse transformation. Our goal now is to obtain the Fourier coefficients of a known perturbing potential using these numerically computed actions (instead of epicyclic).
We proceed in the following way to evaluate Fourier coefficients of the perturbing potential in Eq. (5). The first step is to choose a set of actions within a range representing all the orbits of interest in the axisymmetric background configuration, each triplet of actions representing one of the orbits. For instance, in the solar neighbourhood we consider radial actions ranging from 0 to 220 kpc km s −1 , azimuthal actions ranging from 1200 to 2160 kpc km s −1 , and vertical actions ranging from 0 to 26 kpc km s −1 depending on the height above the Galactic plane. For each orbit, we then define an array of angles (θ R , θ ϕ , θ z ). These actions and angles can then all be converted to positions thanks to the torus machinery in AGAMA. For each triplet of actions, a range of positions (R, ϕ, z) is covered by the angles, and we look for the best-fitting coefficients φ jml (J R , J z , J ϕ ), satisfying the following equation (setting t = 0 for the time being): This is performed with the method of linear least squares using singular value decomposition, as proposed in chapter 15.4 of Press et al. (1992). We then interpolate the value of the coefficient φ jml with cubic splines, also proposed in Press et al. (1992), chapter 3.3. The number of Fourier coefficients is chosen to be high enough to ensure that all orbits passing through a given configuration space point yield the same value of the potential at this point within a relative accuracy of less than 1%. Concretely, we apply this hereafter to the potential of a central bar and of a two-armed spiral pattern. The background axisymmetric potential is chosen to be Model I from Binney & Tremaine (2008). This potential has a bulge described by a truncated oblate spheroidal power law; a gaseous disc with a hole at the centre; a stellar thin disc and a stellar thick disc, both with a scale-length of 2 kpc; and a dark halo with an oblate twopower-law profile. The galactocentric radius of the Sun is set at R 0 = 8 kpc, and the local circular velocity is v 0 = 220 km s −1 .

Bar potential
The potential we choose for the bar is a simple quadrupole potential (Weinberg 1994;Dehnen 2000) with where m = 2, Ω b is the pattern speed of the bar (expressed hereafter in multiples of the angular frequency at the Sun where v 0 is the local circular velocity at the galactocentric radius of the Sun R 0 ), and the azimuth is defined with respect to a line corresponding to the Galactic centre-Sun direction in the Milky Way, ϕ b thus being the angle between the Sun and the long axis of the bar. We also choose where r 2 = R 2 + z 2 is the spherical radius, R b is the length of the bar, and α b represents the maximum ratio between the bar and axisymmetric background radial forces at the Sun's galactocentric radius R = R 0 . We use hereafter, as a representative example, We also consider two typical pattern speeds: The bar potential is quite easy to reproduce using Fourier coefficients since it varies smoothly on orbits. Thus, for a study in the Galactic plane, 41 complex Fourier coefficients for each triplet of actions are sufficient to approximate the value of the potential with an accuracy much better than 1%. Here it should be noted that the potential of the bar oscillates along the azimuth at a given radius and that the relative accuracy can be ill-defined when the potential passes through zero. Therefore, we define here the relative accuracy with respect to the amplitude (i.e. the maximum value) of the bar potential at a given radius. The Fourier coefficients themselves vary smoothly, as illustrated in Fig. 1, which shows the variations of a few Fourier coefficients as J R and J ϕ increase separately, justifying the use of cubic-spline interpolation to get the value of the potential at a specific position. Figure 2 demonstrates the accuracy of our method in reproducing the bar potential in the solar neighbourhood for different values of the local velocities. The potential is estimated at the same configuration space location for the whole range of rele-vant velocities, with a typical accuracy at the per cent level both in the plane and at z = 0.3 kpc. The accuracy remains very good above the plane, although with a slight bias towards lower amplitudes than the true value. More complex Fourier coefficients are needed outside of the plane. This tool is of course not limited to any specific form of the perturbing potential, the only adjustable parameter being the number of Fourier coefficients necessary to recover a given perturbing potential with a per cent-level accuracy.

Spiral potential
The potential we use for the spiral arms is the following (Cox & Gómez 2002;Monari et al. 2016) where m = 2, Ω sp is the pattern speed of the spiral arms, and where Here R sp = 1 kpc is the length parameter of the logarithmic spiral potential, h sp = 0.1 kpc the height parameter, p = −9.9 o the pitch angle, ϕ sp = −26 o the phase, and A = 683.4 km 2 s −2 the amplitude. Hereafter, we adopt a pattern speed Ω sp = 0.84Ω 0 , placing the main resonances away from (or at high azimuthal velocities in) the solar neighbourhood. Within the Galactic plane, the top panel of Fig. 3 shows our reconstruction of the spiral potential at the Sun's position, again with 41 complex Fourier coefficients. The accuracy is again below the per cent level as in the bar case, although in the spiral case there is a slight bias towards higher amplitudes with respect to the input spiral potential. This bias is very small, however, and does not affect our results. At z = 0.3 kpc, more complex Fourier coefficients are again needed (Fig. 3, bottom panel), and the accuracy reaches the per cent level, this time without bias.

Background equilibrium
From here on we work with a background axisymmetric DF f 0 as a sum of two quasi-isothermal DFs (Binney & McMillan 2011) for the thin and thick disc: The form of each DF is where R g , Ω, κ, and ν are all functions of J ϕ , and For the thin disc DF f thin , we choose h R = 2 kpc, z 0 = 0.3 kpc, h σ R = h σ z = 10 kpc,σ R (R 0 ) = 35 km s −1 , andσ z (R 0 ) = 15 km s −1 . For the thick disc DF f thick , we choose h R = 2 kpc, z 0 = 1 kpc, h σ R = 10 kpc, h σ z = 5 kpc,σ R (R 0 ) = 50 km s −1 , andσ z (R 0 ) = 50 km s −1 . Since we normalize the central surface densities of the thin and thick disc to 1, our densities can be multiplied by the appropriate surface density of the relevant stellar population to obtain physical units. The background axisymmetric potential is chosen to be Model I from Binney & Tremaine (2008), in which the above equilibrium DF f 0 is a good representation of the thin and thick disc components. In this model one has R 0 = 8 kpc and v 0 = 220 km s −1 .
The top panels of Fig. 4 display the (u, v)-plane in the solar neighbourhood within the z = 0 plane (and for w = −v z = 0) for this f 0 axisymmetric background, where u = −v R and v = v ϕ −v 0 , obtained by converting velocity-space into action-space through the epicyclic approximation and the Stäckel fudge from AGAMA. The velocity distributions are quite similar. However, as can be seen in the bottom panels of Fig. 4, the epicyclic approximation quickly becomes imprecise outside of the plane as it implies

Resonant zones
In the case of a perturbation with quasi-static amplitude that has reached its plateau, once the Fourier coefficients representing the perturbing potential have been computed (from the epicyclic approximation or from Eq. (10)) the expression for the perturbed DF can be simply expressed away from resonances with Eq. (7) as with n the order of the Fourier series (in this paper, m = 2 in both the bar and spiral cases), and where ω R , ω ϕ , and ω z can be approximated as epicyclic frequencies in the epicyclic case or can be determined with AGAMA.
The denominator of f jml may lead to a divergence in the DF when it approaches zero. Following our notation in Eq. (9), it can be expressed as The amount of the resonances is limited in the epicyclic case because, by construction, indices run only over the values j = {−1, 0, 1} and l = {−2, 0, 2}, but they can be much more numerous in the more accurate AGAMA case. For the bar potential of Eqs. (11) and (12), and choosing a pattern speed Ω p = 1.89Ω 0 as for our fiducial bar model, we explore in Figs. 5 and 6, the values of ω s, jml (J R , J ϕ , J z ) in action space when varying the pair of integer indices ( j, l). The actions are renormalized by the radial velocity dispersion of the thin disc, circular velocity, and vertical velocity dispersion of the thin disc at the Sun, respectively, to only display a relevant range of actions. Exploring indices in the range [−4, +4], it is clear that most combinations do not induce a resonance that is relevant to the dynamics of the solar neighbourhood. We only display in Figs. 5 and 6 the combination of indices (in addition to the corotation) for which a resonant zone appears in the plotted region of action space. It is clear that very few loworder resonances are indeed present in the range of actions that are truly relevant for the solar neighbourhood. To date our method has not been adapted to the projection of the DF on a plane in action space or local velocity space, and therefore works best in 3D. Therefore, we show in Fig. 7 some slices in velocity space at z = 0, denoting the location of the vertical resonances (i.e. resonances involving a non-zero l, hence involving the vertical frequency) either for a fixed value of the azimuthal velocity (and action) or for a fixed value of the radial velocity. Identifying such resonances in the vw-plane and uw-plane should allow new types of constraints to be put on the pattern speed of internal perturbers and the vertical shape of the potential of the Galaxy.  Values of log(ω s, jml ) in the (J R , J ϕ ) plane with fixed J z = 10 kpc km s −1 , for a few combinations of ( j, l) indices giving rise to resonant zones in action space (recalling that m = 2). The pattern speed Ω p here is 1.89Ω 0 . The two actions are renormalized by the radial velocity dispersion of the thin disc and the circular velocity at the Sun, respectively. The deep blue lines correspond to resonant zones. For instance, the (1, 0) case corresponds to the traditional outer Lindblad resonance (OLR) (for a non-zero J z ). Most other low-order combinations of indices did not give rise to any relevant resonant zone in the region of interest.  Interestingly, most of these resonances are very concentrated in w and vary quickly both in u and v as a function of w, making them elusive to find when stacking tracer stars in any 2D plane of velocity space, but in principle they stand out in thin slices of velocity space. Concretely, when considering a change of 10 km s −1 in vertical velocity from 5 to 15 km s −1 , the corresponding change in the location of the vertical resonance in v within the uv-plane is always larger than 10 km s −1 and typically larger (sometimes much larger) than 30 km s −1 .
Moreover, the signature of these vertical resonances in the uv-plane is rather thin, typically of the order of the km s −1 , hence much thinner than the displacement of the resonance with w. This means that, when investigating the uv-plane, vertical resonances should mostly be washed out as soon as the investigated slice is thick enough. Therefore, when investigating the DF in the uv-plane in the next subsection, we limit ourselves to the effect of l = 0 resonances. Figs. 8 and 9, for a lower pattern speed Ω p = 0.84Ω 0 , corresponding to the pattern speed of our fiducial spiral potential, a smaller number of vertical resonances are prominent in the solar neighbourhood.

As displayed in
While a specific treatment is needed in these resonant zones (e.g., Monari et al. 2017a), the signature of the resonances (and thus their location in velocity space) can clearly be identified with our linear perturbation method, and the linear perturbation treatment hereafter should accurately describe the deformations of velocity space outside of these resonant zones.

Comparing the perturbed DF for different action estimates
We are now in a position to compare the linear deformation of local velocity space for different action estimates, namely the epicyclic case used in previous works and the more A50, page 8 of 14  accurate AGAMA action estimates. Since our method works best for now in 3D velocity space, we limit ourselves to slices of zero vertical velocity at different heights and to l = 0 resonances. Figure 10 displays the f 0 + f 1 linearly perturbed distribution function at the position of the Sun in the Galactic plane for the bar potential of Sect. 2.4 and two different pattern speeds, and for the spiral potential of Sect. 2.5. As in Monari et al. (2017b), whenever f 1 > f 0 , we cap f 1 at the value of f 0 to roughly represent the resonant zone. The more rigorous approach, which we leave to further work in 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 resonant zone is not well modelled by the present method, the location and global shape of resonances should be well reproduced. We indeed highlight in Fig. 10 the zone occupied by trapped orbits at the corotation (Ω b = 1.16Ω 0 ) and OLR (Ω b = 1.89Ω 0 ) of the bar, as determined with the method of Monari et al. (2017a) both in the epicyclic and AGAMA cases (Al Kazwini et al., in prep.). While the quantitative enhancement of the DF will be slightly different from our linear treatment in these trapping zones, it is clear that the location of the resonant deformation is well captured by the method, as expected. The linear deformation outside of the resonant zones should be well described by our method as well. Interestingly enough, the linear deformation due to the bar is generally stronger in the AGAMA case, and that due to the spiral is weaker in the AGAMA case. This means that reproducing the effect of spiral arms on the local velocity distribution might require a higher amplitude when considering an accurate estimate of the action-angle variables rather than the epicyclic approximation. We speculate that this is related to the inaccuracy of the reconstruction of the potential in the epicyclic case, which causes different biases in the spiral and bar cases.
The case of the pattern speed of the bar being 1.89Ω 0 would correspond to a configuration where the Hercules stream at negative u and negative v corresponds to the 2:1 outer Lindblad resonance of the bar (e.g., Dehnen 2000;Minchev et al. 2007;Monari et al. 2017b;Fragkoudi et al. 2019). Although this happens in the resonant zone, it is interesting to note that this feature is less squashed in the more realistic AGAMA case. Moreover, a resonance unnoticed within the epicyclic approximation appears at high azimuthal velocities: we can identify this resonance as the outer 1:1 resonance of the bar (Dehnen 2000). In the spiral case, the resonant ridge at large azimuthal velocities can be identified as the corotation of the spiral pattern. Figure 11 displays the linear deformation due to the bar, for the case of pattern speed of 1.89Ω 0 , at different heights above the Galactic plane, both in the epicyclic and AGAMA cases. We again restrict ourselves to a zero vertical velocity slice and l = 0 resonances. As can be seen in this figure, the epicyclic approximation quickly becomes imprecise at large heights because it implies a stronger falloff of the density with height (as already noted in Fig. 4) while not changing the azimuthal velocity distribution (and the location of resonances in v) due to the hypothesis of complete decoupling of vertical motions.
In the AGAMA case the azimuthal velocity distribution is affected by a larger asymmetric drift at large heights, and the location of the outer Lindblad resonance of the bar in the uvplane is also displaced to lower azimuthal v at larger heights. This occurs because at fixed J ϕ the azimuthal and radial frequencies computed with AGAMA are lower at higher z, meaning that one needs to reach lower J ϕ (corresponding to orbits whose guiding radii are in the inner Galaxy) to reach the resonance.
This trend is most clearly visible at z = 1 kpc, where the epicyclic approximation does not accurately represent the location of the Hercules feature compared to the AGAMA case. Interestingly, comparing the displacement with height of the OLR Fig. 10. Distribution function from Fig. 4 in velocity space at the solar position within the Galactic plane, now perturbed to linear order by a bar (perturbing potential from Sect. 2.4) with pattern speeds Ω b = 1.16Ω 0 (left) and Ω b = 1.89Ω 0 (middle), or by a spiral pattern (perturbing potential from Sect. 2.5) with pattern speed Ω sp = 0.84Ω 0 (right). The black dashed contours represent the zones where k is equal to or less than 1, k being a quantity computed in Monari et al. (2017a) that designates the region where the orbits are trapped at the main resonance (the computation used here in the Stäckel case will be presented in detail in Al Kazwini et al., in prep.). Top row: epicyclic approximation. Bottom row: Stäckel fudge. in the case of a bar with pattern speed 1.89Ω 0 with that of the corotation in the case of a 1.16Ω 0 pattern speed, we noted that the corotation location in the uv-plane is more displaced than the OLR. This is because the corotation only depends on the azimuthal frequency, while the OLR depends on a combination of the azimuthal and radial frequencies. On the other hand, with the presently assumed background potential, we found that the displacement with height was rather independent of the pattern speed and therefore of the location of the resonance in local velocity space. We found a gradient in v of 8 km s −1 kpc −1 for the Fig. 12. Same as Fig. 11, 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.89Ω 0 and a spiral pattern (perturbing potential of Sect. 2.5) with pattern speed Ω sp = 0.84Ω 0 . corotation, 6 km s −1 kpc −1 for the OLR, and 4 km s −1 kpc −1 for the 1:1 resonance. This different displacement can also be seen when linearly adding the effect of the bar and spiral in Fig. 12, where the spacing between the 1:1 resonance of the bar and that of the corotation of the spiral increases with height.
Quantitatively, these displacements depend strongly on the background Galactic potential. This means that once the resonances potentially responsible for moving groups in the solar neighbourhood have been identified, studying their position in the uv-plane as a function of z can in principle be a powerful new way to constrain the 3D structure of the Galactic potential. This cannot be done within the epicyclic approximation. We note that marginalizing over vertical velocities instead of taking a zero-velocity slice would not compensate for these variations of the location of resonances with height but would only enhance the effect. In practice, we investigated the displacement of the location of the in-plane OLR with vertical velocities. For w = 50 km s −1 the displacement compared to w = 0 km s −1 in terms of the v-location of the resonance at z = 1 kpc is 8 km s −1 , always towards lower azimuthal velocities; however, the signal will always be dominated by the lowest w values due to the vertical orbital structure of the disc.

Adding the temporal evolution
In the previous sections we always consider a constant amplitude for the perturbing potential in order to determine an analytical expression for the perturbed DF. In this section we investigate the time dependence of the DF by choosing a time-varying amplitude for the perturbing potential.

Time-varying amplitude function
The expression we use for the time-dependent function g controlling the amplitude of the perturbation during its growth is where t f is the time at which the perturbation is completely formed, expressed in Gyr. We consider t f = 0.5 Gyr.
The motivation for this choice of growth function is its analytic simplicity, having a function starting from exactly zero at the origin, and smooth over the whole considered range. The first derivative, [π/(2t f )] sin(πt/t f ), assures the continuity at 0 and t f with both stages, fixed at 0 for t ≤ 0 and at 1 for t ≥ t f (the first derivative is thus equal to 0 at 0 and t f ).

Time-dependent perturbed distribution function
We now take the integral of Eq. (6), restricted to [0, t] (because the g function is equal to 0 on ]−∞, 0]) and integrate by parts. We take φ n (J , t ) = g(t ) h(t ) φ n (J), with h(t ) = e −imΩ p t , and we define η(t) ≡ e iθ s,n (t) iω s,n → dη = e iθ s,n (t) dt, allowing us to rewrite Eq. (6) as We can now integrate by parts and since g(0) = 0, To calculate the second part of the integral, since dg(t)/dt = π/(2t f ) sin(πt/t f ), we write, We look for a primitive G of sin(πt/t f )e iθ s,n (t) of the form Deriving G(t) with respect to t, and equating it to the integrand in Eq. (27) we get which leads to A = π/t f ω 2 s,n − (π/t f ) 2 and B = −iω s,n ω 2 s,n − (π/t f ) 2 . It should be noted that we do not exactly recover the static case at t = t f because not all derivatives of g(t) are strictly zero at the initial and final time, as assumed in M16. If a true plateau is reached after t f in an analytic fashion, the function would nevertheless converge towards the static case. How quickly this would happen is not trivial to compute. We can however compute an upper limit based on the formalism of Monari et al. (2017a). Considering that the most trapped orbits have their slow variables following the behaviour of a harmonic oscillator, and taking 2π over the frequency of this harmonic oscillator as a characteristic time for phase-mixing, we obtain a characteristic time of the order of 2 Gyr. Now we can study analytically how the linear response to a fiducial bar with Ω b = 1.89Ω 0 evolves with time. As before the method is not strictly valid at resonances, where a treatment like that used in Monari et al. (2017a) must be applied (see also Binney 2020a,b). It is nevertheless interesting to see in Fig. 13 how the linear deformation of the velocity plane evolves with time near resonances (in a patch co-moving with the bar, hence at a constant azimuthal angle to the bar), while the amplitude of the perturbation grows. The effect of the OLR appears as soon as the perturbation starts to grow. As it progressively grows, the two linear modes in the DF separate and lead to a velocity plane already very much resembling the stationary form of the perturbed DF after 0.25 Gyr, that is when g(t) = 0.5 and the perturbation is half-formed. In the absence of a pattern speed variation, it is therefore not necessarily obvious to disentangle the effect of a bar whose amplitude is growing from that of a fully formed bar with larger and constant amplitude.

Conclusion
Starting from the formalism exposed in M16, we proposed a more accurate way to determine the DF of the Galactic disc perturbed to linear order by a non-axisymmetric perturbation, using a more accurate action-angle coordinate system. First, we used the torus mapping from AGAMA to numerically compute the perturbing potential in action-angle coordinates as a Fourier series expansion over the angles. We showed that we could estimate typical non-axisymmetric perturbing potentials with an accuracy at the per cent level. The algorithm can be applied to any perturbing potential, including non-plane symmetric vertical perturbations, which will be particularly important when studying the vertical perturbations of the disc with similar methods (Rozier et al., in prep.).
We then computed the DF perturbed to linear order by a typical bar or spiral potential (or a linear combination of both), and computed the local stellar velocity distribution by converting velocities to actions and angles through the Stäckel fudge implemented in AGAMA. The results were compared to those obtained by using the epicyclic approximation. The linear deformation due to the bar is generally stronger in the AGAMA case, and that due to the spiral is weaker in the AGAMA case. This means that reproducing the effect of spiral arms on the local velocity distribution might require a higher amplitude when considering an accurate estimate of the action-angle variables rather than the epicyclic approximation. Most importantly, the epicyclic approximation is inadequate at large heights and does not change the azimuthal velocity location of the resonances due to the hypothesis of complete decoupling of vertical motions. In the AGAMA case instead, the locations of resonances are displaced to lower azimuthal v at larger heights. With the background potential used in this paper, we found a displacement in v of 8 km s −1 kpc −1 for the corotation, 6 km s −1 kpc −1 for the OLR and 4 km s −1 kpc −1 for the 1:1 resonance. Thus, the position of moving groups in the uv-plane as a function of z can be a powerful way to constrain the 3D structure of the Galactic potential.
The key to exploring this will be the DR3 of Gaia (Brown 2019) A50, page 12 of 14 with its ∼35 million radial velocities allowing us to better probe the z-axis above and below the Milky Way plane.
Finally, the temporal treatment is also an improvement over M16. We applied it to the case of a bar of growing amplitude, with an analytic evolution of the amplitude. As the bar progressively grows, the two linear modes in the DF separate, and lead to a velocity plane already very much resembling the stationary form of the perturbed DF once the perturbation is half-formed. In the absence of a pattern speed variation, it is therefore not necessarily obvious to disentangle the effect of a bar whose amplitude is growing from that of a fully formed bar with larger and constant amplitude. We explored here a peculiar form of the growth function motivated by its analytic simplicity. If the perturbation grows by linear instability, exponential growth will be more realistic. Numerical experiments are usually well fitted by a logistic function (exponential growth at the beginning and saturation to the plateau). One problem for our treatment is that the logistic function is never strictly equal to 0. In addition, there is hope that similar analytical simplifications such as those for the amplitude growth 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 the resonant zones with the linear perturbations presented in this paper, the signature of the resonances (and thus their location in velocity space) can clearly be identified with this linear perturbation method. The more rigorous approach, which we leave to further works in the context of AGAMA actions (Al Kazwini et al., in prep.), is to treat the DF with a method like that of Monari et al. (2017a) in these regions, patching these results over the linear deformation computed here. Another caveat is that the torus mapping was used to express the perturbing potential in actions and angles, but for the estimate of the local stellar velocity field, we made use of the less precise Stäckel fudge method. Therefore, another promising way for improvement would be to use the new ACTIONFINDER deep-learning algorithm (Ibata et al. 2021) to make the reverse transformation. Finally, the results presented in this paper were obtained in 3D action and velocity spaces, and were mostly presented in 2D slices: it would therefore be particularly useful to improve our algorithm by including a marginalization over any axis, for instance marginalizing over vertical velocities. This is computationally more intensive but should not, a priori, pose any conceptual problem.
The tools presented in this paper will be useful for a thorough analytical dynamical modelling of the complex velocity distribution of Milky Way disc stars as measured by past and upcoming Gaia data releases. These tools will also be useful for fully self-consistent treatments of the response of the disc to external perturbations. The ultimate goal is to adjust models to the exquisite data from Gaia, which cannot be done properly with N-body simulations due to the vast parameter space to explore.
The theoretical tools and the new code presented in this paper consequently represent a useful step in this direction.