The vertical effects of disc non-axisymmetries from perturbation theory

# The vertical effects of disc non-axisymmetries from perturbation theory: the case of the Galactic bar

Giacomo Monari, Benoit Famaey, Arnaud Siebert
Observatoire astronomique de Strasbourg, Université de Strasbourg, CNRS UMR 7550, 11 rue de l’Université, 67000 Strasbourg, France
Email: giacomo.monari@astro.unistra.fr
Released 2015 Xxxxx XX
###### Abstract

Evidence for non-zero mean stellar velocities in the direction perpendicular to the Galactic plane has been accumulating from various recent large spectroscopic surveys. Previous analytical and numerical work has shown that a “breathing mode” of the Galactic disc, similar to what is observed in the Solar vicinity, can be the natural consequence of a non-axisymmetric internal perturbation of the disc. Here we provide a general analytical framework, in the context of perturbation theory, allowing us to compute the vertical bulk motions generated by a single internal perturber (bar or spiral pattern). In the case of the Galactic bar, we show that these analytically predicted bulk motions are well in line with the outcome of a numerical simulation. The mean vertical motions induced by the Milky Way bar are small (mean velocity of less than ) and cannot be responsible alone for the observed breathing mode, but they are existing. Our analytical treatment is valid close to the plane for all the non-axisymmetric perturbations of the disc that can be described by small-amplitude Fourier modes. Further work should study how the coupling of multiple internal perturbers and external perturbers is affecting the present analytical results.

###### keywords:
Galaxy: evolution – Galaxy: kinematics and dynamics – solar neighborhood – Galaxy: structure – galaxies: spiral
pagerange: The vertical effects of disc non-axisymmetries from perturbation theory: the case of the Galactic barBpubyear: 2015

## 1 Introduction

The kinematics of stars perpendicular to the Galactic plane traditionally had (and still has) great importance, for instance as a probe of the vertical force and mass density in the Solar neighborhood (e.g., Kuijken & Gilmore 1991; Creze et al. 1998; Siebert et al. 2003; Bienaymé et al. 2014; Read 2014). In dynamical modeling, the natural zeroth order approximation is to assume a null mean vertical motion everywhere in the Galaxy (as it should be in a steady-state axisymmetric stellar system, see Binney & Tremaine 2008). However, thanks to the spatial extension and accuracy of recent surveys, it has been discovered that significant non-zero mean vertical motions do exist in the extended Solar neighborhood (Widrow et al. 2012; Williams et al. 2013; Carlin et al. 2013). These consist in patterns looking like “rarefaction-compression” waves or “breathing modes” of the disc. Whilst originally associated solely with excitations by external sources such as a passing satellite galaxy or a small dark matter substructure crossing the Galactic disc (Widrow et al. 2012, Gómez et al. 2013; Yanny & Gardner 2013; Feldmann & Spolyar 2015), it was then shown that if the density perturbation has even parity with respect to the Galactic plane (i.e., is plane-symmetric), and the vertical velocity field has odd parity (i.e., a breathing mode). Then internal non-axisymmetric perturbations, such as spiral arms, could be sufficient to cause this (Faure et al. 2014, hereafter F14). This was shown analytically by solving the linearized Euler equations in a cold fluid toy model, and confirmed through test-particle simulations for quasi-static spiral arms. Subsequently, the same effect was found in self-consistent simulations (Debattista, 2014). Such an effect from internal non-axisymmetries is of course not mutually incompatible with external perturbers playing a role in shaping the velocity field (e.g. Widrow et al., 2014, in which satellite interactions cause both bending and breathing modes depending on the velocity of the satellite), and such perturbers are themselves exciting non-axisymmetric modes such as spiral waves (Widrow & Bonner, 2015).

There is a long history of theoretical studies of the dynamical effects of disc non-axisymmetries in two dimensions in the Galactic plane, dating back to the seminal works of, e.g., Lin & Shu (1964) and Toomre (1964) on spiral arms. From the observational point of view, striking kinematic features related to the bar and spiral arms were for instance found in the Solar neighbourhood in the form of moving groups, i.e. local velocity-space substructures made of stars of very different ages and chemical compositions (Dehnen, 1998; Chereul et al., 1999; Famaey et al., 2005; Antoja et al., 2008). At a less fine-grained level, it is obvious that such local velocity-space substructures will affect the mean motions too. And indeed, it has long been well established that a non-zero mean radial motion is a natural response to disc instabilities such as spiral arms (e.g., Lin & Shu 1964, see also Binney & Tremaine 2008), independently of the exact nature of the perturber (quasi-static or transient). The same is true for the effects of the Galactic bar across the whole Galactic disc (see, e.g., Kuijken & Tremaine 1991).

Such non-zero mean radial motions have recently been detected in the extended Solar neigbourhood with the RAVE survey (Siebert et al., 2011), and shown to be consistent with the effect of spiral arms in the form of Lin-Shu density waves (Siebert et al., 2012). Velocity fluctuations have also been detected on larger scales with the APOGEE survey, and attributed to the effects of the Galactic bar (Bovy et al., 2015). Monari et al. (2014, hereafter M14), found that the bar can also explain the RAVE radial velocity gradient of Siebert et al. (2011), but did not find hints of the observed vertical bulk motions (Williams et al., 2013). On the contrary, the spiral arm model of F14 has been shown to generate a “breathing mode” qualitatively similar to what is observed for vertical motions. Generally speaking, theoretical studies of the dynamical response to disc non-axisymmetries in 3D are much less developed than in 2D (see however Patsis & Grosbol 1996; Cox & Gómez 2002).

In this paper, our aim is to provide a general analytical framework linking the vertical bulk motions generated by disc non-axisymmetries to the horizontal bulk motions. This analytic treatment goes beyond the analytic toy-model of F14 for a pressureless three-dimensional fluid, which was only suitable for an extremely cold stellar population, and which required a few unrealistic assumptions (for instance that the vertical force of the axisymmetric potential was negligible w.r.t. to the maximum vertical force of the perturber). Our reasoning hereafter is based on a linearization of the zeroth order moment of the collisionless Boltzmann equation, i.e., the continuity equation for stellar systems. As an improvement over the analytic treatment in F14, the present calculation is valid for the whole range of velocity dispersions compatible with the epicyclic approximation, and it is valid close to the plane for all non-axisymmetric perturbations of the potential described by small-amplitude Fourier modes. In particular, we show that even the Galactic bar is expected to induce vertical bulk motions in the Galactic disk. We use this particular case to compare our analytic results with a numerical study.

The paper is organized as follows. In Section 2 we linearize the continuity equation, deriving a theoretical prediction for the vertical bulk motion. More specifically, we relate it to the horizontal motions in the case where the potential perturbation is a Fourier mode of small amplitude. In Section 3 we present the outcome of a numerical test-particle simulation with a realistic model of the Galactic disc under the influence of the gravity of the Milky Way and the bar, and compare it with the analytical results. In Section 4 we conclude.

## 2 Vertical effects of perturbations

### 2.1 The linearized continuity equation

Using the cylindrical coordinates and velocities , and integrating the collisionless Boltzmann equation over velocity space, one gets the following form of the continuity equation for stellar systems:

 ∂ρ∂t+1R∂(RρuR)∂R+1R∂(ρuϕ)∂ϕ+∂(ρuz)∂z=0, (1)

where is the density of the system in configuration space, and are the average velocities, all estimated at and at the time .

For an axisymmetric stationary system, symmetric with respect to the plane , all observables depend on only, and one has

 ρ=ρ0,uR=0,uϕ=uϕ,0,uz=0. (2)

Let us now consider the case where the system is perturbed by a small non-axisymmetric perturbation in the potential, i.e.,

 Φ(R,ϕ,z,t)=Φ0(R,z)+ϵΦ1(R,ϕ,z,t), (3)

where is the unperturbed axisymmetric potential, , and has the same order of magnitude as .

Considering a small response to this small perturbation, we can write:

 uR=ϵuR,1,uϕ=uϕ,0+ϵuϕ,1,uz=ϵuz,1,ρ=ρ0+ϵρ1. (4)

Plugging Eq. (4) into Eq. (1) and dropping the terms that are , we obtain the well-known linearized continuity equation:

 ∂ρ1∂t+1R∂(Rρ0uR,1)∂R+ρ0R∂uϕ,1∂ϕ+uϕ,0R∂ρ1∂ϕ=−∂(ρ0uz,1)∂z. (5)

From this equation, it is immediately apparent that, for a given background axisymmetric model, one will be able to relate the vertical bulk motion to the horizontal responses and and the density wake of the perturbation .

### 2.2 Solution for an exponential disc

Disc non-axisymmetries respect the plane symmetry, hence on the galactic plane . This would not be true in the case where a bending mode, due to satellite interactions for example, is present (Xu et al., 2015; Widrow & Bonner, 2015). This will be the subject of further work.

From , the solution of Eq. (5) reads

 ρ0(R,z)uz,1(R,ϕ,z,t)=−∫z0F(R,ϕ,ξ,t)dξ, (6)

where

 F(R,ϕ,z,t)≡∂ρ1∂t+1R∂(Rρ0uR,1)∂R+ρ0R∂uϕ,1∂ϕ+uϕ,0R∂ρ1∂ϕ. (7)

We now assume that the unperturbed stellar system that we consider is an axisymmetric exponential disc, i.e.,

 ρ0(R,z)=ρ0(0,0)exp(−RhR−|z|hz) (8)

where and are the scale length and height respectively. We also assume that and have the same exponential vertical dependence, i.e., is constant with . With these assumptions Eq. (6) becomes

 uz,1(R,ϕ,z,t)=−∫z0e−|ξ|/hzG(R,ϕ,ξ,t)dξe−|z|/hz, (9)

where

 G(R,ϕ,z,t)≡∂~ρ∂t+uϕ,0R∂~ρ∂ϕ+uR,1R−uR,1hR+∂uR,1∂R+1R∂uϕ,1∂ϕ. (10)

### 2.3 Taylor expansion close to the plane

We now indicate only the dependence on of the functions to simplify the notation, and Taylor expand in powers of up to second order close to the plane, i.e., from Eq. (9),

 uz,1(z)≈−∫z0e−|ξ|/hz[G(0)+12∂2G∂z2(0)ξ2]dξe−|z|/hz. (11)

Therefore, from Eq. (4) and Eq. (11),

 uz(z) =ϵ sgn(z)hz[G(0)(1−e|z|/hz) +∂2G∂z2(0)(z22+hz|z|+h2z−e|z|/hzh2z)], (12)

and, averaging over , weighted by , we obtain

 Δ⟨uz⟩(z)≡2∫z0ρ(ξ)uz(ξ) dξ∫z0ρ(ξ) dξ= =2ϵ{[G(0)(hz−z)+∂2G∂z2(0)h2z(3hz−z)] −z[G(0)+∂2G∂z2(0)hz2(6hz+z)]/(ez/hz−1)}, (13)

the difference between the mean vertical velocities of stars within , above and below the Galactic plane, since . A value () implies that the stars of both Galactic hemispheres tend to move away from (towards) the Galactic plane, while for we have , like in an axisymmetric Galaxy.

Note that, whilst being continuous, the function given by Eq. (2.3) is not continuously differentiable () at . This is a consequence of the fact that the vertical density distribution is itself not at . This technical issue disappears for a density distribution, which is at . In Appendix A we compute the formulae corresponding to Eqs. (2.3)-(2.3) in the case, and show that the results are quantitatively similar for the same value of close to the plane.

### 2.4 Computing the horizontal bulk motions

At this point, we will need to estimate and . To do so, we consider the potential perturbation , the density response , and the horizontal mean velocities as Fourier -modes propagating in the disc, i.e.,

 Φ1(R,ϕ,z,t)=Re{Φa(R,z)exp[im(ϕ−Ωbt)]}. (14)
 ~ρ(R,ϕ,t)=Re{~ρa(R)exp[im(ϕ−Ωbt)]}, (15)
 uR,1(R,ϕ,z,t)=Re{uaR(R,z)exp[im(ϕ−Ωbt)]}, (16)
 uϕ,1(R,ϕ,z,t)=Re{uaϕ(R,z)exp[im(ϕ−Ωbt)]}, (17)

where we have reintroduced the , , and explicit dependencies. Dividing Eq. (5) by and integrating over we obtain the 2D continuity equation, i.e.,

 ∂~ρ∂t+Uϕ,0R∂~ρ∂ϕ+UR,1R−UR,1hR+∂UR,1∂R+1R∂Uϕ,1∂ϕ=0, (18)

where

 Ui≡∫∞−∞ρ(z)ui(z) dξ∫∞−∞ρ(z) dz, (19)

and is one of the , , , , and , functions defined above. Using the condition Eq. (18) we can derive as

 ~ρa(R)=mUaϕ/R+i(UaR/hR−UaR/R−∂UaR/∂R)m(Ωb−Uϕ,0/R). (20)

The perturbative regime cannot give accurate estimates of , , and far from the plane. However, typically these functions must decrease quite fast as a function of and tend to as goes to infinity. Therefore, we can approximate

 Ui(R,ϕ,t)≈∫ζ−ζρ(R,ξ)ui(R,ϕ,ξ,t) dξ∫ζ−ζρ(R,ξ) dξ, (21)

where is the height from the plane up to which we are computing the average and at which the integral in the perturbative regime is converging to a constant value111In the rest of the paper we adopt , as it is sufficient for the convergence of the functions in the barred case of Section 3..

We then obtain the value of and close the plane by linearizing Jeans equations, in the same way than Kuijken & Tremaine (1991) in the 2D case. Under the assumptions that the velocity dispersions and are related by the epicyclic approximation (Binney & Tremaine 2008), that the radial scale length of is larger than , and that the mixed terms in the velocity dispersion and are negligible, one obtains the two following equations for and :

 uaR(R,z)=i~Δ[m(Ωb−~Ω)∂Φa∂R−2m~ΩRΦa], (22a) uaϕ(R,z)=−1~Δ⎡⎢ ⎢⎣2~B∂Φa∂R+m(mΩb−m~Ω)RΦa⎤⎥ ⎥⎦, (22b)

where

 ~Ω(R,z)=uϕ,0(R,z)R, (23a) ~κ(R,z)2=R∂~Ω2∂R(R,z)+4~Ω2(R,z), (23b) ~Δ(R,z)=~κ2(R,z)−m2[Ωb−~Ω(R,z)]2, (23c) ~B(R,z)=−~Ω(R,z)−12R∂~Ω∂R(R,z). (23d)

Notice how the quantities Eq. (23) would yield the well known , , , functions (Binney & Tremaine 2008) on the plane for a zero-asymmetric-drift cold fluid.

The assumption that the radial scale length of is larger than is in general not valid for spiral arms, except for very cold populations. Nevertheless, our framework for relating the vertical and horizontal motions remains very general. In the case of spirals we should use the generalized equations for and from Lin & Shu (1964), including reduction factors depending on the velocity dispersion of the background stellar population (Binney & Tremaine, 2008). Hereafter, we will concentrate on the case of the bar, where the above assumption on the radial variation of remains valid even for a realistic background stellar velocity dispersion.

### 2.5 Relating the vertical and horizontal bulk motions

To get the vertical bulk motions, we need and in order to evaluate in Eq. (2.3). To do so, we need and its derivatives. Hence, we further simplify the integrals in Eq. (21) by expanding, up to the second order, the functions in Eq. (22a) and Eq. (22b). These are even function with the respect of , so

 ui(R,z)≈ui(R,0)+12∂2ui∂z2(R,0)z2. (24)

For the specific function coming from the background axisymmetric model, in Appendix B we estimate

 ∂2uϕ,0∂z2(R,0)≈Ruϕ,0(R,0)dν2dR(R), (25)

where is the vertical epicyclic frequency.

We thus now have a fully analytical expression for in Eq. (2.3), which we can apply to the case of the Galactic bar hereafter. Indeed, according to the above calculations, non-zero vertical motions should be induced by any non-axisymmetric perturbations described by small-amplitude Fourier modes, even though M14 had not found them in a bar simulation. This means that they are probably quite small, but they should be present: here, we can explicitly estimate them analytically, and then check for their presence in a numerical simulation similar to that of M14. Note however that the present analytical perturbative approach breaks down close to the corotation and Lindblad resonances where our analytical expression for diverges. Note also that we consider here only the response to the potential perturbation, and do not relate back the density response to the potential perturbation.

## 3 Application to the Galactic bar

### 3.1 Galactic potential and disc stellar population

We take for the background axisymmetric potential the Milky Way Model I by Binney & Tremaine (2008). It consists of a spheroidal dark halo and bulge, and three disc components: thin, thick, and ISM disc. The disc densities decrease exponentially with Galactocentric radius and height from the Galactic plane. The position of the Sun in this model is . In this axisymmetric potential, we consider a background stellar population described by an exponential disc of scale length and scale height , with radial and vertical velocity dispersions varying as , , where , and . The mean tangential velocity has an asymmetric drift described by Stromgren’s formula (Binney & Tremaine 2008) adapted to our case, i.e.,

 uϕ,0(R,0)=vc−σ2R2vc(κ24Ω2−1+75RhR). (26)

### 3.2 Bar potential

We describe the perturbation due to the bar by a 3D extension of the quadrupole bar used by Weinberg (1994) and Dehnen (2000), i.e. as described in Eq. (14) with and the angle measured from the long axis of the bar, with

 Φa(R,z)=V203(R0Rb)3R2r2U(r) (27)

and

 U(r)={−(r/Rb)−3for% r≥Rb,(r/Rb)3−2for r

where is the Galactocentric radius of the Sun, the local circular velocity, the bar length, and .

The amplitude in Eq. (3) then represents the ratio between the bar and axisymmetric background radial forces at .

Here we choose (Dehnen 2000), , , (Antoja et al. 2014), and (Dwek et al. 1995).

### 3.3 Analytical results

With all these inputs, we can then immediately compute from Eq. (2.3), and in Fig. 1 we plot the predicted difference between the mean vertical velocities within 300 pc above and below the Galactic plane, i.e. .

The dashed line represents the orientation of the long axis of the bar, and the dashed circles the position of the corotation and outer Lindblad resonance (OLR), i.e., where , and .

We are interested in studying the response of the background axisymmetric thin disc to the bar potential, hence we focus on the outside of the bar region (). We also limit ourselves to in this plot, since that is the interesting region where there are still non-negligible effects. In this plot the Galaxy and the bar rotate counter-clockwise.

Fig. 1 shows that the mean vertical motions induced by the Milky Way bar are small, of the order of a few tenth of , and up to about . More importantly, it shows how these motions depend on the angle from the long axis of the bar and the distance from the center of the Galaxy. The angular dependence of is that of an mode propagating in the disk while keeping the same configuration w.r.t. the bar as the bar rotates.

As noted earlier, close to the corotation and OLR, diverges, because the linear theory presents poles there, as is evident in Eq. (20), and Eqs. (22a)-(22b). Higher order expansions (and changing variables so that the new variables only vary slowly there) would be required near these resonances, where we leave our predictions as blank. Notice how the correspondence in the figure to the resonances is not exact. This is due to the fact that, since we consider a stellar disk rather than a cold fluid, in the analytical model we use which differs from for the asymmetric drift and the dependence from , and shift the position of the resonances in the model. At the OLR there is a phase shift of in : along a given fixed azimuth , the value of inside and outside the OLR has always opposite signs.

At a given radius, the function also changes sign between the regions ahead of the bar () and behind the bar (). This is due to the fact that is an odd function in , as it is the sum of odd functions, as we can find out by looking at Eqs. (10), (2.3), (15-17), (20), (22a), and (22b). In particular, the compression () is ahead of the bar between the corotation and the OLR, while it is behind the bar outside the OLR, due to the change of sign of in Eq. (23).

We note that the pattern in is actually qualitatively following the pattern, because they are both odd functions in , as can be seen from Eq. (16) and Eq. (22a) for . This is in stark contrast with what happens in the case of spiral arms, for which F14 showed a clear phase shift between the vertical and radial bulk motions. The reason for this difference is that is a pure real function in the case of the bar, but depends on in the case of spirals222Where is the pitch angle, inducing a phase shift in Eq. (16) and Eq. (22a) for spirals.

For visual comparison with the pattern of Fig. 1, we plot the analytical and patterns for the bar model in Fig. 2. In particular, we see that between the corotation and OLR and for there is a net outward motion of stars () with vertical expansion, while for the radial motion is inward () and there is a vertical compression. The same behaviour has been seen in Fig. 12 of Widrow et al. (2014), where a satellite-induced bar has been shown to create the same qualitative pattern. Outside the OLR, the pattern is reversed, both in and .

A heuristic explanation of this phenomenon is that the net radial flows are related to the shape of the orbits induced by the bar inside and outside the resonances (see Binney & Tremaine 2008). Outside the OLR, the orbits in the bar frame are aligned with the bar and retrograde, while they are perpendicular to the bar and retrograde between the corotation and OLR: this naturally explains the shape of the radial flow. Stars are then also adapting to the vertical restoring force, stronger in the inner parts than the outer parts of the Galaxy: when the flow is inward, the force is stronger and the vertical motion corresponds to a compression, and when the flow is outward it corresponds to an expansion-rarefaction. We note that the situation is different in the case of spirals because of the restoring force from the spiral potential itself.

Our analytical predictions for the vertical bulk motions induced by the bar can now be compared to a simulation akin to M14 in the next subsection.

### 3.4 Numerical results

#### 3.4.1 Initial conditions and integration time

A simple and yet very effective way to study numerically the response of a stellar disc to an external perturbation is through test-particle simulations. These consist of the numerical integration of the equations of motion of massless particles (representing the stars in the disc), accelerated by a gravitational field (representing the potential of the galaxy) which is uninfluenced by the particles themselves.

We generate initial positions and velocities for our simulations, as discrete realizations of the Shu-Schwarzschild phase-space distribution function (Shu 1969; Bienayme & Sechaud 1997; Binney & Tremaine 2008) described in Section 2.2 of F14, with the same parameters adopted in our analytical calculations. The surface density of the initial conditions obtained in this way is approximately distributed in space as an exponential disc of scale length . The radial and vertical velocity dispersion on the disc plane vary approximately as , , where , and (Bienayme & Sechaud 1997, F14).

We then integrate forward our initial conditions for a total time , from to . For the bar is absent and the axisymmetric gravitational field is only given by Model I of Binney & Tremaine (2008). This period of time allows the initial conditions to become well mixed in the background potential333The Shu-Schwarzschild distribution function is built on approximate integrals of motion for an axisymmetric galaxy: in particular, the vertical energy is a good approximate integral only very close to the plane. Because of Jeans’ theorem (Binney & Tremaine 2008), the worse the approximation of the integrals, the faster the distribution function will evolve in time.. After the initial , we obtain a stable particle configuration, with velocity dispersions at : . The vertical restoring force determines the density profile, which, at is nicely described by an exponential or profile (see Appendix A), both with scale height close to the plane.

The bar is introduced smoothly in the simulation at with the potential defined in Eq. (3), Eq. (14) and Eq. (27), its amplitude reaching only at . For the amplitude grows with time by a factor (Dehnen 2000)

 η(t)=(316ξ5−58ξ3+1516ξ+12),ξ≡2t3 Gyr−1. (29)

The amplitude then stays constant at for another . The response becomes almost stable in the rotating frame of the bar444Some transient effects are still present in the kinematics at the end of the simulation, resulting in minor asymmetries in the maps of the mean motions. To reach complete stability one should integrate for many more dynamical times, which would be unphysical (see, e.g., Mühlbauer & Dehnen 2003)., and the end of the simulation can be used to compare with our analytical predictions.

#### 3.4.2 Comparison with the analytical model

In Fig. 3 we present the mean vertical kinematics as a function of the cartesian position in the Galaxy, at the end of the simulation. We plot the difference , where () is the average of particles found at () at the end of the simulation. We present it both for particles at (top panel) and at every (bottom panel). The averages are computed inside bins of , also applying a Gaussian smoothing on a scale .

Fig. 3 shows a very good agreement with Fig. 1, with clear trends in the vertical kinematics, depending on the angle from the long axis of the bar and on the distance from the resonances in the same way as on Fig. 1. The quantity appears always to be periodic with , fluctuating twice for , so that and . Notice how the phase of the periodic oscillations changes inside and outside the OLR: the maxima of are at inside the resonance and at outside, as predicted by our analytical model.

The agreement between the analytical and numerical model is not only qualitative, but also quantitative.

In Fig. 4 we show the predictions of the analytical model for (blue line) and (red line), where is measured from the long axis of the bar. In the same plot we show the confidence bands obtained fitting the model to for particles in the simulation with (blue band), and (red band). The concordance between analytical model and simulations shows how our assumptions in the analytical model are convincing. Notice how a small phase shift is present between the analytical predictions and the simulation. This is due to the fact that the kinematics in the simulation did not reach a complete stability w.r.t. the bar disturbance, as previously mentioned.

## 4 Discussion and conclusions

In this paper we demonstrated analytically how the mean vertical motions of stars in disc galaxies are affected by small, non-axisymmetric perturbations of the potential in the form of Fourier modes (e.g., bar or spiral arms). This was previously shown numerically for spiral arms in F14, and supported by an analytical toy-model for a pressureless fluid. The present analytical treatment goes beyond this toy-model, and is valid for the response of realistic disc populations to any non-axisymmetric perturbation described by small amplitude Fourier modes.

Our general analytical solution shows that, depending on the distance from the Lindblad resonances and azimuth in the frame of the perturber, the mean vertical motion points towards or away from the galactic plane, with a bulk velocity depending on the position projected in the galactic plane and on the distance from the galactic plane.

Although M14 did not find hints of significant vertical bulk motions in a Galactic bar simulation, our present analytical treatment indicates that non-zero vertical motions should be induced by the bar too. We thus explicitly estimate them analytically, and then check for their presence in a numerical simulation similar to that of M14. These vertical mean motions induced by the Galactic bar are indeed modest in magnitude, especially if compared with those observed recently in the Solar neighborhood (Widrow et al. 2012,Williams et al. 2013; Carlin et al. 2013), but they are well existing, as we have shown here for the very first time. And the results of our test-particle simulation are well in line with our analytical prediction away from the main resonances. However, it is clear that perturbations exerting vertical forces that change more rapidly than that of the bar with the distance from the galactic plane (e.g., the spiral arms) do create much more significant mean vertical velocities.

Our analytical treatment being valid close to the plane for all the non-axisymmetric perturbations of the disc that can be described by small-amplitude Fourier modes, it will be extremely useful for interpreting the outcome of simulations including any such perturbation, and to estimate how much of the observed breathing modes can be explained by such non-axisymmetries alone. Further work should study how the coupling of multiple internal perturbers (bar+ multiple spirals) and external perturbers (satellite interactions themselves exciting spiral waves, but also bending modes) is affecting the present analytical results.

## acknowledgements

We thank the anonymous referee for useful comments. We are grateful to James Binney for pointing out a mistake in the original version of the manuscript. This work has been supported by a postdoctoral grant from the Centre National d’Etudes Spatiales (CNES) for GM.

## References

• Antoja et al. (2008) Antoja T., Figueras F., Fernández D., Torra J., 2008, A&A, 490, 135
• Antoja et al. (2014) Antoja T. et al., 2014, A&A, 563, A60
• Bienaymé et al. (2014) Bienaymé O. et al., 2014, A&A, 571, A92
• Bienayme & Sechaud (1997) Bienayme O., Sechaud N., 1997, A&A, 323, 781
• Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition, Binney, J. & Tremaine, S., ed. Princeton University Press
• Bovy et al. (2015) Bovy J., Bird J. C., García Pérez A. E., Majewski S. R., Nidever D. L., Zasowski G., 2015, ApJ, 800, 83
• Carlin et al. (2013) Carlin J. L. et al., 2013, ApJ, 777, L5
• Chereul et al. (1999) Chereul E., Crézé M., Bienaymé O., 1999, A&AS, 135, 5
• Cox & Gómez (2002) Cox D. P., Gómez G. C., 2002, ApJS, 142, 261
• Creze et al. (1998) Creze M., Chereul E., Bienayme O., Pichon C., 1998, A&A, 329, 920
• Debattista (2014) Debattista V. P., 2014, MNRAS, 443, L1
• Dehnen (1998) Dehnen W., 1998, AJ, 115, 2384
• Dehnen (2000) Dehnen W., 2000, AJ, 119, 800
• Dwek et al. (1995) Dwek E. et al., 1995, ApJ, 445, 716
• Famaey et al. (2005) Famaey B., Jorissen A., Luri X., Mayor M., Udry S., Dejonghe H., Turon C., 2005, A&A, 430, 165
• Faure et al. (2014) Faure C., Siebert A., Famaey B., 2014, MNRAS, 440, 2564, (F14)
• Feldmann & Spolyar (2015) Feldmann R., Spolyar D., 2015, MNRAS, 446, 1000
• Gómez et al. (2013) Gómez F. A., Minchev I., O’Shea B. W., Beers T. C., Bullock J. S., Purcell C. W., 2013, MNRAS, 429, 159
• Kuijken & Gilmore (1991) Kuijken K., Gilmore G., 1991, ApJ, 367, L9
• Kuijken & Tremaine (1991) Kuijken K., Tremaine S., 1991, in Dynamics of Disc Galaxies, Sundelius B., ed., p. 71
• Lin & Shu (1964) Lin C. C., Shu F. H., 1964, ApJ, 140, 646
• Monari et al. (2014) Monari G., Helmi A., Antoja T., Steinmetz M., 2014, A&A, 569, A69, (M14)
• Mühlbauer & Dehnen (2003) Mühlbauer G., Dehnen W., 2003, A&A, 401, 975
• Patsis & Grosbol (1996) Patsis P. A., Grosbol P., 1996, A&A, 315, 371
• Read (2014) Read J. I., 2014, Journal of Physics G Nuclear Physics, 41, 063101
• Shu (1969) Shu F. H., 1969, ApJ, 158, 505
• Siebert et al. (2003) Siebert A., Bienaymé O., Soubiran C., 2003, A&A, 399, 531
• Siebert et al. (2012) Siebert A. et al., 2012, MNRAS, 425, 2335
• Siebert et al. (2011) Siebert A. et al., 2011, MNRAS, 412, 2026
• Toomre (1964) Toomre A., 1964, ApJ, 139, 1217
• Weinberg (1994) Weinberg M. D., 1994, ApJ, 420, 597
• Widrow et al. (2014) Widrow L. M., Barber J., Chequers M. H., Cheng E., 2014, MNRAS, 440, 1971
• Widrow & Bonner (2015) Widrow L. M., Bonner G., 2015, MNRAS, 450, 266
• Widrow et al. (2012) Widrow L. M., Gardner S., Yanny B., Dodelson S., Chen H.-Y., 2012, ApJ, 750, L41
• Williams et al. (2013) Williams M. E. K. et al., 2013, MNRAS, 436, 101
• Xu et al. (2015) Xu Y., Newberg H. J., Carlin J. L., Liu C., Deng L., Li J., Schönrich R., Yanny B., 2015, ApJ, 801, 105
• Yanny & Gardner (2013) Yanny B., Gardner S., 2013, ApJ, 777, 91

## Appendix A Results for a sech2 vertical density distribution

We report here the analytical results for and obtained with a different unperturbed vertical density distribution, namely

 ρ0(R,z)=ρ0(0,0)exp(−RhR)sech2(zhz). (30)

With this continuously differentiable choice of density, and following the same steps as in Section 2, Eqs. (2.3)-(2.3) become

 uz(z) =−ϵhzcosh2(zhz){G(0)tanh(zhz) +∂2G∂z2(0)[h2zπ224+h2z2Li2(−e−2z/hz) −hzln(1+e−2z/hz)z+12(tanh(zhz)−1)z2]}, (31)

and,

 Δ⟨uz⟩ =−2ϵcoth(zhz){−hzG(0)ln[cosh(zhz)] +∂2G∂z2(0)[−h3z4Li3(−e−2zhz)−h3z2Li3(−e2zhz) −9h3z16ζ(3)−h2z2Li2(−e−2zhz)z+h2z2Li2(−e2zhz)z −π2h2z24z+hz2log(e−2zhz+1)z2+z33]}, (32)

where the “polylogarithm” function is defined by the power series

 Lis(z)=∞∑k=1zkks, (33)

and the Riemann function by

 ζ(s)=∞∑k=1k−s. (34)

Using these new formulae Eqs. (A)-(A), we obtain Fig. 5, corresponding to Fig. 1, and Fig. 6, corresponding to Fig. 4. These figures show how in the bar case the relative difference between the predicted bulk motions for the case of a , and density distributions are tiny (of the order of at most). We note that the initial axisymmetric stellar population in our test particle simulation is well fit close to the plane by and exponential distributions with the same scale height . However, we note that the profile is actually a better fit, as expected from the form of the Shu-Schwarzschild distribution function used for the initial conditions.

## Appendix B Estimating uϕ,0 near the plane

Jeans’ equation for , in the case of an axisymmetric system, and neglecting mixed terms can be rewritten as (Binney & Tremaine 2008)

 (35)

Differentiating twice with the respect of , assuming an exponential disc, and that the velocity dispersion near the plane is approximately constant with , we obtain

 u2ϕ,0(R,z)+uϕ,0(R,z)∂2uϕ,0∂z2(R,z)=R2∂3Φ0∂R∂z2(R,z), (36)

which, estimated at , is

 ∂2uϕ,0∂z2(R,0)=R2uϕ,0(R,0)dν2dR(R), (37)

and

 ν2(R)≡∂2Φ0∂z2(R,0), (38)

is the square of the vertical frequency.

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters