© ESO 2021

1. Introduction

Determining the gravitational potential of the Milky Way is important for constraining its composition, history, and dynamical properties (Dehnen & Binney 1998; Klypin et al. 2002; Widrow et al. 2008; Weber & de Boer 2010; McMillan 2011, 2017; Kafle et al. 2014; Cole & Binney 2017; Nitschai et al. 2020; Cautun et al. 2020; Li et al. 2020). Weighing the Galaxy is especially important for understanding dark matter and its spatial distribution. The local density of dark matter influences the sensitivity of direct and indirect detection experiments (Jungman et al. 1996; Klasen et al. 2015) and could potentially depend on dark substructures, which would inform us of its particle nature (Read et al. 2008; Purcell et al. 2009; Fan et al. 2013; Ruchti et al. 2014).

The first dynamical mass measurements of the Galactic disk were carried out roughly a century ago, by Kapteyn (1922), Jeans (1922), and Oort (1932). They were able to obtain remarkably accurate estimates of the total matter density in the solar neighbourhood, under the assumption that the local population of stars was in a stable configuration. Modern methods use the same basic principles: By assuming a steady state, the stellar number density distribution and velocity distribution are interrelated via the gravitational potential (although there are also recent developments in terms of probing the gravitational potential with direct acceleration measurements, for example using pulsar timings; Chakrabarti et al. 2021). With the advent of the Gaia mission (Gaia Collaboration 2016), which recently published an early instalment of its third data release (EDR3; Gaia Collaboration 2020a), the stars of our Galaxy are observed with greater depth and precision than ever before.

One great discovery directly enabled by the Gaia survey is that of the now famous phase-space spiral (Antoja et al. 2018). This spiral is visible in the plane of vertical position and vertical velocity (where vertical refers to the direction perpendicular to the Galactic disk), either as a function of the median velocity in the azimuthal direction or as a perturbation to the stellar number density. In terms of the stellar number density in the plane of vertical position and velocity, the phase-space spiral can be seen by eye only at higher vertical energies; however, Laporte et al. (2019) demonstrated that the spiral, when plotted as a relative over- and under-density with respect to a bulk density component, is also visible at lower vertical energies.

The phase-space spiral is a clear example that the Galaxy is actually not in a steady state; other examples are ridges in stellar number counts as a function of azimuthal velocity and galactocentric radius (Antoja et al. 2018) and pre-Gaia studies of the solar neighbourhood (Minchev et al. 2009; Widrow et al. 2012) and the Milky Way disk-halo interface (Newberg et al. 2002; Rocha-Pinto et al. 2003; Xu et al. 2015; Price-Whelan et al. 2015; Bergemann et al. 2018). It remains an open question as to what extent these structures can bias dynamical mass measurements (see for example Haines et al. 2019 and Salomon et al. 2020). However, because Gaia is precise enough to resolve and identify such features, it seems reasonable that we have already entered an era where time-varying dynamical structures induce significant systematic biases that in many cases surpass the strictly statistical uncertainty of traditional methods. While most local mass measurements quote rather small statistical uncertainties, there are significant discrepancies between studies, depending for example on the choice of method, and sometimes even within studies, depending for example on the choice of stellar tracer population (see reviews by Read 2014 and de Salas & Widmark 2020 on recent determinations of the local dark matter density). Furthermore, in a study by Widmark et al. (2021), in which the total matter density is measured in subregions of the solar neighbourhood, the authors argue that their inferred matter density distributions can only be explained by biases induced by time-varying dynamics.

This work is an attempt to move beyond the traditional methods that are based on the assumption of a steady state. We demonstrate that the gravitational potential can be inferred from the shape of a phase-space spiral in the plane of vertical position and vertical velocity. Our method, as formulated in this paper, uses the relative stellar number density perturbation of the spiral with respect to a smooth background component. This background component, which we refer to as the bulk, is the quantity that is traditionally used to perform dynamical mass measurement. Our method of inference extracts information only from the shape of the spiral, while the bulk density distribution is disregarded and does not influence, nor is influenced by, the inferred gravitational potential. We tested our method on one-dimensional simulations, which were initiated in a steady state and then perturbed by applying an external force similar to that of a passing satellite. The gravitational potential was inferred with high accuracy, illustrating that time-varying structures are not necessarily obstacles to be overcome, but can be regarded as valuable assets containing useful information.

This article is structured in the following way. In Sect. 2 we describe our analytical model for the phase-space spiral, as well as our statistical model of inference. In Sect. 3 we present our simulation on which we test our model of inference. The results of those tests are found in Sect. 4. In Sects. 5 and 6 we discuss our results and conclude.

2. Analytical model of the spiral

In this work we considered a simple analytical model for the Galactic disk phase-space spiral. This analytical model was based on the three following simplifying assumptions: (i) We reduced the problem to a single spatial dimension and considered motion only in the direction perpendicular to the Galactic plane. (ii) The phase-space spiral was modelled as a first order perturbation of the stellar phase-space density distribution, inhabiting a stationary gravitational potential. As such, our analytical model neglects the self-gravity of the spiral structure. (iii) The initial perturbation is not in the shape of a spiral to begin with; it is modelled as a function that depends on, and is separable with respect to, the total vertical energy and angle of vertical oscillation (see Sect. 2.1 for the mathematical definition of these quantities).

The actual phase-space spiral of the Milky Way is of course not subject to these assumptions. Rather, it inhabits the Galaxy’s three-dimensional potential and six-dimensional phase-space structure, it is subject to self-gravity, and the initial perturbation could have a more complex structure (the Galactic disk might be subject to several or even continuous spiral-producing perturbations; see for example Hunt et al. 2018; Binney & Schönrich 2018; Laporte et al. 2019; Khoperskov et al. 2019). A discussion on the validity of these assumption is provided in Sect. 5.

2.1. Coordinate system and stellar kinematics

In order to describe the analytical model for a phase-space spiral, we begin by making the following definitions.

The vertical position, written as z and also referred to as height, is defined in relation to the disk’s mid-plane, such that the gravitational potential, written as Φ(z), has a zero-valued minimum at z = 0 pc. The vertical velocity, written as w, is defined in relation to the disk’s rest frame. Hence, a star’s vertical energy¹ is equal to

$$\begin{array}{c}\hfill E_z=\mathrm{\Phi }(z)+\frac{w^2}{2}.\end{array}$$(1)

Given a fixed value for E_z, the absolute value of w is fully determined by the absolute value of z, and vice versa.

In reality, stars are not observed from the rest frame of the disk itself, but rather from the vantage point of the Solar system. The disk frame coordinates (z, w) and solar frame coordinates (Z, W) have the relation

$$\begin{array}{cc}\hfill z& =Z+Z_{\odot },\hfill \\ \hfill w& =W+W_{\odot },\hfill \end{array}$$(2)

where Z_⊙ is the height of the Sun with respect to the disk’s mid-plane and W_⊙ is the vertical velocity of the Sun in the disk’s rest frame.

A star with a certain vertical energy will reach a maximum height, z_max, which fulfils that

$$\begin{array}{c}\hfill E_z=\mathrm{\Phi }(z_{\max }).\end{array}$$(3)

The star’s vertical oscillation has a total time period of

$$\begin{array}{c}\hfill P(E_z)={\displaystyle \oint }\frac{\mathrm{d}z}{w}=4{\displaystyle {\int }_0^{z_{\max }}\frac{\mathrm{d}z}{\sqrt{2[E_z-\mathrm{\Phi }(z)]}},}\end{array}$$(4)

where the factor 4 on the right-hand side is due to the integral only covering one of four quadrants of the oscillation in the (z, w)-plane.

We define the starting point of the vertical oscillation as passing through the disk’s mid-plane with a positive vertical velocity (z = 0, w > 0). Finally, we define an angle of oscillation (φ) as the fraction of elapsed time relative to the total period, where a complete period amounts to 2π, such that [z, w](φ) = [z, w](φ + 2πn), where n is an integer. The angle, here chosen to lie in the range [0, 2π), can be calculated according to

$$\begin{array}{c}\hfill \phi (z,w)=\{\begin{array}{cc}2\pi P^{-1}{\displaystyle {\int }_0^{|z|}}{\displaystyle \frac{\mathrm{d}{z}^{\prime }}{\sqrt{2[E_z-\mathrm{\Phi }(z^{\prime })]}}}\hfill & \mathrm{if}z\ge 0\mathrm{and}w\ge 0,\hfill \\ \pi -2\pi P^{-1}{\displaystyle {\int }_0^{|z|}}{\displaystyle \frac{\mathrm{d}{z}^{\prime }}{\sqrt{2[E_z-\mathrm{\Phi }(z^{\prime })]}}}\hfill & \mathrm{if}z\ge 0\mathrm{and}w<0,\hfill \\ \pi +2\pi P^{-1}{\displaystyle {\int }_0^{|z|}}{\displaystyle \frac{\mathrm{d}{z}^{\prime }}{\sqrt{2[E_z-\mathrm{\Phi }(z^{\prime })]}}}\hfill & \mathrm{if}z<0\mathrm{and}w<0,\hfill \\ 2\pi -2\pi P^{-1}{\displaystyle {\int }_0^{|z|}}{\displaystyle \frac{\mathrm{d}{z}^{\prime }}{\sqrt{2[E_z-\mathrm{\Phi }(z^{\prime })]}}}\hfill & \mathrm{if}z<0\mathrm{and}w\ge 0,\hfill \end{array}\end{array}$$(5)

where these four cases correspond to the four quadrants of the (z, w)-plane and it is implicit that E_z is given by Eq. (1) and P by Eq. (4)².

In our analytical model, we parametrise the gravitational potential according to

$$\begin{array}{c}\hfill \mathrm{\Phi }(z|{\rho }_h)={\displaystyle \sum _{h=1}^4\frac{4\pi G{\rho }_h}{{(2^{h-1}\times 100\mathrm{pc})}^2}log[cosh({\displaystyle \frac{z}{2^{h-1}\times 100\mathrm{pc}}})].}\end{array}$$(6)

Via the one-dimensional Poisson equation,

$$\begin{array}{c}\hfill \frac{{\partial }^2\mathrm{\Phi }(z)}{\partial z^2}=4\pi G\rho (z),\end{array}$$(7)

this corresponds to a total matter density distribution equal to

$$\begin{array}{c}\hfill \rho (z)={\displaystyle \sum _{h=1}^4{\rho }_h{cosh}^{-2}({\displaystyle \frac{z}{2^{h-1}\times 100\mathrm{pc}}}),}\end{array}$$(8)

where the four separate components labelled by the index h have scale heights of {100, 200, 400, 800} pc, respectively. This functional form allows us to emulate, with flexibility, the general characteristics of the total matter density of the Galactic disk in the solar neighbourhood (in terms of total mid-plane density, scale height, and shape; see for example McKee et al. 2015 and Schutz et al. 2018).

2.2. Spiral winding

In our analytical model, the initial perturbation was assumed to be an over-density along some initial angle (${\stackrel{\sim }{\phi }}_0$). For example, the over-density could lie along a vertical or horizontal semi-axis in the (z, w)-plane.

The angle of a star vertically oscillating in the disk evolves with time according to

$$\begin{array}{c}\hfill \stackrel{\sim }{\phi }(t,E_z)={\stackrel{\sim }{\phi }}_0+2\pi \frac{t}{P(E_z)},\end{array}$$(9)

where the period P(E_z) has an implicit dependence on Φ.

Because the initial over-density is assumed to be described by a single initial angle, which we write as ${\stackrel{\sim }{\phi }}_0$, this equation also describes the winding of the phase-space spiral and how it evolves with time. The winding of the spiral arises from variations in the period with respect to E_z. If the gravitational potential of the disk were harmonic, such that the period were invariant with respect to E_z, no winding would occur. However, the matter density of the Galactic disk decreases with distance from the mid-plane, such that the period of oscillation increases with E_z. As a result, an initial perturbation will wind into a spiral that is dragging if considered from the inside out – in other words, winding in the opposite direction of the rotation of the stars (which evolve clockwise) in the (z, w)-plane.

2.3. Spirals represented as line graphs

The shape of the phase-space spiral in the (z, w)-plane is determined by: the initial angle, ${\stackrel{\sim }{\phi }}_0$; the time since the initial perturbation, t; and the gravitational potential, Φ(z). Here, we consider a simplistic case where the spiral is a one-dimensional curve in the (z, w)-plane. This curve can be written in terms of its angle as a function of vertical energy, $\stackrel{\sim }{\phi }(E_z)$, which is a smooth and strictly decreasing function. The angle can in turn be translated into phase-space coordinates, $z[\stackrel{\sim }{\phi }(E_z)]$ and $\mathit{w}[\stackrel{\sim }{\phi }(E_z)]$, according to the equations presented in Sects. 2.1 and 2.2.

An example of such a spiral is presented in Fig. 1. For this spiral, the time since the perturbation is t = 600 Myr, the initial angle is chosen such that the spiral cuts the z-axis at 500 pc, and the gravitational potential is determined via Eq. (6) with parameters ρ_h = {0,  0.06,  0.03,  0} M_⊙pc⁻³.

Fig. 1. Phase-space spiral represented as a line graph, shown in solid blue. Two points, labelled A and B, are connected to the origin point via dashed black lines. Points A and B are both associated with a total vertical energy: Ez, A = w2/2 is given by the vertical velocity where the spiral passes through the Galactic mid-plane; and Ez, B = Φ(z) is given by gravitational potential at a point of the spiral with zero velocity. In the limit of high winding, the difference Ez, A − Ez, B would approach zero from above.

Any point on the spiral in Fig. 1 is associated with a specific vertical energy. We have highlighted points A and B, for which the vertical energy takes the form of either purely kinematic or potential energy:

$$\begin{array}{cc}\hfill E_{z,A}& =\frac{w^2}{2}=\frac{{(31.5\mathrm{km}{\mathrm{s}}^{-1})}^2}{2}=496{\mathrm{km}}^2{\mathrm{s}}^{-2},\hfill \\ \hfill E_{z,B}& =\mathrm{\Delta }\mathrm{\Phi }=\mathrm{\Phi }(500\mathrm{pc})=400{\mathrm{km}}^2{\mathrm{s}}^{-2}.\hfill \end{array}$$(10)

In the limit of high winding, the difference in energy between these two points will be small, allowing us to make measurements of the disk’s mass knowing only the (z, w) coordinates of the two points. In the specific example above, the numerical values for the two energies do have a sizeable difference. However, even in a case like this, where the spiral’s winding is not very high, the gravitational potential can be robustly and accurately inferred by considering the spiral curve over a longer segment, for example over half a period. The shape of the spiral will be strongly constrained by the gravitational potential and from the fact that $\stackrel{\sim }{\phi }(E_z)$ is a smooth function.

In order to illustrate how the shape of the spiral can change when varying the gravitational potential, we show three different examples in Fig. 2. In the upper panel of this figure, the gravitational potentials are illustrated in terms of their vertical force per mass (K_z = −∂Φ/∂z). The gravitational potential of spiral A corresponds to a matter density distribution with parameters ρ_h, A = {0,  0.05,  0.02,  0.02} M_⊙pc⁻³. The matter density distribution of spiral B is identical to that of spiral A in terms of shape, but multiplied by a factor of 1.5, giving ρ_h, B = {0,  0.075,  0.03,  0.03} M_⊙pc⁻³. Finally, the matter density distribution of spiral C has parameters ρ_h, C = {0.1,  0,  0.02,  0.02} M_⊙pc⁻³. This is identical to that of spiral A for z ≳ 400 pc (also in terms of the gravitation force) but differs at lower heights.

Fig. 2. Three phase-space spirals (labelled A, B, and C) inhabiting three different gravitational potentials. Top panel: vertical force per mass (Kz = −∂Φ/∂z) as a function of height for the respective gravitational potentials. Bottom panel: three spirals in the (z, w)-plane, where their winding time (t) and initial angle (${\stackrel{\sim }{\phi }}_0$) are chosen such that the spirals have the same phase ($\stackrel{\sim }{\phi }$) and phase derivative ($\partial \stackrel{\sim }{\phi }/\partial z_{\max }$) when their vertical energy is equal to Ez = Φ(800 pc).

The three spirals, shown in the bottom panel of Fig. 2, are individually normalised in terms of t and ${\stackrel{\sim }{\phi }}_0$ to have the same winding behaviour for E_z = Φ(800 pc). This makes, in some sense, the three spirals maximally similar to one another. Despite this, the three spirals can be differentiated from one another as a result of inhabiting different gravitational potentials. In terms of shape, spirals A and B are identical; they differ solely by a factor of $\sqrt{1.5}$ in their velocity coordinate, w. Spirals A and C are practically identical when z ≳ 400 pc and E_z ≳ Φ(400 pc). They differ slightly when z ≃ 0 and E_z ≳ Φ(400 pc) due to spiral C inhabiting a potential that is slightly steeper close to the mid-plane. In the innermost region of the plot, where E_z ≲ Φ(400 pc), spiral C also differs from spiral A in terms of its winding behaviour because the oscillation period differs significantly for low z_max and E_z.

The examples presented above illustrate the general principles of how the gravitational potential affects the phase-space spiral shape, and in effect how the latter can be used to infer the former. The total amount of winding is affected by the total mass of the disk as well as t, so in this sense these quantities are degenerate. However, as we see in the comparison between spirals A and B, for which the total mass differs by a factor 1.5 and t differs by $1/\sqrt{1.5}$, the shape of the spirals are dramatically different despite having the same amount of winding in terms of $\partial \stackrel{\sim }{\phi }/\partial z_{\max }$. By comparing spirals A and C, we also see that the spiral shape in a certain segment, which we can define by some range in E_z, is sensitive to the gravitational potential at the associated energy. For example, the gravitational potential close to the mid-plane affects the spiral’s winding behaviour, especially at low vertical energies.

2.4. Complete spiral model

In this work we fitted an analytical phase-space spiral to data coming from a one-dimensional simulation of the Galactic disk. In our complete analytical model of inference, the phase-space spiral was not modelled as the idealised one-dimensional structure discussed in Sect. 2.3, but instead as a continuous function with respect to (z, w). This full model is described here.

The free parameters of our model of inference are listed in Table 1. The free parameters, encapsulated in the quantity Ψ, are split into three subgroups that determine the bulk phase-space density (Ψ_bulk), the relative phase-space density perturbation of the spiral (Ψ_spiral), and the current phase-space position of the Sun (Ψ_⊙). These model components are described in detail below.

The spiral is modelled as a relative perturbation with respect to a bulk tracer stellar density. We modelled the bulk density in the (z, w)-plane as a sum of bivariate Gaussians labelled by the index k, which are all centred on (z, w) = (0, 0) and have diagonal covariance matrices, according to

$$\begin{array}{c}\hfill B(z,w|{\mathbf{\Psi }}_{\mathrm{bulk}})={\displaystyle \sum _k{a}_k\mathcal{N}(\left[\begin{array}{c}z\\ w\\ \end{array}\right],\left[\begin{array}{cc}{\sigma }_{z,k}^2& 0\\ 0& {\sigma }_{w,k}^2\\ \end{array}\right]),}\end{array}$$(11)

where Ψ_bulk includes the Gaussian weights a_k as well as the dispersions σ_z, k and σ_w, k. The function

$$\begin{array}{c}\hfill \mathcal{N}(\mathit{p},\mathbf{\Sigma })\equiv \frac{exp(-{\displaystyle \frac{1}{2}}{\mathit{p}}^{\top }{\mathbf{\Sigma }}_{\mathit{p}}^{-1}\mathit{p})}{2\pi \sqrt{|{\mathbf{\Sigma }}_{\mathit{p}}|}}\end{array}$$(12)

is a bivariate Gaussian.

The bulk stellar density is the quantity that is usually used in order to perform dynamical mass measurements, under the assumption that this distribution is in a steady state. In the method used in this work, the bulk is merely a background that needs to be subtracted in order to extract the spiral structure. The bulk density, as expressed in Eq. (12), is not required to fulfil the stationary Boltzmann equation or the Jeans equations; it is completely free to vary in terms of its free parameters (meaning weights and dispersions), regardless of the gravitational potential. In this manner, the gravitational potential is completely uninformed and independent of the bulk stellar density and will be inferred only from the shape of the spiral.

The phase-space spiral was modelled as a sum of an anti-symmetric and a symmetric component, corresponding to a single- or double-armed spiral, parametrised by amplitudes α and β according to

$$\begin{array}{cc}\hfill S(z,w|{\mathbf{\Psi }}_{\mathrm{spiral}})=\alpha & cos[\phi (z,w|{\rho }_h)-\stackrel{\sim }{\phi }(z,w,{\rho }_h,t,{\stackrel{\sim }{\phi }}_0)]+\hfill \\ \hfill \beta & cos\{2[\phi (z,w|{\rho }_h)-\stackrel{\sim }{\phi }(z,w,{\rho }_h,t,{\stackrel{\sim }{\phi }}_0)]\}.\hfill \end{array}$$(13)

Here, φ(z, w | ρ_h) is the intrinsic phase of the (z, w) coordinates according to Eq. (5), although here it is written with an explicit dependence on the parameters that determine the gravitational potentia, and $\stackrel{\sim }{\phi }(z,\mathit{w},\mathrm{\Phi },t,{\stackrel{\sim }{\phi }}_0)$ is the phase associated with the spiral according to Eq. (9). Because the spiral density is a relative density, the quantities α and β are unit-less and constrained to fulfil that (α, β) > 0 and α + β < 1. The parameters ρ_h determine the gravitational potential according to Eq. (6).

In the limit of low vertical energies, the self-gravity of the phase-space spiral is not negligible, and the spiral structure in this innermost region does not form (as found by Darling & Widrow 2019). For this reason, we imposed a smooth inner boundary to the relative spiral density, corresponding to E_z = Φ(400 pc), within which there is no spiral. This inner boundary is defined by

$$\begin{array}{c}\hfill m(z,w|{\rho }_h)=\mathrm{sigm}[\frac{E_z(z,w|{\rho }_h)-\mathrm{\Phi }(400\mathrm{pc}|{\rho }_h)}{\mathrm{\Phi }(400\mathrm{pc}|{\rho }_h)-\mathrm{\Phi }(380\mathrm{pc}|{\rho }_h)}],\end{array}$$(14)

where

$$\begin{array}{c}\hfill \mathrm{sigm}(x)\equiv \frac{1}{1+exp(-x)}\end{array}$$(15)

is a sigmoid function that outputs values in the range (0, 1).

The total phase-space density of our analytical model is equal to

$$\begin{array}{cc}\hfill f(z,w|{\mathbf{\Psi }}_{\mathrm{bulk}},{\mathbf{\Psi }}_{\mathrm{spiral}})& =B(z,w|{\mathbf{\Psi }}_{\mathrm{bulk}})\hfill \\ \hfill & \times [1+m(z,w|{\rho }_h)S(z,w|{\mathbf{\Psi }}_{\mathrm{spiral}})],\hfill \end{array}$$(16)

where Ψ_bulk = {a_k, σ_z, k, σ_w, k} is the subset of parameters that affect the bulk density and ${\mathbf{\Psi }}_{\mathrm{spiral}}=\{{\rho }_h,t,{\stackrel{\sim }{\phi }}_0,\alpha ,\beta \}$ are the parameters that determine the relative phase-space perturbation of the spiral.

2.5. Data histogram and mask function

In the model of inference, the bulk and spiral densities are fitted to data, where the data have the form of a two-dimensional histogram in the (Z, W)-plane of the solar frame coordinates. This histogram is the number count of observed stars in bins of size (20 pc)×(1 km s⁻¹), written as d_i, j, where (i, j) labels the respective bins.

The analytical model is fitted to a circular area of this two-dimensional histogram, defined by a mask function in the (Z, W)-plane with smooth outer boundaries, roughly corresponding to a constant E_z. The mask function is defined as

$$\begin{array}{c}\hfill M(Z,W)=\mathrm{sigm}\{10[(\frac{Z}{Z_{\lim .}}{)}^2+(\frac{W}{W_{\lim .}}{)}^2-1]\},\end{array}$$(17)

where sigm is the sigmoid function defined in Eq. (15). The mask function is visible in Fig. 3.

Fig. 3. Mask function M(Z, W) in the (Z, W)-plane, with boundary values Zlim. = 700 pc and Wlim. = 40 km s−1. See Eq. (17) for details.

The reason for applying such a mask has to do with the outer boundaries for which the analytical spiral model can be assumed to be a valid description. In the high energy limit, the number density of stars decreases and the spiral becomes less pronounced. Furthermore, there are additional issues that could cause systematic errors in the case of the actual Milky Way. The simplifying assumption of only considering the vertical dimension will be problematic, especially for stars with high vertical energies, which have more complex dynamical behaviour.

2.6. Likelihood and fitting procedure

The likelihood is given by the Poisson count comparison of the model and data of the respective bins, in the circular area defined by the mask function. The logarithm of the likelihood is equal to

$$\begin{array}{cc}\hfill ln& \mathcal{L}({\mathit{d}}_{i,j}|\mathbf{\Psi })=-{\displaystyle \sum _{i,j}M(Z_i,W_j){\displaystyle \frac{{[{\mathit{d}}_{i,j}-f(Z_i+Z_{\odot },W_j+W_{\odot },\mathbf{\Psi })]}^2}{2f(Z_i+Z_{\odot },W_j+W_{\odot },\mathbf{\Psi })}}}\hfill \\ \hfill & +\{\mathrm{constant}\mathrm{term}\}.\hfill \end{array}$$(18)

We sought to minimise this likelihood in our fit, which would allow us to neglect any constant term in the above expression.

Our model was fitted in two separate steps. In the first step we performed a joint fit of the bulk density and solar parameters to the data while disregarding the spiral; equivalently stated, we minimised Eq. (18) with respect to Ψ_bulk and Ψ_⊙, while (α, β) = 0 remained fixed. In a second step, the bulk density and solar parameters remained fixed, while we fitted the relative density of the phase-space spiral; in other words, we minimised Eq. (18) with respect to Ψ_spiral, while Ψ_bulk and Ψ_⊙ remained fixed. Only in the second step of this process, where the spiral is fitted, does the gravitational potential vary and affect the likelihood.

In order to avoid any fitting artefacts with regards to the boundary of the mask function, defined in Eq. (17), the two steps of the fitting procedure use slightly different masks. For the first step, where the bulk density is fitted, we used a larger mask, with Z_lim. = 800 pc and W_lim. = 44 km s⁻¹. For the second step, where the relative stellar density of the spiral is fitted, we used a smaller mask, with Z_lim. = 700 pc and W_lim. = 40 km s⁻¹.

2.7. Model implementation

The first step of the minimisation procedure is straightforward and computationally fast as it only requires fitting a two-dimensional Gaussian mixture model to data. However, the second step, where the spiral is fitted, is significantly more expensive. Calculating the spiral likelihood of Eq. (18) requires calculating the intrinsic angle φ(Z_i, W_j) as a numerical integral for each bin of the two-dimensional histogram, iterating over (i, j). This is done for each individual step of the minimisation algorithm.

In order to make the algorithm computationally tractable, the method was implemented in TENSORFLOW, allowing for efficient minimisation using the Adam optimiser (Kingma 2014). Minimising the spiral likelihood function is still computationally demanding and requires several hundred CPU hours.

3. Simulations

In this section we describe the simulations on which we tested our model of inference. The simulations are one-dimensional, constrained to the (z, w)-plane. The phase-space distribution of the simulations were initiated in a steady state and then perturbed by an external force, creating phase-space over- and under-densities that wind into a spiral with time. We tested our method on two simulations, labelled simulation A and simulation B.

The simulations were run with 10⁵ equally massive particles, labelled by the index i, representing stars and gas of the Galactic disk. Because the simulations are only in one dimension, the particles can be thought of as sheets of mass extending in the directions parallel to the Galactic disk. In the simulations, the gravitational acceleration is equal to

$$\begin{array}{c}\hfill K_z(z)=\frac{\mathrm{d}w}{\mathrm{d}t}=-4\pi G[{\rho }_{\mathrm{DM}}z+{\displaystyle \sum _i}m\mathrm{sign}(z-z_i)],\end{array}$$(19)

where ρ_DM is a fixed constant representing a constant dark matter density, m is the surface mass per particle (in units M_⊙pc⁻²), z_i is the height of the ith particle, and

$$\begin{array}{c}\hfill \mathrm{sign}(x)=\{\begin{array}{cc}1\hfill & \mathrm{if}x>0,\hfill \\ 0\hfill & \mathrm{if}x=0,\hfill \\ -1\hfill & \mathrm{if}x<0.\hfill \end{array}\end{array}$$(20)

3.1. Initial state

The simulation particles were drawn from three separate distributions, representing stars of the thin and thick disk, as well as a component of cold gas. The initial configuration of the particles was in a steady state, and the three respective components were assumed to be isothermal. The component of halo dark matter was not represented by dynamical particles in the simulation, but by the fixed constant ρ_DM. The three dynamical components are defined by their respective velocity dispersions (σ_w, thin, σ_w, thick, σ_w, gas) and mid-plane matter densities (ρ_0, thin, ρ_0, thick, ρ_0, gas).

Given these specifications, the matter density of the respective components, here labelled by x = {thin, thick, gas}, follows the relation

$$\begin{array}{c}\hfill \rho {(z)}_x={\rho }_{0,x}exp[-\frac{\mathrm{\Phi }(z)}{{\sigma }_{w,x}^2}],\end{array}$$(21)

where Φ(z) is found via the Poisson equation of Eq. (7) and the total matter density, which is equal to

$$\begin{array}{c}\hfill \rho (z)=\rho {(z)}_{\mathrm{thin}}+\rho {(z)}_{\mathrm{thick}}+\rho {(z)}_{\mathrm{gas}}+{\rho }_{\mathrm{DM}}.\end{array}$$(22)

By solving this system of equations, assuming boundary conditions Φ(0 pc) = 0 km² s⁻² and ∂Φ/∂z|_{z = 0 pc} = 0 km s⁻², the phase-space distributions of the respective components are fully determined.

The 10⁵ particles of our simulations were randomly drawn realisations of the stellar and gas steady state phase-space distributions, where the surface mass per particle is equal to 10⁻⁵ times the total surface mass of the gas and stellar components. When generating the data, meaning the two-dimensional histogram d_i, j representing observed stars, only the particles drawn from the thin and thick stellar disk components were used; particles belonging to the gas component were ignored because they would not be observed by the Gaia survey.

An example of the initial configuration of the matter density distributions is shown in Fig. 4. The values for ρ_DM, σ_w, x, and ρ_0, x correspond to those of simulation A, which are listed in Table 2. The initial matter density distributions are mirror-symmetric with respect to the mid-plane (z → −z).

Fig. 4. Example of the matter density distributions from which the particles of our one-dimensional simulations are drawn. In this plot, the distribution corresponds to that of simulation A; see Sect. 3 and Table 2 for details.

3.2. Perturbation

After the initial set-up, we perturbed the simulation by applying an external force for a limited amount of time. The external force that created the perturbation is of the form

$$\begin{array}{c}\hfill K_{z,\mathrm{sat}.}(t,z)=-4\pi G\mathrm{\Sigma }exp(-\frac{t^2}{2{\sigma }_t^2})tanh[\frac{z-(z_{0,\mathrm{sat}.}+w_{\mathrm{sat}.}t)}{H}],\end{array}$$(23)

where Σ is a surface mass, σ_t corresponds to the duration of the external force, H is a softening length, z_0, sat. is the position of the satellite at time t = 0, and w_sat. is its constant vertical velocity. The perturbation is active for a range in time of ( − 2.5σ_t,  2.5σ_t). The time t is defined with respect to the instance when this external force was maximal (t = 0 yr). The specific numerical values we used for the simulations in this work are provided in Sect. 3.3 and Table 2.

This perturbation can be thought of as a massive satellite that passes through the Galactic disk at some angle of incidence. Due to its velocity parallel to the disk, it only affects a disk region for a limited amount of time. The softening length would correspond to the physical scale of the satellite.

We also added a scattering component, written as s_n, that is unique for each separate particle and is randomly distributed according to a Gaussian centred on unity and with a standard deviation of 1/2. The particles are affected to different extents by the external force, according to

$$\begin{array}{c}\hfill K_{z,\mathrm{sat}.,n}(t,z)=s_n{K}_{z,\mathrm{sat}.}(t,z).\end{array}$$(24)

If such a scattering component is not included, the perturbation produces definite holes in the (z, w)-plane and creates a spiral structure that is too clearly defined. The scattering component was introduced in order to create a smoother distribution of over- and under-densities, similar to the actual phase-space spiral of the Milky Way.

The external force we used in this work was not meant to represent a fully realistic perturbation due to a passing satellite or other source; a fully realistic perturbation of the disk can only be simulated in a complete six-dimensional phase space. Rather, the aim of the simulations was to create phase-space spirals that are qualitatively similar to the actual phase-space spiral of the Milky Way. Furthermore, the relatively low vertical speed of the passing satellite³ was chosen in order to achieve a perturbation that affects stars at lower vertical energies and also produces a perturbation that is asymmetrical (see Widrow et al. 2014 for how a satellite’s properties are correlated with the resulting disk perturbation). This gives rise to a spiral that is similar to that of the Milky Way in the following sense: While the spiral structure is not at all visible close to the origin of the (z, w)-plane, outside this innermost region the relative density of the spiral structure is close to constant in amplitude.

The main difficulty we had in reproducing the Milky Way phase-space spiral in the (z, w)-plane is related to its asymmetry. The actual phase-space spiral of the Milky Way has a single arm. Despite our efforts, we could only produce double-armed spirals, albeit with some asymmetry in terms of the relative amplitudes of the two arms. Even though the applied external force of our simulations was asymmetrical and produced clearly asymmetrical initial perturbations, the subsequent evolution of the perturbations evolved into double-armed spirals. This is further demonstrated in Appendix A, where we run a simulation with a more idealised and clearly asymmetric initial perturbation; even in this case, a second arm develops spontaneously due to self-gravity after roughly 200 Myr. Perhaps a less idealised, noisier environment is enough to quench the formation of the secondary arm, possibly as a consequence of occupying a six-dimensional phase space or other effects of dynamical diffusion. Even though our simulations did not reproduce the single-armed spiral structure, we are still confident that this test on mock data is valuable as a proof of concept and demonstration of a completely novel method for inferring the Galactic gravitational potential.

3.3. Simulations A and B

In this work we ran two separate simulations, referred to as simulation A and simulation B, on which we tested our analytical model. The parameters of these two simulations are listed in Table 2.

In terms of the initial configuration of simulations A and B, their respective matter density distributions differ, both in terms of stars and gas. For the perturbation, simulations A and B differ in all parameters. The external force applied to simulation B was weaker (mainly due to a smaller surface mass Σ) and produced a spiral with a smaller relative amplitude. The perturbing satellites’ total ranges in height are equal to σ_t × w_sat. ≃ {2, 3.6} kpc for simulations A and B, respectively.

For simplicity, we set Z_⊙ = 0 pc and W_⊙ = 0 km s⁻¹ when producing the two-dimensional data histogram d_i, j. In reality, the Sun’s vertical position and velocity are non-zero, but this would be trivial to correct for as these values are known to sufficient precision (see for example Gaia Collaboration 2020b) compared to the resolution at which we resolve the phase-space spiral in our method. In our method of inference, the parameters Z_⊙ and W_⊙ are still free parameters, which were fitted in the first step of our minimisation procedure.

4. Results

In this section we present the results of simulations A and B. Our method was applied to each simulation at two separate times. For simulation A, these times are t_A, 1 = 400 Myr and t_A, 2 = 600 Myr; for simulation B, they are t_B, 1 = 400 Myr and t_B, 2 = 500 Myr. For these instances in time, the spiral winding is roughly one complete period in the spiral region of our analytical model, where the outer boundary is defined by the mask function M(Z, W) and the inner boundary is defined by the model’s lower energy bound m(z, w | Ψ_spiral) (see Eqs. (17) and (14)). For the bulk density distribution B(z, w,  | Ψ_bulk), which is minimised in the first step of the fitting procedure, we used a Gaussian mixture model consisting of K = 6 Gaussians; fitting a higher number of Gaussians did not significantly improve the maximum likelihood.

The results of the four different cases are presented in Figs. 5–8. The figures each contain four panels that show the true and inferred gravitational potential, a two-dimensional histogram of the data, the spiral extracted directly from the data, and the best fit spiral of our analytical model. The spiral as extracted from the data, seen in panel (c), is defined as

$$\begin{array}{c}\hfill M(Z_i,W_j)\times [\frac{{\mathit{d}}_{i,j}-B(Z_i+Z_{\odot },W_j+W_{\odot }|{\mathbf{\Psi }}_{\mathrm{bulk}})}{B(Z_i+Z_{\odot },W_j+W_{\odot }|{\mathbf{\Psi }}_{\mathrm{bulk}})}].\end{array}$$(25)

Fig. 5. Data and results of simulation A at time tA, 1 = 400 Myr. The four panels show: (a) the true and inferred gravitational potentials; (b) the two-dimensional data histogram di, j in the (Z, W)-plane; (c) the phase-space spiral as extracted from the data (see Eq. (25) for the precise definition) and smoothed to an effective bin size of (40 pc)×(2 km s−1); and (d) the best fit spiral density S(Z + Z⊙, W + W⊙ | Ψspiral), as defined in Eq. (26). Panel a: the height is plotted in the range Z ∈ [ − 700, 700] pc, corresponding to the outer limit of the mask function; in all other panels, the range is Z ∈ [ − 900, 900] pc. Panels c and d: the outer boundary corresponds to that of the mask function of Eq. (17); panel d: the inner boundary corresponds to the lower limit in vertical energy for the spiral model, according to Eq. (14).

Fig. 6. Same as Fig. 5, but for simulation A at time tA, 2 = 600 Myr. (a) Gravitational potential. (b) Data. (c) Mask × (Data−Bulk)/Bulk. (d) Best fit spiral.

Fig. 7. Same as Fig. 5, but for simulation B at time tB, 1 = 400 Myr. (a) Gravitational potential. (b) Data. (c) Mask × (Data−Bulk)/Bulk. (d) Best fit spiral.

Fig. 8. Same as Fig. 5, but for simulation B at time tB, 2 = 500 Myr. (a) Gravitational potential. (b) Data. (c) Mask × (Data−Bulk)/Bulk. (d) Best fit spiral.

In panel (d), the spiral of our analytical model also includes the lower energy bound of Eq. (14), according to

$$\begin{array}{cc}& M(Z_i,W_j)\times m(Z_i+Z_{\odot },W_j+W_{\odot }|{\rho }_h)\hfill \\ \hfill & \times S(Z_i+Z_{\odot },W_j+W_{\odot }|{\mathbf{\Psi }}_{\mathrm{spiral}}).\hfill \end{array}$$(26)

We only show our inferred results in terms of the best fit, meaning the maximised likelihood. Running the minimisation algorithm is already very computationally intensive (see Sect. 2.7 for details), and computing a full Markov chain Monte Carlo would be very costly and without much return on investment because our results are dominated by systematic rather than statistical uncertainties; because of the high number of stars and the simplicity of our analytical model, the statistical uncertainty is very small (below the 1 percent level). When applying our model to the real phase-space spiral of the Milky Way, the total amount of available stars will be even higher and there will be plenty of other potential sources of systematic bias; in other words, applying this method to real data will most probably only make systematic effects dominate even more.

The inferred potential at greater heights (|z|≳600 pc) is less robustly inferred and in fact largely an extrapolation with regards to its shape at smaller heights. The main reason for this is that only a very small portion of the two-dimensional histogram, after applying the mask of Eq. (17), covers such heights. There seems to be a slight bias towards a steeper gravitational potential in the regime of larger |z|. In terms of accuracy, the most robustly inferred quantity is the gravitational potential value at heights around 400 to 500 pc, for which the relative error is only a few percent (strictly less than 5%). The precise shape of the gravitational potential close to the Galactic mid-plane is less robustly inferred; as was demonstrated in Sect. 2.3 and Fig. 1, the shape of the gravitational potential at low heights gives rise to quite small differences in the spiral’s shape.

In terms of the time that has passed since the perturbation (the parameter t), the best fit results are t_A, 1 = 432 Myr, t_A, 2 = 506 Myr, t_B, 1 = 334 Myr, and t_B, 2 = 405 Myr. The true values are 400 Myr, 600 Myr, 400 Myr, and 500 Myr, respectively, where the true time is defined in terms of when the external force of Eq. (23) was maximal. However, the force is actually active for roughly 100 Myr both before and after t = 0 Myr. Hence, we did not expect the inferred value for t to be accurate within such time frames.

In order to demonstrate the limits of our model of inference, we also show results for simulation B at times t_B, 3 = 600 Myr and t_B, 4 = 700 Myr in Appendix C. At these times, the spiral arms can no longer be seen by eye as clear continuous structures and our method starts to lose accuracy in terms of the inferred gravitational potential (the relative error is as high as 10% for the inferred potential at |z| = 400 pc). We also ran other similar simulations, for example with more complicated initial stellar phase-space distributions, for which we obtained similar results.

5. Discussion

When applying our method of inference to one-dimensional simulations, we were able to infer the gravitational potential to high accuracy. Our method extracts information from the time-varying structure of a phase-space spiral; this is a novel approach, constituting a break from traditional methods that are based on the assumption of a steady state. For this reason, our method is also complementary to such traditional methods and likely subject to different sources of systematic bias. We note that there are also other methods for extracting information from the phase-space spiral. For example, Li & Widrow (2021) infer the gravitational potential of the disk using the bulk stellar density in the (z, w)-plane and then fit the spiral to the residual of the inferred steady state distribution. The fitted spiral constitutes a consistency check of the inferred potential and is seen as a straight line in the plane of angle and frequency.

While our method is accurate in terms of one-dimensional simulations, the actual Milky Way is significantly more complex. To begin with, stars in the Galactic disk do not only oscillate in the vertical direction, but also move in the radial direction as they orbit the Galactic centre. Indeed, the phase-space spiral of the Milky Way has a different form depending on cuts in angular momentum or the galactocentric radius (see for example Fig. 14 in Laporte et al. 2019, as well as Khanna et al. 2019). When applying our method to Gaia data for the first time, we plan to select stars of the solar neighbourhood that have approximately circular orbits, for which the assumptions of our model are most likely to be valid. It would be interesting, although probably quite complicated, to see what kind of information about the global structure of the Milky Way could be extracted from the shape of the spiral in different bins of Galactic angular momentum. Ultimately, we would want to test and refine our method using a high-resolution simulation of a galaxy in three spatial dimensions. Up until now, simulations of Milky Way-like systems (such as Laporte et al. 2018 and Bland-Hawthorn & Tepper-García 2021) have not had sufficient resolution to produce well-resolved phase-space spirals in stellar density in the (z, w)-plane. To do so requires billion-particle simulations, which have only recently become feasible (Asano et al. 2020, Hunt et al. in prep.). Furthermore, our modelling of the gravitational potential needs to be expanded somewhat; the Poisson equation is no longer one-dimensional but also has a smaller contribution from the radial direction (the so-called rotation curve term).

A useful test in terms of systematic uncertainties would be to perform separate analyses of different segments of the phase-space spiral. The assumptions of our analytical model are probably only valid over a certain range in vertical energy; fitting the spiral over a large area in the z, w-plane would give a small statistical uncertainty, but potentially at the price of large systematic errors. It would be possible, at least to some extent, to quantify systematic errors by comparing the results of the smaller, separate spiral segments. Another test in a similar vein would be to abandon the form of the function $\stackrel{\sim }{\phi }(t,E_z)$ as defined in Eq. (9), which might not be completely accurate due to effects of, for example, self-gravity; instead, we could simply assume that the spiral angle $\stackrel{\sim }{\phi }$ is a smooth function with respect to E_z and model it as a Gaussian process.

Observational uncertainties are probably not a significant factor when applying this method to Gaia data and our own Galaxy. The most significant observational uncertainties are associated with the parallax measurements. With Gaia EDR3, those uncertainties are of the order of 0.02–0.03 mas for stars with absolute G-band magnitudes smaller than 15 (Lindegren et al. 2020). Within kiloparsec distances, this corresponds to a relative uncertainty of at most a few percent (meaning at most a few tens of parsecs in distance). The uncertainties with respect to velocities are typically only a few hundred metres per second (Katz et al. 2019). Neither distance- nor velocity-related uncertainties should be significant issues given the scale of the phase-space spiral and the large amount of available statistics.

Observing the phase-space spiral using Gaia requires radial velocity measurements, which are not available for all stars. Cleaning the data to exclude stars with missing radial velocities introduces a significant selection effect. However, selection effects should not be a hindrance to our method, at least not as long as the selection function is fairly smooth. What is crucial is that the shape of the phase-space spiral be extracted with accuracy; any selection effects that vary smoothly with distance will be accounted for when subtracting the bulk density component.

6. Conclusion

In this work we have developed a method for inferring the gravitational potential of the Galactic disk from the shape of a phase-space spiral in the (z, w)-plane. Our analytical model of the spiral, which is fitted to data, is based on a few simplifying assumptions, such as neglecting the self-gravity of the spiral and assuming that it evolves in a static potential. We tested our method on one-dimensional simulations (which do not adhere to the model assumptions) and were able to recover the gravitational potential with high accuracy. In cases where the phase-space spiral could be seen as a smooth continuous structure in the modelled region of the (z, w)-plane (similar to the actual spiral of the Milky Way), the relative error for the inferred potential at heights |z|≃500 pc was as small as a few percent; poor accuracy, meaning relative errors of 6–10%, was obtained only in cases where the phase-space spiral lacked these characteristics.

This article is the first in a series in which we plan to apply our method to the real Milky Way phase-space spiral, as well as test and refine our method using more complex Galaxy simulations. This is a step in the direction of modelling the time-dependent dynamics of the Milky Way and demonstrates that a time-varying structure is not necessarily an obstacle to dynamical mass measurements, but can in fact be an asset from which it is possible to extract useful information. The method in this paper completely disregards the bulk density distribution and is therefore complementary to traditional methods for which the bulk is the key quantity. Eventually, we would like to combine these separate methods (as also discussed by Li & Widrow 2021) and infer the gravitational potential of the Galactic disk using the joint information of the spiral and the bulk phase-space density distribution.

Acknowledgments

We would like to thank the anonymous referee for a useful and instructive review. AW acknowledges support from the Carlsberg Foundation via a Semper Ardens grant (CF15-0384). CL acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 852839). PFdS acknowledges support by the Vetenskapsrådet (Swedish Research Council) through contract No. 638-2013-8993 and the Oskar Klein Centre for Cosmoparticle Physics. This work made use of an HPC facility funded by a grant from VILLUM FONDEN (projectnumber 16599). This work was supported in part by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. This research utilised the following open-source Python packages: MATPLOTLIB (Hunter 2007), NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), PANDAS (McKinney et al. 2010), TENSORFLOW (Abadi et al. 2015).

References

Appendix A: Second arm forms via self-gravity

In this section we demonstrate that a phase-space spiral develops double arms via self-gravity, even when only a single arm is present to begin with. We ran a simulation whose initial phase-space configuration was the same as for simulation A (see Sect. 3.3) and then perturbed it by hand at time t = 0 by adding a kick to the vertical velocity, written as w⁺, for a subset of the simulation particles. This was done according to the formula

$$\begin{array}{c}\hfill w^{+}=\mathrm{\Theta }(z)\times \mathrm{sigm}[\frac{w-15\mathrm{km}{\mathrm{s}}^{-1}}{10\mathrm{km}{\mathrm{s}}^{-1}}]\times 10\mathrm{km}{\mathrm{s}}^{-1},\end{array}$$(A.1)

where Θ is the Heaviside step function. In this manner, only simulation particles with positive z and positive w were affected by the perturbation.

After this initial perturbation, the phase-space distribution evolves under self-gravity according to the vertical acceleration of Eq. (19). The distribution of stars in the (z, w)-plane is shown at different times in Fig. A.1. The histograms on the left-hand side, labelled N_i, j, show the stellar number density in bins with size (20 pc)×(1 km s⁻¹). The panels on the right-hand side show the difference between the stellar number density histogram and the average value of its neighbouring area; this is written as $N_{i,j}-{\overline{N}}_{i,j}$, where ${\overline{N}}_{i,j}$ is equal to N_i, j convoluted with a bivariate Gaussian distribution with standard deviations of 80 pc and 4 km s⁻¹. For better visibility, the panels on the right-hand side are illustrated with a resolution of (40 pc)×(2 km s⁻¹).

Fig. A.1. Stellar number density distribution in the (z, w)-plane shown at times t = {0, 50, 100, 200, 400} Myr (going from top to bottom) for a simulation affected by an asymmetric perturbation at time t = 0 Myr. Left and right panels: (a) a histogram of the total stellar density with bin size (20 pc)×(1 km s−1); and (b) the difference between the stellar density histogram and the average of the surrounding histogram bins. All panels have the same coordinate axis. The colour scheme in the panels on the right-hand side is not linear, but in arcsinh scale.

Even though the initial perturbation is asymmetric and only a single spiral arm is visible at t = 100 Myr, a second arm develops; it is visible at t = 200 Myr and becomes even more visible at t = 400 Myr. This indicates that a symmetric spiral component (i.e. a double-armed spiral) develops spontaneously as an effect of self-gravity in our one-dimensional simulations. We tried many other initial simulation conditions, for example: varying the matter components’ respective scale heights and mid-plane densities; varying the satellite perturbation’s properties (most importantly, the vertical speed) in order to perturb stars with either higher or lower vertical energies and to produce more of a bending or breathing mode; and trying a variety of more idealised and strictly asymmetrical perturbations (similar to the one shown in Fig. A.1). Despite these efforts, we always found similar results, which tended towards a symmetric phase-space spiral.

Appendix B: Inferred vertical acceleration and total matter density

In Fig. B.1 we show the inferred vertical acceleration (K_z = −∂Φ/∂Z) and the inferred total matter density (ρ, fulfilling the Poisson equation; see Eq. (7)). These quantities are proportional to the first and second derivatives of the gravitational potential. For this reason, they are less robustly inferred, and the matter density is especially sensitive to the precise shape of the inferred gravitational potential.

Fig. B.1. Inferred vertical acceleration, Kz(Z), and total matter density, ρ(Z), for simulations A and B. All panels share the same horizontal axis.

In the inferred matter density profiles, the scale height, shape, and mid-plane matter density are not very accurately inferred quantities. However, the gravitational acceleration for heights |Z| in the range of 250–300 pc (or, equivalently, the total surface matter density integrated over heights |Z|< 300 pc) is accurately inferred for all four examples, with a relative error strictly smaller than 4%. The gravitational acceleration and matter density at greater heights (|Z|≳400 pc) are not very accurately inferred, which is in large part due to the small amount of statistics at these heights (connected to how the gravitational potential is poorly inferred for |Z|≳600 pc; see Sect. 4 for further details).

Appendix C: Further results for simulation B

In this section we present some further results for simulation B, at times t_B, 3 = 600 Myr and t_B, 4 = 700 Myr. We do so in order to demonstrate the limits of our method. Simulation B was subjected to a relatively weak perturbation, and the resulting spiral has a low amplitude (roughly 10% in terms of relative stellar number density, compared to roughly 20% for simulation A). With time and further winding, the spiral structure becomes all the more phase-mixed and increasingly hard to detect, eventually disappearing completely. In Figs. C.1 and C.2, we show results in the interim, where the spiral begins to be more difficult to see.

Fig. C.1. Same as Fig. 5, but for simulation B at time tB, 3 = 600 Myr. (a) Gravitational potential. (b) Data. (c) Mask × (Data−Bulk)/Bulk. (d) Best fit spiral.

Fig. C.2. Same as Fig. 5, but for simulation B at time tB, 4 = 700 Myr. (a) Gravitational potential. (b) Data. (c) Mask × (Data−Bulk)/Bulk. (d) Best fit spiral.

At these later times, when the spiral structure is less pronounced, our method loses accuracy. In terms of the gravitational potential at heights in the range 400 to 500 pc, the relative error is roughly 6% for t_B, 3 and 10% for t_B, 4. The inferred times since the perturbation are t_B, 3 = 451 Myr and t_B, 4 = 725 Myr (true values are 600 Myr and 700 Myr).

All Tables

All Figures

Fig. 1. Phase-space spiral represented as a line graph, shown in solid blue. Two points, labelled A and B, are connected to the origin point via dashed black lines. Points A and B are both associated with a total vertical energy: Ez, A = w2/2 is given by the vertical velocity where the spiral passes through the Galactic mid-plane; and Ez, B = Φ(z) is given by gravitational potential at a point of the spiral with zero velocity. In the limit of high winding, the difference Ez, A − Ez, B would approach zero from above.

In the text

Fig. 2. Three phase-space spirals (labelled A, B, and C) inhabiting three different gravitational potentials. Top panel: vertical force per mass (Kz = −∂Φ/∂z) as a function of height for the respective gravitational potentials. Bottom panel: three spirals in the (z, w)-plane, where their winding time (t) and initial angle (${\stackrel{\sim }{\phi }}_0$) are chosen such that the spirals have the same phase ($\stackrel{\sim }{\phi }$) and phase derivative ($\partial \stackrel{\sim }{\phi }/\partial z_{\max }$) when their vertical energy is equal to Ez = Φ(800 pc).

In the text

	Fig. 3. Mask function M(Z, W) in the (Z, W)-plane, with boundary values Zlim. = 700 pc and Wlim. = 40 km s−1. See Eq. (17) for details.
In the text

	Fig. 4. Example of the matter density distributions from which the particles of our one-dimensional simulations are drawn. In this plot, the distribution corresponds to that of simulation A; see Sect. 3 and Table 2 for details.
In the text

Fig. 5. Data and results of simulation A at time tA, 1 = 400 Myr. The four panels show: (a) the true and inferred gravitational potentials; (b) the two-dimensional data histogram di, j in the (Z, W)-plane; (c) the phase-space spiral as extracted from the data (see Eq. (25) for the precise definition) and smoothed to an effective bin size of (40 pc)×(2 km s−1); and (d) the best fit spiral density S(Z + Z⊙, W + W⊙ | Ψspiral), as defined in Eq. (26). Panel a: the height is plotted in the range Z ∈ [ − 700, 700] pc, corresponding to the outer limit of the mask function; in all other panels, the range is Z ∈ [ − 900, 900] pc. Panels c and d: the outer boundary corresponds to that of the mask function of Eq. (17); panel d: the inner boundary corresponds to the lower limit in vertical energy for the spiral model, according to Eq. (14).

In the text

	Fig. 6. Same as Fig. 5, but for simulation A at time tA, 2 = 600 Myr. (a) Gravitational potential. (b) Data. (c) Mask × (Data−Bulk)/Bulk. (d) Best fit spiral.
In the text

	Fig. 7. Same as Fig. 5, but for simulation B at time tB, 1 = 400 Myr. (a) Gravitational potential. (b) Data. (c) Mask × (Data−Bulk)/Bulk. (d) Best fit spiral.
In the text

	Fig. 8. Same as Fig. 5, but for simulation B at time tB, 2 = 500 Myr. (a) Gravitational potential. (b) Data. (c) Mask × (Data−Bulk)/Bulk. (d) Best fit spiral.
In the text

Fig. A.1. Stellar number density distribution in the (z, w)-plane shown at times t = {0, 50, 100, 200, 400} Myr (going from top to bottom) for a simulation affected by an asymmetric perturbation at time t = 0 Myr. Left and right panels: (a) a histogram of the total stellar density with bin size (20 pc)×(1 km s−1); and (b) the difference between the stellar density histogram and the average of the surrounding histogram bins. All panels have the same coordinate axis. The colour scheme in the panels on the right-hand side is not linear, but in arcsinh scale.

In the text

	Fig. B.1. Inferred vertical acceleration, Kz(Z), and total matter density, ρ(Z), for simulations A and B. All panels share the same horizontal axis.
In the text

	Fig. C.1. Same as Fig. 5, but for simulation B at time tB, 3 = 600 Myr. (a) Gravitational potential. (b) Data. (c) Mask × (Data−Bulk)/Bulk. (d) Best fit spiral.
In the text

	Fig. C.2. Same as Fig. 5, but for simulation B at time tB, 4 = 700 Myr. (a) Gravitational potential. (b) Data. (c) Mask × (Data−Bulk)/Bulk. (d) Best fit spiral.
In the text