© H. Al Kazwini et al. 2022

1. 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 dynamical 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 (Antoja et al. 2018), 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 action-angle 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.

2. Perturbing potential and perturbed DF

2.1. 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_idx_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 quasi-circular 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₀(J) be the Hamiltonian of the axisymmetric and time-independent zeroth order gravitational potential Φ₀ of the Galaxy. The equations of motion connecting actions J and the canonically conjugate angles θ are then simply

$$\begin{array}{c}\hfill \dot{\mathit{\theta }}=\frac{\partial H_0}{\partial \mathit{J}}=\mathit{\omega }(\mathit{J}),\dot{\mathit{J}}=-\frac{\partial H_0}{\partial \mathit{\theta }}=0,\end{array}$$(1)

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)=θ₀ + ωt. The Jeans theorem then tells us that the phase-space distribution function (DF) of an axisymmetric potential f = f₀(J) is in equilibrium. In other words, f₀ is a solution of the collisionless Boltzmann equation:

$$\begin{array}{c}\hfill \frac{\mathrm{d}f}{\mathrm{d}t}=0.\end{array}$$(2)

2.2. Perturbed distribution functions

We now let Φ₁ be the potential of a small perturbation to the axisymmetric background potential Φ₀ of the Galaxy, with an amplitude relative to this axisymmetric background ϵ ≪ 1. Then the total potential is Φ = Φ₀ + Φ₁ and the DF becomes, to first order in ϵ, f = f₀ + f₁, which is still a solution of the collisionless Boltzmann equation. Inserting f = f₀ + f₁ in Eq. (2), and keeping only the terms of order ϵ, leads to the linearized collisionless Boltzmann equation

$$\begin{array}{c}\hfill \frac{\mathrm{d}{f}_1}{\mathrm{d}t}+[f_0,{\mathrm{\Phi }}_1]=0,\end{array}$$(3)

where [,] is the Poisson bracket. Integrating Eq. (3) over time, from −∞ to the time t, leads to

$$\begin{array}{c}\hfill f_1(\mathit{J},\mathit{\theta },t)={\displaystyle {\int }_{-\infty }^t}\mathrm{d}{t}^{\prime }\frac{\partial f_0}{\partial {\mathit{J}}^{\prime }}({\mathit{J}}^{\prime })·\frac{\partial {\mathrm{\Phi }}_1}{\partial {\mathit{\theta }}^{\prime }}({\mathit{J}}^{\prime },{\mathit{\theta }}^{\prime },t^{\prime }),\end{array}$$(4)

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 Φ₁ can be expanded in a Fourier series as

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_1(\mathit{J},\mathit{\theta },t)=\text{Re}\{{\displaystyle \sum _{\mathit{n}}}{\varphi }_{\mathit{n}}(\mathit{J},t){\mathrm{e}}^{\mathrm{i}\mathit{n}·\mathit{\theta }}\}.\end{array}$$(5)

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Ω_pt, 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

$$\begin{array}{cc}\hfill f_1(\mathit{J},\mathit{\theta },t)=\text{Re}\{& \mathrm{i}\frac{\partial f_0}{\partial \mathit{J}}(\mathit{J})·{\displaystyle \sum _{\mathit{n}}}\mathit{n}{\displaystyle {\int }_{-\infty }^t}\mathrm{d}{t}^{\prime }{\varphi }_{\mathit{n}}({\mathit{J}}^{\prime },t^{\prime }){\mathrm{e}}^{\mathrm{i}\mathit{n}·{\mathit{\theta }}^{\prime }(t^{\prime })}\}.\hfill \end{array}$$(6)

In M16, assuming ϕ_n(J′,t′) = g(t′) h(t′) ϕ_n(J), with h(t′) = e^{−imΩ_pt′}, 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₁(J, θ, t),

$$\begin{array}{c}\hfill f_1(\mathit{J},\mathit{\theta },t)=\text{Re}\{\frac{\partial f_0}{\partial \mathit{J}}(\mathit{J})·{\displaystyle \sum _{\mathit{n}}}\mathit{n}{\varphi }_{\mathit{n}}(\mathit{J})\frac{{\mathrm{e}}^{\mathrm{i}{\theta }_{\text{s},\mathit{n}}}}{{\omega }_{\text{s},\mathit{n}}}\},\end{array}$$(7)

where we defined

$$\begin{array}{c}\hfill {\theta }_{\text{s},\mathit{n}}=\mathit{n}·\mathit{\theta }-m{\mathrm{\Omega }}_{\mathrm{p}}t,\end{array}$$(8)

$$\begin{array}{c}\hfill {\omega }_{\text{s},\mathit{n}}=\mathit{n}·\mathit{\omega }-m{\mathrm{\Omega }}_{\mathrm{p}}.\end{array}$$(9)

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 Action-based Galaxy Modelling Architecture (AGAMA; Vasiliev 2019, 2018).

2.3. 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⁻¹, azimuthal actions ranging from 1200 to 2160 kpc km s⁻¹, and vertical actions ranging from 0 to 26 kpc km s⁻¹ 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):

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_1(R,\phi ,z)& =\text{Re}\{{\displaystyle \sum _{j,l}}{\varphi }_{\mathit{jml}}(J_R,J_z,J_{\phi }){\mathrm{e}}^{\mathrm{i}(j{\theta }_R+m{\theta }_{\phi }+l{\theta }_z)}\}.\hfill \end{array}$$(10)

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 two-power-law profile. The galactocentric radius of the Sun is set at R₀ = 8 kpc, and the local circular velocity is v₀ = 220 km s⁻¹.

2.4. Bar potential

The potential we choose for the bar is a simple quadrupole potential (Weinberg 1994; Dehnen 2000) with

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{1,\mathrm{b}}(R,z,\phi ,t)=\text{Re}\{{\mathrm{\Phi }}_{\mathrm{a},\mathrm{b}}(R,z){\mathrm{e}}^{\mathrm{i}m(\phi -{\phi }_{\mathrm{b}}-{\mathrm{\Omega }}_{\mathrm{b}}t)}\},\end{array}$$(11)

where m = 2, Ω_b is the pattern speed of the bar (expressed hereafter in multiples of the angular frequency at the Sun Ω₀ = v₀/R₀, where v₀ is the local circular velocity at the galactocentric radius of the Sun R₀), 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

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{a},\mathrm{b}}(R,z)=-{\alpha }_{\mathrm{b}}{\displaystyle \frac{v_0^2}{3}}{\left({\displaystyle \frac{R_0}{R_{\mathrm{b}}}}\right)}^3{\left({\displaystyle \frac{R}{r}}\right)}^2\{\begin{array}{cc}{\left({\displaystyle \frac{r}{R_{\mathrm{b}}}}\right)}^{-3}& R\ge R_{\mathrm{b}},\\ 2-{\left({\displaystyle \frac{r}{R_{\mathrm{b}}}}\right)}^3& R<R_{\mathrm{b}},\end{array}\end{array}$$(12)

where r² = R² + z² 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₀. We use hereafter, as a representative example, R_b = 0.625 R₀, φ_b = 25^o, and α_b = 0.01. We also consider two typical pattern speeds: Ω_b = 1.89Ω₀ and Ω_b = 1.16Ω₀.

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 relevant 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.

Fig. 1. Variations in a few Fourier coefficients ϕjml(J) of the bar potential from Sect. 2.4 as JR or Jφ increase separately at Jz = 0. The actions on the abscissa axis are in kpc km s−1. The curves are very smooth, which justifies our use of the cubic splines method to interpolate.

Fig. 2. Accuracy of the reconstruction of the bar potential. Left panel (top): estimate of the bar potential from Sect. 2.4 at the Sun’s position in the plane with 41 complex Fourier coefficients in Eq. (10) and the reconstruction using cubic splines. The vertical line denotes the true value. The value in the top right inset denotes the true value in physical units. The potential is always estimated at the same configuration space location (within the plane, at the Sun’s position) but for different velocities. Left panel (bottom): relative accuracy compared to the maximum value of the bar potential at the Sun’s radius denoted in the top right inset. Right panel: Same, but at z = 0.3 kpc with 231 complex Fourier coefficients. The typical accuracy is well below the per cent level, although with a slight bias towards lower amplitudes (i.e. |Φb estimated|< |Φb input|) above the plane.

2.5. Spiral potential

The potential we use for the spiral arms is the following (Cox & Gómez 2002; Monari et al. 2016)

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{1,\mathrm{sp}}(R,z,\phi ,t)=\text{Re}\left\{{\mathrm{\Phi }}_{\mathrm{a},\mathrm{sp}}(R,z){\mathrm{e}}^{\mathrm{i}m(\phi -{\phi }_{\mathrm{sp}}-{\mathrm{\Omega }}_{\mathrm{sp}}t)}\right\},\end{array}$$(13)

where m = 2, Ω_sp is the pattern speed of the spiral arms, and

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{a},\mathrm{sp}}(R,z)=-\frac{A}{R_{\mathrm{sp}}KD}{\mathrm{e}}^{\mathrm{i}m\frac{\ln (R/R_{\mathrm{sp}})}{\tan (p)}}{\left[\mathrm{sech}\left(\frac{\mathit{Kz}}{\beta }\right)\right]}^{\beta },\end{array}$$(14)

where

$$\begin{array}{cc}& K(R)=\frac{2}{Rsin(p)},\beta (R)=K(R)h_{\mathrm{sp}}[1+0.4K(R)h_{\mathrm{sp}}],\hfill \\ \hfill & D(R)=\frac{1+K(R)h_{\mathrm{sp}}+0.3{[K(R)h_{\mathrm{sp}}]}^2}{1+0.3K(R)h_{\mathrm{sp}}}.\hfill \end{array}$$(15)

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² s⁻² the amplitude. Hereafter, we adopt a pattern speed Ω_sp = 0.84Ω₀, 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.

Fig. 3. Accuracy of the reconstruction of the spiral arms potential. Left panel (top): estimate of the spiral arms potential from Sect. 2.5 at the Sun’s position in the plane with 41 complex Fourier coefficients in Eq. (10) and the reconstruction using cubic splines. The vertical line denotes the true value. The value in the top right inset denotes the true value in physical units. The potential is always estimated at the same configuration space location (within the plane, at the Sun’s position) but for different velocities. Left panel (bottom): relative accuracy compared to the maximum value of the spiral arms potential at the Sun’s radius denoted in the top right inset. Right panel: Same, but at z = 0.3 kpc with 231 complex Fourier coefficients. The typical accuracy is again well below the per cent level in the plane, while it is around the per cent level above the plane.

3. Results and comparison with the epicyclic approximation

3.1. Background equilibrium

From here on we work with a background axisymmetric DF f₀ as a sum of two quasi-isothermal DFs (Binney & McMillan 2011) for the thin and thick disc:

$$\begin{array}{c}\hfill f_0(J_R,J_z,J_{\phi })=f_{\mathrm{thin}}+0.075f_{\mathrm{thick}}.\end{array}$$(16)

The form of each DF is

$$\begin{array}{c}\hfill f(J_R,J_z,J_{\phi })=\frac{\mathrm{\Omega }\exp (-R_{\mathrm{g}}/h_R)}{2{(2\pi )}^{3/2}\kappa {\stackrel{\sim }{\sigma }}_R^2{\stackrel{\sim }{\sigma }}_z{z}_0}\exp (-\frac{J_R\kappa }{{\stackrel{\sim }{\sigma }}_R^2}-\frac{J_z\nu }{{\stackrel{\sim }{\sigma }}_z^2}),\end{array}$$(17)

where R_g, Ω, κ, and ν are all functions of J_φ, and

$$\begin{array}{cc}& {\stackrel{\sim }{\sigma }}_R(R_{\mathrm{g}})={\stackrel{\sim }{\sigma }}_R(R_0)\exp (-\frac{R_{\mathrm{g}}-R_0}{h_{{\sigma }_R}}),\hfill \\ \hfill & {\stackrel{\sim }{\sigma }}_z(R_{\mathrm{g}})={\stackrel{\sim }{\sigma }}_z(R_0)\exp (-\frac{R_{\mathrm{g}}-R_0}{h_{{\sigma }_z}}).\hfill \end{array}$$(18)

For the thin disc DF f_thin, we choose h_R = 2 kpc, z₀ = 0.3 kpc, h_σR = h_σz = 10 kpc, ${\stackrel{\sim }{\sigma }}_R(R_0)=35\text{km}{\text{s}}^{-1}$, and ${\stackrel{\sim }{\sigma }}_z(R_0)=15\text{km}{\text{s}}^{-1}$. For the thick disc DF f_thick, we choose h_R = 2 kpc, z₀ = 1 kpc, h_σR = 10 kpc, h_σz = 5 kpc, ${\stackrel{\sim }{\sigma }}_R(R_0)=50\text{km}{\text{s}}^{-1}$, and ${\stackrel{\sim }{\sigma }}_z(R_0)=50\text{km}{\text{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₀ is a good representation of the thin and thick disc components. In this model one has R₀ = 8 kpc and v₀ = 220 km s⁻¹.

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₀ axisymmetric background, where u = −v_R and v = v_φ − v₀, 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 a sharper falloff of the density compared to the better Stäckel action estimates.

Fig. 4. Local uv-plane stellar velocity distribution at axisymmetric equilibrium and for w = 0, from Eq. (16) at (R, z, φ) = (R0, 0, 0). Left panel: epicyclic approximation in the plane (top) and at z = 0.3 kpc (bottom). Right panel: Stäckel fudge with AGAMA in the plane (top) and at z = 0.3 kpc (bottom).

3.2. 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

$$\begin{array}{c}\hfill f_1(\mathit{J},\mathit{\theta },t)=\text{Re}\{{\displaystyle \sum _{j,l=-n}^n}{f}_{\mathit{jml}}{\mathrm{e}}^{\mathrm{i}[j{\theta }_R+m({\theta }_{\phi }-{\mathrm{\Omega }}_{\mathrm{p}}t)+l{\theta }_z]}\},\end{array}$$(19)

with n the order of the Fourier series (in this paper, m = 2 in both the bar and spiral cases), and

$$\begin{array}{c}\hfill f_{\mathit{jml}}={\varphi }_{\mathit{jml}}\times \frac{j\frac{\partial f_0}{\partial J_R}+m\frac{\partial f_0}{\partial J_{\phi }}+l\frac{\partial f_0}{\partial J_z}}{j{\omega }_R+m({\omega }_{\phi }-{\mathrm{\Omega }}_{\mathrm{p}})+l{\omega }_z},\end{array}$$(20)

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

$$\begin{array}{c}\hfill {\omega }_{\text{s},jml}(J_R,J_{\phi },J_z)=j{\omega }_R+m({\omega }_{\phi }-{\mathrm{\Omega }}_{\mathrm{p}})+l{\omega }_z.\end{array}$$(21)

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Ω₀ 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 low-order resonances are indeed present in the range of actions that are truly relevant for the solar neighbourhood.

Fig. 5. Values of log(ωs, jml) in the (JR, Jφ) plane with fixed Jz = 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 Jz). Most other low-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 (JR, Jz) plane with fixed Jφ = 1759 kpc km s−1 for different (j, l) resonances. The pattern speed Ωp is that of our fiducial central bar fixed at 1.89Ω0. The two actions are renormalized by the radial velocity dispersion and the vertical velocity dispersion of the thin disc at the Sun, respectively. The deep blue lines correspond to resonance zones. Most combinations of indices explored did not give rise to any relevant resonant zone in the region of interest.

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.

Fig. 7. Values of log(ωs, jml) in the uw-plane and vw-plane. Top row: values of log(ωs, jml) at z = 0 in the uw-plane with fixed Jφ = 1759 kpc km s−1, for the various vertical resonances relevant in the solar neighbourhood (the l = 0 resonances are treated in detail in Sect. 3.3). They all appear at relatively large values of w and are very concentrated in w, varying very quickly in u as a function of w. Bottom row: values of log(ωs, jml) in the vw-plane with fixed u = 0 km s−1. The pattern speed Ωp is that of our fiducial central bar fixed at 1.89Ω0.

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⁻¹ in vertical velocity from 5 to 15 km s⁻¹, the corresponding change in the location of the vertical resonance in v within the uv-plane is always larger than 10 km s⁻¹ and typically larger (sometimes much larger) than 30 km s⁻¹.

Moreover, the signature of these vertical resonances in the uv-plane is rather thin, typically of the order of the km s⁻¹, 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.

As displayed in Figs. 8 and 9, for a lower pattern speed Ω_p = 0.84Ω₀, corresponding to the pattern speed of our fiducial spiral potential, a smaller number of vertical resonances are prominent in the solar neighbourhood.

Fig. 8. Same as Fig. 5, but with some combinations of indices giving rise to resonant zones for Ωp = 0.84Ω0.

Fig. 9. Same as Fig. 6, but with some combinations of indices giving rise to resonant zones for Ωp = 0.84Ω0.

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.

3.3. 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 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₀ + f₁ 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₁ >  f₀, we cap f₁ at the value of f₀ 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Ω₀) and OLR (Ω_b = 1.89Ω₀) 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.

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.

The case of the pattern speed of the bar being 1.89Ω₀ 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Ω₀, 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.

Fig. 11. Local stellar velocity distribution perturbed to linear order at the solar galactocentric radius and azimuth at three different heights (left: z = 0 kpc, middle: z = 0.3 kpc, right: z = 1 kpc), when perturbed by a bar (perturbing potential of Sect. 2.4) with pattern speed Ωb = 1.89Ω0. Top row: epicyclic approximation. Bottom row: Stäckel fudge. The scale of the colour bar is different in the upper and lower panels for z = 1 kpc.

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 uv-plane 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 in the case of a bar with pattern speed 1.89Ω₀ with that of the corotation in the case of a 1.16Ω₀ 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⁻¹kpc⁻¹ for the corotation, 6 km s⁻¹ kpc⁻¹ for the OLR, and 4 km s⁻¹ kpc⁻¹ 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.

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.

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⁻¹ the displacement compared to w = 0 km s⁻¹ in terms of the v-location of the resonance at z = 1 kpc is 8 km s⁻¹, 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.

4. 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.

4.1. Time-varying amplitude function

The expression we use for the time-dependent function g controlling the amplitude of the perturbation during its growth is

$$\begin{array}{c}\hfill g(t)=\frac{1-cos(\pi t/t_{\mathrm{f}})}{2},\end{array}$$(22)

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).

4.2. 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Ω_pt′}, and we define

$$\begin{array}{c}\hfill \eta (t)\equiv \frac{{\mathrm{e}}^{\mathrm{i}{\theta }_{\text{s},\mathit{n}}(t)}}{\mathrm{i}{\omega }_{\text{s},\mathit{n}}}\to \mathrm{d}\eta ={\mathrm{e}}^{\mathrm{i}{\theta }_{\text{s},\mathit{n}}(t)}\mathrm{d}t,\end{array}$$(23)

allowing us to rewrite Eq. (6) as

$$\begin{array}{cc}\hfill f_1(\mathit{J},\mathit{\theta },t)=\text{Re}\{& \mathrm{i}\frac{\partial f_0}{\partial \mathit{J}}(\mathit{J})·{\displaystyle \sum _{\mathit{n}}}\mathit{n}{\varphi }_{\mathit{n}}(\mathit{J}){\displaystyle {\int }_0^t}g(t^{\prime })\frac{\mathrm{d}\eta }{\mathrm{d}{t}^{\prime }}(t^{\prime })\mathrm{d}{t}^{\prime }\}.\hfill \end{array}$$(24)

We can now integrate by parts

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^t}g(t^{\prime })\frac{\mathrm{d}\eta }{\mathrm{d}{t}^{\prime }}(t^{\prime })\mathrm{d}{t}^{\prime }={\left[g(t^{\prime })\eta (t^{\prime })\right]}_0^t-{\displaystyle {\int }_0^t}\frac{\mathrm{d}g}{\mathrm{d}{t}^{\prime }}(t^{\prime })\eta (t^{\prime })\mathrm{d}{t}^{\prime },\end{array}$$(25)

and since g(0)=0,

$$\begin{array}{c}\hfill {\left[g(t^{\prime })\eta (t^{\prime })\right]}_0^t=g(t)\eta (t).\end{array}$$(26)

To calculate the second part of the integral, since dg(t)/dt = π/(2t_f)sin(πt/t_f), we write,

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^t}\frac{\mathrm{d}g}{\mathrm{d}{t}^{\prime }}(t^{\prime })\eta (t^{\prime })\mathrm{d}{t}^{\prime }=\frac{\pi }{2t_{\mathrm{f}}}\frac{1}{\mathrm{i}{\omega }_{\text{s},\mathit{n}}}{\displaystyle {\int }_0^t}sin\left(\frac{\pi t^{\prime }}{t_{\mathrm{f}}}\right){\mathrm{e}}^{\mathrm{i}{\theta }_{\text{s},\mathit{n}}(t^{\prime })}\mathrm{d}{t}^{\prime }.\end{array}$$(27)

We look for a primitive G of sin(πt/t_f)e^{iθ_s, n(t)} of the form

$$\begin{array}{c}\hfill G(t)=[Acos\left(\frac{\pi t}{t_{\mathrm{f}}}\right)+Bsin\left(\frac{\pi t}{t_{\mathrm{f}}}\right)]{\mathrm{e}}^{\mathrm{i}{\theta }_{\text{s},\mathit{n}}(t)}.\end{array}$$(28)

Deriving G(t) with respect to t, and equating it to the integrand in Eq. (27) we get

$$\begin{array}{c}\hfill B\frac{\pi }{t_{\mathrm{f}}}+A\mathrm{i}{\omega }_{\text{s},\mathit{n}}=0\mathrm{and}B\mathrm{i}{\omega }_{\text{s},\mathit{n}}-A\frac{\pi }{t_{\mathrm{f}}}=1,\end{array}$$(29)

which leads to

$$\begin{array}{c}\hfill A=\frac{\pi /t_{\mathrm{f}}}{{\omega }_{\text{s},\mathit{n}}^2-{(\pi /t_{\mathrm{f}})}^2}\mathrm{and}B=\frac{-\mathrm{i}{\omega }_{\text{s},\mathit{n}}}{{\omega }_{\text{s},\mathit{n}}^2-{(\pi /t_{\mathrm{f}})}^2}.\end{array}$$(30)

Substituting the Eqs. (26) and (28) into Eq. (25) results in the following expression for the perturbed DF

$$\begin{array}{cc}& f_1(\mathit{J},\mathit{\theta },t)=\text{Re}\{\frac{\partial f_0}{\partial \mathit{J}}(\mathit{J})·{\displaystyle \sum _{\mathit{n}}}\mathit{n}{\varphi }_{\mathit{n}}(\mathit{J})\hfill \\ \hfill & \times [\frac{1}{2}(1-cos\left(\frac{\pi t}{t_{\mathrm{f}}}\right))\frac{{\mathrm{e}}^{\mathrm{i}{\theta }_{\text{s},\mathit{n}}}}{{\omega }_{\text{s},\mathit{n}}}-\frac{\pi }{2t_{\mathrm{f}}}\frac{1}{{\omega }_{\text{s},\mathit{n}}}\frac{1}{{\omega }_{\text{s},\mathit{n}}^2-{(\pi /t_{\mathrm{f}})}^2}\hfill \\ \hfill & \times ((\frac{\pi }{t_{\mathrm{f}}}cos\left(\frac{\pi t}{t_{\mathrm{f}}}\right)-\mathrm{i}{\omega }_{\text{s},\mathit{n}}sin\left(\frac{\pi t}{t_{\mathrm{f}}}\right)){\mathrm{e}}^{\mathrm{i}{\theta }_{\text{s},\mathit{n}}}-\frac{\pi }{t_{\mathrm{f}}}{\mathrm{e}}^{\mathrm{i}({\theta }_{\text{s},\mathit{n}}-{\omega }_{\text{s},\mathit{n}}t)})]\}.\hfill \end{array}$$(31)

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Ω₀ 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.

Fig. 13. Local stellar velocity distribution perturbed to linear order in the Galactic plane by the bar of Sect. 2.4 with Ωb = 1.89Ω0 with the Stäckel fudge, and an amplitude of the bar growing as described in Sect. 4.1. The first 11 panels display the temporal evolution of the perturbation. The last panel displays the stationary case. The amplitude of the bar goes from 0 at t = 0 to its plateau (g(t)=1) at t = 0.5 Gyr.

5. 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⁻¹ kpc⁻¹ for the corotation, 6 km s⁻¹ kpc⁻¹ for the OLR and 4 km s⁻¹ kpc⁻¹ 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) 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.

Acknowledgments

We thank the referee, Dr Paul McMillan, for a thoughtful and constructive report which has helped improving the paper and correcting a small bug in the first version. A. S., B. F., G. M., S. R., P. R. and R. I. acknowledge funding from the Agence Nationale de la Recherche (ANR project ANR-18-CE31-0006 and ANR-19-CE31-0017) and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 834148).

References

Appendix A: PERDIGAL code presentation

All the results shown in this article were obtained with a single code, written in C++ and Python, that calculates both the perturbing potential in action-angle coordinates and the perturbed DF. The code is named the PERturbed DIstribution functions for the GALactic disc (PERDIGAL) and will be made available on request, although it may eventually be embedded in a larger Galactic dynamics toolkit. Launching the code without argument gives the explanations shown in Fig. A.1.

The calculation of the perturbed DF consists of five steps: the creation of a file storing the positions of the DF (DFpos), the creation (via AGAMA) of several files containing positions from the many orbits (FCorb) required for the following step, the determination of Fourier coefficients for the perturbing potential to create the grid (FCgrid), the determination of the Fourier coefficients for each position which the DF will be calculated at (FCforDF), and finally the calculation of the perturbed DF (PertDF).The mode ALL processes the three steps (FCgrid, FCforDF, and PertDF) at one time, without saving Fourier coeffcients from FCgrid and FCforDF in files. However, ALL should not be used without certainty about the decomposition of the potential. This process strongly depends on the initial conditions fed to the code, and a check after each part of the calculation is recommended. Finally, PDFdisp displays the perturbed DF using Matplotlib.

All useful parameters are stored in two particular files and can be modified without compiling the code. They contain parameters for the calculation of Fourier coefficients, parameters for the perturbing potentials, and others depending on whether the DF is determined at first or second order, or with a time-dependent amplitude for the perturbing potential. Figure A.2 is a logigram recalling the procedure of the calculation of the perturbed DF.

Fig. A.1. Guide for the use of the code. The first three modes deal with the determination of the Fourier coefficients when expanding 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 these modes into one, but it is not recommended because of the need to check the Fourier decomposition. EPI directly calculates the perturbed distribution function under the epicyclic approximation, as explained in M16. The FUNCTIONs correspond to different operations required before using the MODEs. DFpos creates the file storing the positions at which the distribution function will be determined. FCorb creates several files containing positions from the many orbits needed to calculate the Fourier coefficients in the FCgrid step.

Fig. A.2. Logigram recalling the procedure of the calculation of the perturbed DF. The different steps 1-5 correspond respectively to the modes DFpos, FCorb, FCgrid, FCforDF, and PertDF. There are two cases where the result of the code needs to be verified: after the third step (to check that the decomposition is correct) and after the fourth step (to check that the perturbing potential is well reproduced).

All Figures

	Fig. 1. Variations in a few Fourier coefficients ϕjml(J) of the bar potential from Sect. 2.4 as JR or Jφ increase separately at Jz = 0. The actions on the abscissa axis are in kpc km s−1. The curves are very smooth, which justifies our use of the cubic splines method to interpolate.
In the text

Fig. 2. Accuracy of the reconstruction of the bar potential. Left panel (top): estimate of the bar potential from Sect. 2.4 at the Sun’s position in the plane with 41 complex Fourier coefficients in Eq. (10) and the reconstruction using cubic splines. The vertical line denotes the true value. The value in the top right inset denotes the true value in physical units. The potential is always estimated at the same configuration space location (within the plane, at the Sun’s position) but for different velocities. Left panel (bottom): relative accuracy compared to the maximum value of the bar potential at the Sun’s radius denoted in the top right inset. Right panel: Same, but at z = 0.3 kpc with 231 complex Fourier coefficients. The typical accuracy is well below the per cent level, although with a slight bias towards lower amplitudes (i.e. |Φb estimated|< |Φb input|) above the plane.

In the text

Fig. 3. Accuracy of the reconstruction of the spiral arms potential. Left panel (top): estimate of the spiral arms potential from Sect. 2.5 at the Sun’s position in the plane with 41 complex Fourier coefficients in Eq. (10) and the reconstruction using cubic splines. The vertical line denotes the true value. The value in the top right inset denotes the true value in physical units. The potential is always estimated at the same configuration space location (within the plane, at the Sun’s position) but for different velocities. Left panel (bottom): relative accuracy compared to the maximum value of the spiral arms potential at the Sun’s radius denoted in the top right inset. Right panel: Same, but at z = 0.3 kpc with 231 complex Fourier coefficients. The typical accuracy is again well below the per cent level in the plane, while it is around the per cent level above the plane.

In the text

	Fig. 4. Local uv-plane stellar velocity distribution at axisymmetric equilibrium and for w = 0, from Eq. (16) at (R, z, φ) = (R0, 0, 0). Left panel: epicyclic approximation in the plane (top) and at z = 0.3 kpc (bottom). Right panel: Stäckel fudge with AGAMA in the plane (top) and at z = 0.3 kpc (bottom).
In the text

Fig. 5. Values of log(ωs, jml) in the (JR, Jφ) plane with fixed Jz = 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 Jz). Most other low-order combinations of indices did not give rise to any relevant resonant zone in the region of interest.

In the text

Fig. 6. Values of log(ωs, jml) in the (JR, Jz) plane with fixed Jφ = 1759 kpc km s−1 for different (j, l) resonances. The pattern speed Ωp is that of our fiducial central bar fixed at 1.89Ω0. The two actions are renormalized by the radial velocity dispersion and the vertical velocity dispersion of the thin disc at the Sun, respectively. The deep blue lines correspond to resonance zones. Most combinations of indices explored did not give rise to any relevant resonant zone in the region of interest.

In the text

Fig. 7. Values of log(ωs, jml) in the uw-plane and vw-plane. Top row: values of log(ωs, jml) at z = 0 in the uw-plane with fixed Jφ = 1759 kpc km s−1, for the various vertical resonances relevant in the solar neighbourhood (the l = 0 resonances are treated in detail in Sect. 3.3). They all appear at relatively large values of w and are very concentrated in w, varying very quickly in u as a function of w. Bottom row: values of log(ωs, jml) in the vw-plane with fixed u = 0 km s−1. The pattern speed Ωp is that of our fiducial central bar fixed at 1.89Ω0.

In the text

	Fig. 8. Same as Fig. 5, but with some combinations of indices giving rise to resonant zones for Ωp = 0.84Ω0.
In the text

	Fig. 9. Same as Fig. 6, but with some combinations of indices giving rise to resonant zones for Ωp = 0.84Ω0.
In the text

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 text

Fig. 11. Local stellar velocity distribution perturbed to linear order at the solar galactocentric radius and azimuth at three different heights (left: z = 0 kpc, middle: z = 0.3 kpc, right: z = 1 kpc), when perturbed by a bar (perturbing potential of Sect. 2.4) with pattern speed Ωb = 1.89Ω0. Top row: epicyclic approximation. Bottom row: Stäckel fudge. The scale of the colour bar is different in the upper and lower panels for z = 1 kpc.

In the text

	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.
In the text

Fig. 13. Local stellar velocity distribution perturbed to linear order in the Galactic plane by the bar of Sect. 2.4 with Ωb = 1.89Ω0 with the Stäckel fudge, and an amplitude of the bar growing as described in Sect. 4.1. The first 11 panels display the temporal evolution of the perturbation. The last panel displays the stationary case. The amplitude of the bar goes from 0 at t = 0 to its plateau (g(t)=1) at t = 0.5 Gyr.

In the text

Fig. A.1. Guide for the use of the code. The first three modes deal with the determination of the Fourier coefficients when expanding 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 these modes into one, but it is not recommended because of the need to check the Fourier decomposition. EPI directly calculates the perturbed distribution function under the epicyclic approximation, as explained in M16. The FUNCTIONs correspond to different operations required before using the MODEs. DFpos creates the file storing the positions at which the distribution function will be determined. FCorb creates several files containing positions from the many orbits needed to calculate the Fourier coefficients in the FCgrid step.

In the text

Fig. A.2. Logigram recalling the procedure of the calculation of the perturbed DF. The different steps 1-5 correspond respectively to the modes DFpos, FCorb, FCgrid, FCforDF, and PertDF. There are two cases where the result of the code needs to be verified: after the third step (to check that the decomposition is correct) and after the fourth step (to check that the perturbing potential is well reproduced).

In the text