Recent studies have shown that the passage of a massive satellite through the disk of a spiral galaxy can induce vertical wobbles in the disk and produce features such as in-plane rings and phase-space spirals. Here we analyze a high-resolution N-body simulation of a live stellar disk perturbed by the recent passage of a massive dwarf galaxy that induces such perturbations. We study the implications of the phase-space structures for the estimate of the matter density through traditional Jeans modelling. The dwarf satellite excites rapid time-variations in the potential, leading to a significant bias of the local matter surface density determined through such a method. In particular, while the Jeans modelling gives reasonable estimates in the most overdense regions of the disk, we show that it tends to overestimate the dynamical surface density in underdense regions. In these regions, the phase-space spiral is indeed more marked in the surface of section of height vs. vertical velocity. This prediction can be verified with future Gaia data releases. Our finding is highly relevant for future attempts at determining the dynamical surface density of the outer Milky Way disk as a function of radius. The outer disk of the Milky Way is indeed heavily perturbed, and Gaia DR2 data have clearly shown that such phase-space perturbations are even present locally. While our results show that traditional Jeans modelling should give reliable results in overdense regions of the disk, the important biases in underdense regions call for the development of non-equilibrium methods to estimate the dynamical matter density locally and in the outer disk.
Subject headings:Galaxy: kinematics and dynamics - stars: kinematics and dynamics
Vertically perturbed Galactic disk]Implications of a time-varying Galactic potential
for determinations of the dynamical surface density
Thanks to the second Gaia data release (Gaia DR2) (Gaia Collaboration et al., 2018), a map of the kinematics of disk stars in the Galaxy can now be provided not only for the solar vicinity but also over a larger spatial extent out to a few kpc from the Sun, allowing us to explore the phase-space of the stellar disk. This will become even more topical when the Gaia DR2 will be combined with future large spectroscopic surveys such as 4MOST (de Jong et al., 2012) and WEAVE (Dalton et al., 2012; Famaey et al., 2016). In principle, this should allow us to map the dynamical surface density as a function of height at various radii in the disk, and compare this with the baryonic density to infer the distribution of dark matter around the disk (e.g. Bovy & Rix, 2013). Locally the dynamical surface mass density estimates inferred from stars oscillating vertically above and below the plane (Kuijken & Gilmore, 1991; Garbari et al., 2012; Bovy & Tremaine, 2012; Bienaymé et al., 2014; Widmark & Monari, 2019; Hagen & Helmi, 2018; Kuijken & Gilmore, 1989a, b; Siebert et al., 2003; Holmberg & Flynn, 2004; Bovy & Tremaine, 2012) have been used to probe the local dark matter density, the dark matter vertical distribution and its possible deviations from a spherical halo profile, as well as the presence of a dark disk, thereby providing potential constraints on the nature of dark matter, and on alternatives (Bienaymé et al., 2009). In particular, if the dark matter density is very low at the solar vicinity, around the midplane, the relative contributions from visible and dark matter can be disentangled by trying to constrain the gravitational potential out to larger heights (Bovy & Tremaine, 2012; Holmberg & Flynn, 2004).
However, the various values reported for the dark matter density estimates in the solar neighborhood are affected by uncertainties that are associated with the large errors in the estimate of the gas mass density and the local stellar densities (see Silverwood et al., 2016). Furthermore, all these analyses assume that the disk potential is static and in dynamical equilibrium. In particular potential collective modes in the Galactic disk are not accounted for. For example Banik et al. (2017) recently studied models of a Galactic disk out of dynamical equilibrium to show that when a single tracer population is adopted to measure the vertical force at the solar vicinity, the uncertainty in the inferred value may be of the order of 20%.
Gaia DR2 data indicate that the Milky Way stellar disk in the solar neighborhood is far from being smooth but shows evidence of substructure in the kinematic space that manifests in the form of ridges and arches. While such structures had long been known for planar motions (e.g. Dehnen, 1998; Famaey et al., 2005; Monari et al., 2017), the finding of a phase-space spiral structure reported in the plane (Antoja et al., 2018) could potentially have important implications for determinations of the dynamical surface density. This feature appears to extend well beyond the solar neighborhood cylinder in which it was originally discovered and traced out to (Laporte et al., 2018b) using Gaia DR2 (see also Bland-Hawthorn et al., 2019, for a GALAH exploration).
Previous observational facts already indicated that the stellar disk is out of dynamical equilibrium, namely the vertical bulk motion of stars discovered in the solar neighborhood. These motions manisfest as vertical oscillations of the stellar disk as measured by SEGUE, RAVE and LAMOST (e.g., Widrow et al., 2012; Williams et al., 2013; Xu et al., 2015; Carrillo et al., 2018) and may be caused by the passage of a massive satellite such as the Sagittarius dwarf galaxy (Gómez et al., 2012, 2013; Widrow et al., 2014; de la Vega et al., 2015; D’Onghia et al., 2016), although Faure et al. (2014), Debattista (2014), and Monari et al. (2016) showed that vertical oscillations in the form of breathing modes are also obtained by the passage of the stars through spiral arms. Several independent studies showed that rings, wobbles and in-plane velocity anisotropies may be the dynamical response of the stellar disk to the gravitational disturbance of a satellite galaxy (Minchev et al., 2009; Purcell et al., 2011; D’Onghia et al., 2016) or their collective effects (Kazantzidis et al., 2009; Chequers et al., 2018). In the case of a perturbation by the Sagittarius dwarf, not only do such disturbances affect the inner part of the galaxy, but they can impact the outer disk in regions where the stellar surface density decreases and the disk is more sensitive to external perturbations. In particular, the recent pre-Gaia DR2 Sgr impact models of Laporte et al. (2018a) demonstrated that it could simultaneously account for the morphology of the outer disk and amplitude of density and streaming motion fluctuations in the solar neighbourhood.
Recently, the phase-space spiral was also attributed to the buckling of the Galactic bar (Khoperskov et al., 2019) although Laporte et al. (2018b) showed that it affects stars of all ages, which is a priori more in line with a recent perturbation. It is however possible that perturbations linked to the last pericentric passage of the Sagittarius dwarf and to the buckling of the bar could co-exist in the Milky Way disk.
The extension of the spatial scale of kinematic studies to larger regions of the Galactic disk by Gaia DR2 indicated the evidence of streaming motions in the radial, azimuthal and vertical velocity as well as small-amplitude fluctuations in the velocity dispersions in the region between 3 and 8 kpc far from the Solar location. In particular, the data suggest an increase in median vertical velocity, from the inner to the outer disk, that shows a complex dependence on the azimuthal component(Katz et al., 2018).
In this Letter we use an N-body simulation of a Milky Way-sized galaxy to investigate the implications of time-varying potentials for the phase-space structure and the surface density of the stellar disk induced by the impact of a massive satellite galaxy. We show that the assumptions of dynamical equilibrium and a static potential used to infer the disk dynamical density for different azimuthal angles, and particularly at large distances from galaxy center, based on stellar kinematics may not be robust unless the time-varying effects are taken into account. This paper starts with a brief description of the numerical experiment used in this study and its peculiar phase-space structure (Sec. 2). We then describe the classical Jeans approach that one would be allowed to follow in dynamical equilibrium (Sec. 3), and show how it could mislead an observer in deducing the incorrect dynamical surface density of the disk when the potential is actually time-varying.
2. The Numerical Simulation
We analyse a high-resolution N-body simulation of a Milky Way sized galaxy performed with the parallel TreePM code GADGET-3. A detailed description of the simulation is available in the literature (Laporte et al., 2018a).
In this study the galaxy model consists of a live dark matter halo, a rotationally supported disk of stars, and a spherical stellar bulge. The dark matter mass distribution of the galaxy is initially modeled with a Hernquist profile (Hernquist, 1990). The total mass of the halo is and is sampled with 40 million particles. The radial scale length of the halo is prior to applying adiabatic contraction according to (Blumenthal et al., 1986). Each dark particle has a mass of and a gravitational softening length of . For the stellar disk we adopt an exponential radial profile with an isothermal vertical distribution:
where is the stellar disk mass, is the exponential disk scale length and is the scale-height of the vertical distribution. The disk mass is set to , , and . The stellar disk is discretized with 5 million particles. The stellar particles have masses of and gravitational softening lengths of . The stellar bulge is described by the Hernquist model, with a total mass and scale radius .
2.2. Interaction with a massive satellite
The satellite adopted in the simulation (Model L2 of Laporte et al., 2018a) has a virial mass of and sinks through dynamical friction for 6 Gyr into the host galaxy halo. It excites a wake in the dark halo, creating torques acting on the disk, before transitioning to acting directly on the disk through tides at late times. As shown in Laporte et al. (2018b), each pericentric passage excites a new phase-space spiral in the plane, similar to the one observed in the Solar neighbourhood with Gaia. Each impact induces a prompt response of the disk after the collision, with the stars around the impact region moving inward first and then outwards. A wave of enhanced density is the outcome of this contraction and expansion. Such a wave that propagates through the disk tends to form a two-armed ring-like pattern with over-dense and under-dense regions together with strong bending vertical velocities (see e.g. Gómez et al., 2013; D’Onghia et al., 2016).
Hereafter, we are studying the problem of performing Jeans modelling on such a galaxy which has generically reacted to a satellite perturbation, in a general fashion, and have chosen a random time of t=6 Gyr for the analysis. The mass model of the reacting galaxy is close to that of the Milky Way in its gross properties, but our Jeans modelling should not be quantitatively compared to Gaia data, since the model is not meant to be a perfect realization of the present-day Milky Way. However, as already shown in Laporte et al. (2018a, b), a number of qualitative and quantitative agreements with the data exist. Here, we divide the disk into six sextants, and display in Fig.1 the breathing and bending velocities in layers of 400 pc height along the vertical component of the disk. Following Katz et al. (2018), these velocities are computed as the half-sum (mean) and half-difference of the median vertical velocities in symmetric layers with respect to the Galactic mid-plane.
With the above caveat in mind, if we locate the Sun at the sextant with the azimuth angle between (which corresponds also to an under-dense region in the disk) the bending velocities are negative (i.e. they are oriented towards the south Galactic pole) for between [0,400] pc and [400, 800] pc but become positive , in qualitative agreement with the Gaia DR2 data. Furthermore, the breathing velocity changes from positive to negative values at larger heights from the midplane, again in qualitative agreement with the current observations.
At the Solar vicinity, the unprecedented precision and the high-sampling offered by the Gaia DR2 data showed that the stars are distributed in the plane defined by the vertical coordinate and the vertical velocity, with a curled spiral-shape whose density increases towards the leading edge of the feature (Antoja et al. 2018). Such a spiral shape is also clearly visible when mapping the average azimuthal velocity as a function of position in the surface of section.
This feature is interpreted as the incomplete phase-mixing phenomenon occurring in the stellar disk and a manifestation of the disk being locally out of dynamical equilibrium. Recent studies using analytic models (Antoja et al., 2018; Binney & Schönrich, 2018) or numerical simulations (D’Onghia et al., 2016; Laporte et al., 2018a; Darling & Widrow, 2019) interpreted this feature as likely caused by the coupling of vertical and in-plane motions induced by the passage of a satellite galaxy through the Galactic disk. In fact in the presence of rings induced by the impact of a satellite, the intra-ring underdense regions coincide with a local increase of the characteristic scale height (see Fig. 1 in D’Onghia et al., 2016). This is consistent also with the fact that the feature is not traced by the OB stars, which are younger than 100 Myrs and did not complete one orbital period around the Galactic center (Cheng et al., 2019). In normal conditions for disk stars, the vertical motion of the stars is almost fully decoupled from the in-plane motion with the vertical frequency of stars being usually much greater than the radial epicyclic frequency, with the stars oscillating vertically much faster than their radial motions. However, in under-dense regions, stars feel a weaker gravitational restoring force, and have their vertical frequency reduced. Thus, we expect that stars located in the under-dense regions of the disk will have the horizontal and vertical motions coupled. Hence the spiral-shape feature in the plane should be visible preferentially in the under-dense regions and should vanish more quickly in locations with enhanced density like a density wave.
Fig. 2 displays the azimuthal velocity in the plane in one of the most underdense regions of the simulated disk at t=6 Gyr, defined by the sextan with the azimuthal angle ranging (right panel) as compared to the same analysis reported in the enhanced density region in the sextant characterized by (left panel). Indeed, as expected, the spiral-shape feature is visible in the under-dense region, but has a much less well-defined shape in the over-dense region.
3. The Jeans Analysis
We now proceed with analyzing the simulation at two different times, namely (i) Gyr, after the disk has settled from the initial conditions, but before the massive satellite has had a significant impact on it, and (ii) Gyr, an arbitrary time-step 100 Myr before the last pericentric passage of the massive satellite.
In both cases, we make a vertical Jeans analysis at kpc by constructing an annulus 200 pc in width at this radius. This annulus is then broken into equally-sized azimuthal patches that approximately correspond to the view used by Katz et al. (2018) to analyze the vertical velocity field of the Gaia DR2 data (see Fig. 1).
As is traditionally done, we then combine the vertical Jeans equation with the Poisson equation to obtain an estimate of the dynamical surface density as a function of height at kpc. This Jeans estimate can then directly be compared with the actual surface density of stars and dark matter in the simulation.
The equation we use for the Jeans estimate of the dynamical surface density at kpc is (see, e.g., Hagen & Helmi, 2018)
where the three first terms correspond to the estimate of the vertical force from the vertical Jeans equation and the last term accounts for a non-flat rotation curve at the radius of interest, assuming that it does not vary significantly within the height of interest. Here, corresponds to the exponential scale-height at the height of interest, not to be confused with the scale-height of Eq.(1).
Within each azimuthal patch, we recenter the stars by subtracting the position of the peak vertical density. The stars are then split into vertical bins of equal counts to provide an approximately equal signal-to-noise in each bin when measuring the velocity dispersion. The vertical density is fit by a distribution, which allows us to compute . At Gyr, the simulated disk is nicely plane-symmetric, which allows us to fit a single to the whole distribution and test our framework by computing the total density inferred from Eq.(3), , and the actual total density, and compare it to . The results are shown in Fig. 3, and demonstrate that the method works reasonably well when the simulated disk is close to equilibrium.
At Gyr, even after correcting for the position of the peak vertical density, the important North-South asymmetries in the vertical distribution of stars force us (as an observer confronted to such asymmetries would be) to consider a separate Jeans analysis in the Northern and Southern Galactic hemispheres.
Fig. 4 shows the cumulative stellar surface mass density (SMD) (blue), the cumulative total matter SMD (stars + DM; orange), the SMD derived from Eq. 2 using the stellar kinematics (green at and red at ) as described in Section 2.1 as a function of height above the plane.
We note that, in the most over-dense regions of the stellar disk (the upper left and bottom right panels), the matter density dynamically inferred by the Jeans analysis tends to be consistent in both the Northern and Southern hemispheres, and is consistent with the total matter density present in the simulations (albeit with a slight underestimation at large heights in the sextant ). On the other hand, in the most under-dense regions of the stellar disk (the left middle and bottom panels) the matter surface density inferred from the Jeans analysis is also quite consistent in both hemispheres, but systematically overestimates the total surface density. This systematic overestimation is in line with the previous analysis of a simpler model of a disk affected by breathing modes (see Figure 6 in Banik et al. 2017). Finally, in the regions of average density, the Northern and Southern estimates tend to disagree with each other, but their average is close to the true value.
In this work, we extended the classical Jeans analysis to a simulated stellar disk impacted by a massive satellite. This collision induces a vertical wobble in the disk and in-plane rings, features that are not accounted for in the estimate of the local matter density inferred from the observed kinematics of the stars.
Our main finding is that, while the Jeans modelling gives reasonable estimates of the total surface density in the most over-dense regions of the disk, it tends to systematically overestimate the surface density, by a factor as large as 1.5, in the most under-dense regions. This validates our interpretation that the vertical motion of the stars is coupled to their in-plane motion in these under-dense regions, thereby rendering the equilibrium analysis obsolete. As a side note, we also point out that, in the regions of average density, the Jeans analysis tends to give very discrepant results in the Northern and Southern Galactic hemispheres, whose average value is close to the true one.
Our findings are highly relevant for future attempts at determining the dynamical surface density of the outer Milky Way disk as a function of radius, by combining Gaia data with future spectroscopic surveys. The important biases of Jeans modelling, especially in under-dense regions call for the development of non-equilibrium methods to estimate the dynamical matter density in the outer disk. For instance, the method recently proposed by Binney & Schönrich (2018) could make it possible to take advantage of the presence of a phase-space spiral to constrain the potential once the timing of the main perturber (for instance the last pericentric passage of the Sgr dwarf) is known. However, we note that such methods rely on the impulse approximation, whereas yet again indicated by our results the self-gravity of the disk is relevant, meaning that the method should be improved to take it into account.
Acknowledgements and References
We are grateful to Justin Read and George Lake for their valuable advice. E.D.O acknowledges support from the Vilas Associate Research Fellowship and thanks the Institute for Theory and Computation (ITC) at Harvard and the Center for Computational Astrophysics at the Flatiron Institute for the hospitality during the completion of this work. CL is supported by a CITA National Fellow award. This work is supported by ATP NASA Grant No NNX144AP53G. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number OCI-1053575.
- Antoja et al. (2018) Antoja, T. et al. 2018, Nature, 561, 360
- Banik et al. (2017) Banik, N., Widrow, L. M., and Dodelson, S. 2017, MNRAS, 464, 3775
- Bienaymé et al. (2014) Bienaymé, O. et al. 2014, Astronomy & Astrophysics, 571, A92
- Bienaymé et al. (2009) Bienaymé, O. et al. 2009, A&A, 500, 801
- Binney & Schönrich (2018) Binney, J. and Schönrich, R. 2018, MNRAS, 481, 1501
- Bland-Hawthorn et al. (2019) Bland-Hawthorn, J. et al. 2019, MNRAS
- Blumenthal et al. (1986) Blumenthal, G. R., Faber, S. M., Flores, R., and Primack, J. R. 1986, ApJ, 301, 27
- Bovy & Rix (2013) Bovy, J. and Rix, H.-W. 2013, The Astrophysical Journal, 779, 115
- Bovy & Tremaine (2012) Bovy, J. and Tremaine, S. 2012, The Astrophysical Journal, 756, 89
- Carrillo et al. (2018) Carrillo, I. et al. 2018, MNRAS, 475, 2679
- Cheng et al. (2019) Cheng, X., Liu, C., Mao, S., and Cui, W. 2019, ApJ, 872, L1
- Chequers et al. (2018) Chequers, M. H., Widrow, L. M., and Darling, K. 2018, MNRAS, 480, 4244
- Dalton et al. (2012) Dalton, G. et al. 2012, in Proc. SPIE, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, 84460P
- Darling & Widrow (2019) Darling, K. and Widrow, L. M. 2019, MNRAS, 484, 1050
- de Jong et al. (2012) de Jong, R. S. et al. 2012, in Proc. SPIE, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, 84460T
- de la Vega et al. (2015) de la Vega, A. et al. 2015, MNRAS, 454, 933
- Debattista (2014) Debattista, V. P. 2014, MNRAS, 443, L1
- Dehnen (1998) Dehnen, W. 1998, AJ, 115, 2384
- D’Onghia et al. (2016) D’Onghia, E. et al. 2016, ApJ, 823, 4
- Famaey et al. (2016) Famaey, B. et al. 2016, in SF2A-2016: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. C. Reylé, J. Richard, L. Cambrésy, M. Deleuil, E. Pécontal, L. Tresse, & I. Vauglin, 281–286
- Famaey et al. (2005) Famaey, B. et al. 2005, A&A, 430, 165
- Faure et al. (2014) Faure, C., Siebert, A., and Famaey, B. 2014, Monthly Notices of the Royal Astronomical Society, 440, 2564
- Gaia Collaboration et al. (2018) Gaia Collaboration et al. 2018, A&A, 616, A1
- Garbari et al. (2012) Garbari, S., Liu, C., Read, J. I., and Lake, G. 2012, Monthly Notices of the Royal Astronomical Society, 425, 1445
- Gómez et al. (2013) Gómez, F. A. et al. 2013, MNRAS, 429, 159
- Gómez et al. (2012) —. 2012, MNRAS, 419, 2163
- Hagen & Helmi (2018) Hagen, J. H. J. and Helmi, A. 2018, A&A, 615, A99
- Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359
- Holmberg & Flynn (2004) Holmberg, J. and Flynn, C. 2004, Monthly Notices of the Royal Astronomical Society, 352, 440
- Katz et al. (2018) Katz, D. et al. 2018, A&A, 616, A11
- Kazantzidis et al. (2009) Kazantzidis, S. et al. 2009, ApJ, 700, 1896
- Khoperskov et al. (2019) Khoperskov, S. et al. 2019, A&A, 622, L6
- Kuijken & Gilmore (1989a) Kuijken, K. and Gilmore, G. 1989a, Monthly Notices of the Royal Astronomical Society, 239, 571
- Kuijken & Gilmore (1989b) —. 1989b, Monthly Notices of the Royal Astronomical Society, 239, 605
- Kuijken & Gilmore (1991) Kuijken, K. and Gilmore, G. 1991, ApJ, 367, L9
- Laporte et al. (2018a) Laporte, C. F. P. et al. 2018a, MNRAS, 481, 286
- Laporte et al. (2018b) Laporte, C. F. P., Minchev, I., Johnston, K. V., and Gómez, F. A. 2018b, arXiv e-prints
- Minchev et al. (2009) Minchev, I. et al. 2009, MNRAS, 396, L56
- Monari et al. (2016) Monari, G. et al. 2016, MNRAS, 461, 3835
- Monari et al. (2017) Monari, G., Kawata, D., Hunt, J. A. S., and Famaey, B. 2017, MNRAS, 466, L113
- Purcell et al. (2011) Purcell, C. W. et al. 2011, Nature, 477, 301
- Siebert et al. (2003) Siebert, A., Bienaymé, O., and Soubiran, C. 2003, Astronomy & Astrophysics, 399, 531
- Silverwood et al. (2016) Silverwood, H. et al. 2016, MNRAS, 459, 4191
- Widmark & Monari (2019) Widmark, A. and Monari, G. 2019, MNRAS, 482, 262
- Widrow et al. (2014) Widrow, L. M., Barber, J., Chequers, M. H., and Cheng, E. 2014, MNRAS, 440, 1971
- Widrow et al. (2012) Widrow, L. M. et al. 2012, ApJ, 750, L41
- Williams et al. (2013) Williams, M. E. K. et al. 2013, Monthly Notices of the Royal Astronomical Society, 436, 101
- Xu et al. (2015) Xu, Y. et al. 2015, ApJ, 801, 105