Skip to Main Content

Article Navigation

Journal Article

Hamiltonians of spherical Galaxies in action-angle coordinates

Abstract

We present a simple formula for the Hamiltonian in terms of the actions for spherically symmetric, scale-free potentials. The Hamiltonian is a power law or logarithmic function of a linear combination of the actions. Our expression reduces to the well-known results for the familiar cases of the harmonic oscillator and the Kepler potential. For other power laws, as well as for the singular isothermal sphere, it is exact for the radial and circular orbits, and very accurate for general orbits. Numerical tests show that the errors are always very small, with mean errors across a grid of actions always <1 per cent and maximum errors <2.5 per cent. Simple first-order corrections can reduce mean errors to <0.6 per cent and maximum errors to <1 per cent. We use our new result to show that: (1) the misalignment angle between debris in a stream and a progenitor is always very nearly zero in spherical scale-free potentials, demonstrating that streams can sometimes be well approximated by orbits, (2) the effects of an adiabatic change in the stellar density profile in the inner region of a galaxy weaken any existing 1/r dark matter density cusp, which is reduced to 1/r^1/3. More generally, we derive the full range of adiabatic cusp transformations and show that a |$1/r^{\gamma _0}$| density cusp may be changed to |$1/r^{\gamma _1}$| only if γ₀/(4 − γ₀) ≤ γ₁ ≤ (9 − 2γ₀)/(4 − γ₀). It follows that adiabatic transformations can never completely erase a dark matter cusp.

methods: analytical, Galaxy: kinematics and dynamics, galaxies: kinematics and dynamics

1 INTRODUCTION

Action–angle coordinates are an extremely useful set of coordinates with which to describe a stellar system. Actions are integrals of motion given by

\begin{equation} J_{i} = \frac{1}{2\pi }\oint _{\gamma _{i}} \boldsymbol {p}\,\mathrm{d}\boldsymbol {q}, \end{equation}

(1)

where γ_i is a closed path around an orbital torus, and |$(\boldsymbol {p},\boldsymbol {q})$| are canonically conjugate phase-space coordinates. Since actions are integrals of motion, they satisfy |$\dot{J_{i}}=0$|⁠, greatly simplifying Hamilton's equations for the canonically conjugate angle variables. These follow the trivial time evolution

\begin{equation} \theta _{i}(t)=\theta _{i}(0)+\Omega _{i}t,\qquad \qquad \Omega _{i} = \frac{\mathrm{\partial} H}{\mathrm{\partial} J_{i}}. \end{equation}

(2)

Already, we can see the appeal of the action–angle coordinate space. In the physical phase space |$(\boldsymbol {x},\boldsymbol {v})$|⁠, each of the six phase-space coordinates generally evolves in a complicated fashion. In |$(\boldsymbol {J},\boldsymbol {\theta })$| space, three of the coordinates are constant and the other three evolve linearly with time. Action–angle space is the natural arena in which to view stellar orbits.

The benefit of action–angles extends beyond their simple time evolution, however. Despite the fact that galaxies are never truly in a steady state, equilibrium models are hugely important tools in developing our understanding of galaxies. Using equilibrium models, we can infer the gravitational potential of a galaxy and predict all observable quantities. Under most circumstances, we are justified in assuming that a galaxy is always close to an equilibrium state, and we can model non-equilibrium effects such as accretion or stellar evolution as small perturbations to the system.

Given the assumption of dynamical equilibrium, the most efficient way to proceed is via the use of the Jeans theorem (see Binney & Tremaine 2008), which tells us that the phase-space distribution function (hereafter DF) of a galaxy or stellar system in equilibrium may be written as a function of just three isolating integrals. This reduces the dimensionality of the space under consideration from six to three, and so invoking the Jeans theorem leads to an elegant simplification of the modelling problem. Since actions are adiabatic invariants, if we choose to write the DF as a function of the actions, its functional form is preserved under slow changes to the galaxy. This property can be exploited to model various processes in galaxy evolution. For example, Goodman & Binney (1984) showed that the adiabatic growth of mass at the centre of a stellar system causes the velocity dispersion to become tangentially distended, a process invoked by Agnello et al. (2014) to explain the kinematics of M87's globular cluster populations. Similarly, Katz et al. (2014) recently modelled the adiabatic compression of dark matter haloes in galaxies from the THINGS survey to reproduce their rotation curves. Another example is provided by Gondolo & Silk (1999), who showed that the adiabatic growth of a black hole can cause a dark matter spike at the centre of a cusped dark halo.

Unfortunately, although action–angles are very useful when known, they are not usually available analytically. In fact, Evans, de Zeeuw & Lynden-Bell (1990) showed that the most general potential in which the actions are available as elementary functions is the isochrone (which reduces to the Kepler and harmonic potentials in limiting cases). Recent efforts by Binney (2012), Bovy (2014) and Sanders & Binney (2014) mean that actions can now be numerically recovered in many potentials, but analytical approximations are still sparse. Though numerical methods are ultimately the most widely applicable, analytical approximations are important in so far as they give us valuable insight into the properties of a system, not to mention drastically reduced computational time when they are sufficiently accurate.

In this paper, we present a new method for finding accurate approximations to the Hamiltonians of spherically symmetric, scale-free systems and apply it to general power-law potentials and the isothermal sphere. We then use the method in two applications: finding the misalignment angle of a stellar stream, and calculating the effect of a slowly changing central stellar density profile on the slope of a dark matter cusp.

2 METHOD

In a spherical potential, the radial action is given by

\begin{equation} J_{r} = \frac{1}{2\pi }\oint \sqrt{2E-2\Phi (r)-L^{2}/r^{2}}\,\mathrm{d}r. \end{equation}

(3)

The angular actions are related to the angular momentum components via J_ϕ = |L_z| and J_θ = L − |L_z|, where L = J_θ + J_ϕ is the total angular momentum and L_z is the angular momentum about the z-axis. The Hamiltonian can then be regarded as a function of just two quantities: H ≡ H(L, J_r). The most general method of finding H(L, J_r) is to solve the integral equation (3) for E but, as previously mentioned, this cannot be done for the general case.

We now present a prescription for finding approximations |$\mathcal {H}(L,J_{r})$| to the exact Hamiltonian H(L, J_r) of a spherical, scale-free potential Φ(r). In terms of the usual phase-space coordinates, the Hamiltonian is given by

\begin{equation} H(r,\boldsymbol {v}) = \frac{1}{2}\boldsymbol {v}^{2} + \Phi (r). \end{equation}

(4)

We now look to find an approximate transformation by considering only the two most basic types of orbit: radial and circular orbits. For the radial orbits, we must evaluate

\begin{equation} J_{r} = \frac{1}{2\pi }\oint \sqrt{2E(0,J_{r})-2\Phi (r)}\,\mathrm{d}r \end{equation}

(5)

and invert the resulting expression for E(0, J_r). The energy of the circular orbits is given by the system of simultaneous equations

\begin{equation} L(r) = rv_{\rm c}(r), \qquad \qquad E(L,0) = \frac{1}{2}v_{\rm c}(r)^{2} + \Phi (r), \end{equation}

(6)

where v_c(r) is the circular speed. For a scale-free potential, simply as a consequence of dimensional analysis, both results will have the same functional form

\begin{equation} E(J_{i},J_{j}=0) = f(a_{i}J_{i}), \end{equation}

(7)

where f (x) is some function and a_i is a constant. This means that a natural choice of interpolation for a general orbit is

\begin{equation} \mathcal {H}(L,J_{r}) = f(aL+bJ_{r}), \end{equation}

(8)

where a and b are found using the solutions to equations (5) and (6). Therefore, as long as a and b are available, one has a very simple analytical approximation to the true Hamiltonian of any scale-free spherical system.

If our Hamiltonian is not scale free, we can often break the problem up into scale-free regimes, and devise an interpolation formula that smoothly matches the limits. A specific example of this is provided by Evans & Williams (2014), who provided a Hamiltonian |$\mathcal {H}(L,J_r)$| for a halo model with a 1/r density cusp at the centre but a flat rotation curve at large radii.

3 POWER-LAW POTENTIALS

3.1 Results

Let us now consider power-law potentials of the form

\begin{equation} \Phi (r) = Ar^{\alpha },\qquad A = {v_0^2\over \alpha r_0^\alpha }. \end{equation}

(9)

The choice of constant A is motivated by the rotation law, which takes the form

\begin{equation} v_{\rm c}^2(r) = v_0^2\left({r\over r_0}\right)^{\alpha }. \end{equation}

(10)

So, the rotation curve has the value v₀ at reference radius r₀. Models with α > 0 (or α < 0) have a rising (or falling) rotation curve. The model with α = 0 has a flat rotation curve, and is recognized as the singular isothermal sphere. To recover its familiar logarithmic potential, we must add a constant |$-v_0^2/\alpha$| to equation (9) and use L'Hôpital's rule to take carefully the limit α → 0. The physical range of α is −1 ≤ α ≤ 2, with the extremes representing a Keplerian potential (point mass) and a harmonic potential (homogenous matter distribution), respectively.

Applying our method to case when A, α > 0 (the rising rotation curve case), we find

\begin{eqnarray} E(0,J_{r})&=&\bigg [\frac{\sqrt{2\pi }\alpha A^{1/\alpha }\Gamma (3/2+ 1/\alpha )}{\Gamma (1/\alpha )}\bigg ]^{\beta }J_{r}^{\beta }, \nonumber \\ E(L,0)&=&\frac{(A\alpha )^{\beta /\alpha }}{\beta } L^\beta , \end{eqnarray}

(11)

where we have set β = 2α/(α + 2). This implies that f(x) ∝ x^β. This gives us the following very simple approximation for the Hamiltonian

\begin{equation} \mathcal {H}(L,J_{r}) = C(L + DJ_{r})^{\beta }, \end{equation}

(12)

where the constants C and D are

\begin{eqnarray} C &=& \frac{v_0^{2\beta /\alpha }}{r_0^\beta \beta }, \nonumber \\ D &=& \frac{\sqrt{2\pi }\Gamma (3/2+ 1/\alpha )\alpha ^{-1/\alpha }\beta ^{1/\beta }}{\Gamma (1+1/\alpha )}. \end{eqnarray}

(13)

This prescription may be tested on the special case of the harmonic oscillator (α = 2). Setting |$A=\frac{1}{2}\Omega ^{2}$|⁠, we recover the well-known exact analytical result (see, for example, Goldstein 1980)

\begin{equation} \mathcal {H}(L,J_{r}) = 2\Omega J_{r} + \Omega L. \end{equation}

(14)

So, |$\mathcal {H}$| perfectly coincides with the true Hamiltonian H for the harmonic oscillator.

When A, α < 0 (the falling rotation curve case), the result is the same as equation (12), except

\begin{equation} D = \frac{\sqrt{2\pi }\Gamma (1-1/\alpha )(-\alpha )^{1-1/\alpha }(-\beta )^{1/\beta }}{\Gamma (-1/\alpha -1/2)}. \end{equation}

(15)

The obvious test in this case is the Keplerian potential, where A = −GM and α = −1, which gives

\begin{equation} \mathcal {H}(L,J_{r}) = -\frac{1}{2}\bigg (\frac{GM}{L+J_{r}}\bigg )^{2}. \end{equation}

(16)

This coincides with the exact result for the Kepler potential, derived by the residue calculus in Born (1927) and repeated in Goldstein (1980).

The case of the flat rotation curve (or the singular isothermal sphere) can be derived by taking the careful limit α → 0 and using L'Hopital's rule. The potential becomes logarithmic |$\Phi =v_{0}^{2}\log r/r_0$|⁠, and so we find that f(x) ∝ log (x). Our approximation then gives

\begin{equation} \mathcal {H}(L,J_{r}) = v_{0}^{2}\log \bigg (\frac{\sqrt{e}L+\sqrt{2\pi }J_{r}}{v_{0}r_0}\bigg ), \end{equation}

(17)

as stated in Evans & Williams (2014).

So, we have found a surprisingly simple and compact formula (12) for the Hamiltonian in terms of the actions valid for any power-law potential with −1 ≤ α ≤ 2. It is exact for the two classical cases – the Kepler potential and the harmonic oscillator – which have been known for over a century. Note too that the frequencies Ω_i are rationally related only in the case of the Kepler potential and the harmonic oscillator. Thus, we have recovered Bertrand's theorem (see e.g. Goldstein 1980), which states that the only central force laws with everywhere closed bound orbits are the linear and inverse square force laws. For all other cases, our formula (12) gives frequencies that are incommensurable and so the motion is conditionally periodic.

3.2 Errors and first-order corrections

Given that the approximation is only exact for the cases α = −1 and α = 2, we tested the effectiveness of the approximation for several values of α: −0.5, 0 (the isothermal sphere), 0.5, 1 and 1.5. In each case, we evaluate a set of 10 000 orbits on a grid of pericentre and apocentre radii: r_p = 0 → 50, r_a = 0 → 100. For each orbit, we evaluate the energy and angular momentum, which are available as algebraic expressions in terms of pericentre and apocentre radii (Lynden-Bell 2010), and the radial action, which requires the evaluation of a numerical integral. This results in an irregularly spaced grid in action space.

Before discussing numerical values, it is necessary to explain how we define errors across the sequence of models. We consider correction terms of the form

\begin{equation} \mathcal {H} = f(aL+bJ_{r}+\epsilon \sqrt{LJ_{r}}), \end{equation}

(18)

where f(x) ∝ x^β (if β ≠ 0) or f(x) ∝ log x (if β = 0). Due to lack of dimensional constants in the problem, the correction must possess the dimensions of actions. Its form is suggested by examing the error distribution along rays in action space. A measure of the fractional error between the corrected values and the original approximation is then

\begin{equation} \frac{\Delta \mathcal {H}^{1/\beta }}{\mathcal {H}^{1/\beta }} =\frac{\epsilon \sqrt{LJ_{r}}}{aL+bJ_{r}+\epsilon \sqrt{LJ_r}} \approx \frac{\epsilon \sqrt{LJ_{r}}}{aL+bJ_{r}}. \end{equation}

(19)

In the isothermal case (β = 0), this becomes

\begin{equation} \frac{\Delta \mathcal {\exp H_{\mathrm{iso}}}}{\exp \mathcal {H}_{\mathrm{iso}}} \approx \frac{\epsilon \sqrt{LJ_{r}}}{aL+bJ_{r}} \end{equation}

(20)

and so the measure of error is a continuous function of power-law index α.

We find that the errors in all cases are very small. The mean absolute percentage error on the action grid is <1 per cent in every case, and the maximum absolute percentage error is always < 2.5 per cent. These results are given in Fig. 1. The distribution of errors is very similar in all cases: it is scale free and maximized along a particular ray in action space. A typical example is given in Fig. 2.

Absolute percentage errors against power-law index (α = 0 corresponds to the isothermal sphere). The dashed lines correspond to the original approximation, and the solid lines include first-order corrections. Blue lines are maximum absolute percentage errors across the action grid, red lines are the mean absolute percentage errors. One can see that the first-order corrections reduce the maximum error significantly, and also have a noticeable effect on the mean error.

Figure 1.

Absolute percentage errors against power-law index (α = 0 corresponds to the isothermal sphere). The dashed lines correspond to the original approximation, and the solid lines include first-order corrections. Blue lines are maximum absolute percentage errors across the action grid, red lines are the mean absolute percentage errors. One can see that the first-order corrections reduce the maximum error significantly, and also have a noticeable effect on the mean error.

Open in new tab Download slide

Distribution of percentage error with actions for the case α = 1.5. One can see a definite systematic trend that could be removed with an astute choice of first-order correction.

Figure 2.

Distribution of percentage error with actions for the case α = 1.5. One can see a definite systematic trend that could be removed with an astute choice of first-order correction.

Open in new tab Download slide

Given the simplicity of the error distribution, we can calculate a value for ϵ in equation (18) very easily. Assuming that the errors are maximized along the line L = χ²J_r, we can substitute this into equations (19) or (20) to give

\begin{equation} \epsilon = \frac{\delta _{\mathrm{max}}(\chi ^2a+b)}{\chi }, \end{equation}

(21)

where δ_max is the fractional error with the largest absolute value on the grid. We carried out this analysis for each of our test cases and found that it had a significant effect on the maximum error on the grid, as well as reducing the average error in every case. These results are also given in Fig. 1. If these first-order corrections are included, the error never exceeds 1 per cent for any of our test cases.

On account of its astrophysical importance, we give the expression for the first-order corrected formula for the singular isothermal sphere

\begin{eqnarray} \mathcal {H} &=& v_{0}^2\log \bigg (\frac{\sqrt{e}L+\sqrt{2\pi }J_{r}+\epsilon \sqrt{LJ_{r}}}{v_{0}r_0}\bigg ), \nonumber \\ \epsilon &=& -0.1009. \end{eqnarray}

(22)

For the case α = 1, which corresponds to the cosmologically motivated 1/r or Navarro, Frenk & White (1996, NFW) density cusp, the corrected expression is

\begin{eqnarray} \mathcal {H} &=& {3v_0^{4/3}\over 2r_0^{2/3}} \bigg (L+ {\pi \over \sqrt{3}} J_{r}+\epsilon \sqrt{LJ_{r}}\bigg )^{2/3}, \nonumber \\ \epsilon &=& -0.0412. \end{eqnarray}

(23)

4 APPLICATIONS

Here, we present two simple applications of our results. The first is to stellar streams. We demonstrate that within this approximation streams are well delineated by orbits. The second is to cusps in the centre of galaxies. We show that the central slope of a dark matter halo is reduced when the stellar populations begin to dominate the potential.

4.1 Stellar streams

Sanders & Binney (2013) showed that a stream can only be successfully modelled by an orbit if the misalignment angle is very small. This angle is given by

\begin{equation} \varphi = \arccos (\hat{\boldsymbol {\Omega }}.\hat{\boldsymbol {e}}_{1}), \end{equation}

(24)

where |$\hat{\boldsymbol {\Omega }}$| is the normalized frequency vector of the progenitor of the stellar stream and |$\hat{\boldsymbol {e}}_{1}$| is the leading eigenvector of the Hessian matrix

\begin{equation} D_{ij}(\hat{J})=\frac{\mathrm{\partial} ^{2}H}{\mathrm{\partial} J_{i}\mathrm{\partial} J_{j}}. \end{equation}

(25)

The misalignment angle quantifies the offset in angle space between stripped material and the progenitor, and encodes the direction in which stars spread in action–angle space. For a long, thin stream to form, the Hessian must have one dominant eigenvalue. We now evaluate D_ij within our approximate framework. Any approximate Hamiltonian produced by our method is of the form

\begin{equation} \mathcal {H}(L,J_{r})=f(aL+bJ_{r}) \end{equation}

(26)

giving rise to a Hessian matrix

\begin{equation} D_{ij}(\hat{J}) = G(L,J_{r})\left(\begin{array}{cc}a^{2} &\quad ab\\ ab &\quad b^{2}\end{array}\right). \end{equation}

(27)

Such a matrix has only one non-zero eigenvalue, λ = G(a² + b²), with corresponding eigenvector

\begin{equation} \hat{\boldsymbol {e}}_{1}=\frac{1}{\sqrt{a^2/b^2+1}}\left(\begin{array}{c}a/b\\ 1\end{array}\right).\ \end{equation}

(28)

Now, if we evaluate a general frequency vector using |$\mathcal {H}$|⁠, we find

\begin{equation} \hat{\boldsymbol {\Omega }}\propto \hat{\boldsymbol {e}}_{1}. \end{equation}

(29)

This implies that the misalignment angle vanishes, and so the stream will always develop along an orbit under this approximation. This suggests that, at least in the scale-free spherical case, the misalignment angle must be very small and therefore the stream should be very well approximated by an orbit. This provides a measure of justification for the assumption used in Koposov, Rix & Hogg (2010), in which the GD1 stream (Grillmair & Dionatos 2006) was assumed to lie along an orbit in a scale-free, mildly axisymmetric potential.

4.2 Adiabatic transformation of density cusps

Navarro et al. (1996) famously discovered that, in dissipationless simulations, dark haloes form a 1/r density cusp in the centre of galaxies. However, there is a discrepancy between the predicted density of dark matter at the centre of spiral galaxies and the inferred dark matter density consistent with observations. For example, in the case of the Milky Way, Binney & Evans (2001) show that cuspy dark matter density haloes would violate constraints on the rotation curve, given the known contributions from the stellar and ISM discs and Galactic bar. Clearly, the 1/r cusp, if it originally formed in the Milky Way galaxy, has been weakened with the passage of time.

We assume that the dark halo of the galaxy formed first, and that the dark matter density in the centre of a galaxy originally followed a 1/r cusp. With the build-up of the disc and bulge stellar populations, we assume that the density profile of the dominant stellar populations becomes roughly constant in the very central parts. In this instance, the potential is initially linear (α₀ = 1), corresponding to the NFW cusp, and becomes harmonic, (α₁ = 2), corresponding to a constant density core. Assuming that this process occurs sufficiently slowly that the adiabatic invariance of the actions holds good, then we can find the final profile of the dark matter.

For the collisionless dark matter, we use the self-consistent power-law DF derived by Evans (1994, equation 5.4), and assume an initially isotropic population so that

\begin{equation} f(E) = \mathcal {N}E^{- \frac{\alpha _{0} + 4}{2\alpha _{0}}}, \end{equation}

(30)

where |$\mathcal {N}$| is a normalization factor. We then insert our approximation for H(L, J_r) from equation (13) for the case α = α₀ = 1 to find

\begin{equation} f(L,J_{r}) \simeq \mathcal {M}\bigg (L+\frac{\pi J_{r}}{\sqrt{3}}\bigg )^{-5/3}, \end{equation}

(31)

where |$\mathcal {M}$| is a new normalization factor. Since the changes in the potential are presumed to occur adiabatically, the DF is always of this form, even at the endpoint. We now eliminate J_r by again using equation (13) but for the final potential, |$\Phi _{1}=\frac{1}{2}\Omega ^{2}r^2$| (which is exact in this case), giving

\begin{equation} f(L,J_{r})\propto \bigg (L + \frac{\pi }{2\sqrt{3}}\bigg [\frac{H}{\Omega } - L\big )\bigg ]\bigg )^{-5/3} \end{equation}

(32)

and substitute for H using equation (4), which finally leaves

\begin{eqnarray} f(L,J_{r})\propto \bigg (\big (1-\frac{\pi }{2\sqrt{3}}\big )L + \frac{\pi }{2\Omega \sqrt{3}}\big (\frac{1}{2}\boldsymbol {v}^{2} + \frac{1}{2}\Omega ^{2}r^{2}\big )\bigg )^{-5/3}.\nonumber \\ \end{eqnarray}

(33)

We can simplify this rather messy looking expression. The first term is proportional to L, which will be small for particles that penetrate the centre of the galaxy. This term also has a considerably smaller numerical pre-factor than the isotropic term. On this basis, we consider only the second term and integrate over velocity space to find the spatial dependence of the final density profile

\begin{equation} \rho _{\mathrm{dm}}(r) \propto \int _{0}^{\infty } \frac{v^2\mathrm{d}v}{\left (v^{2} + \Omega ^{2}r^{2}\right )^{5/3}}. \end{equation}

(34)

This integral can be evaluated using the substitution v² = Ω²r²tan ²θ, and gives the final result

\begin{equation} \rho _{\mathrm{dm}}(r) \propto r^{-1/3}. \end{equation}

(35)

So, we have found that the density profile of dark matter particles that initially followed an r⁻¹ density cusp has developed a shallower density profile of r^−1/3 as the density becomes homogeneous in the very central parts.

This calculation can be carried out completely, generally. Suppose, the initial dark matter cusp has form |$r^{-\gamma _0}$|⁠. Let the gravitational potential evolve so that the final end-state has |$\phi \propto r^{\alpha _{1}}$|⁠. Then the final dark matter density cusp behaves like |$r^{-\gamma _1}$|⁠, where

\begin{equation} \gamma _1(\gamma _0,\alpha _1) = {\gamma _0 - \alpha _1(\gamma _0-3) -6\over \gamma _0 -4}. \end{equation}

(36)

Our earlier result (35) corresponds to the special case of an NFW cusp responding to a potential that slowly changes to harmonic, that is γ₁(1, 2) = 1/3.

This equation enables us to map out the domain of accessible adiabatic changes of cusp slope, as shown in Fig. 3. The upper line corresponds to γ₁(γ₀, −1), in which the underlying potential is slowly changed to Keplerian and the matter distribution becomes strongly concentrated. The lower line corresponds to γ₁(γ₀, 2), as the potential is changed to harmonic and the matter distribution becomes homogeneous. The starting cusp γ₀ may be adiabatically transformed into γ₁ provided

\begin{equation} {\gamma _0 \over 4 - \gamma _0} \le \gamma _1 \le {9- 2\gamma _0 \over 4-\gamma _0}. \end{equation}

(37)

This region is shown unshaded in Fig. 3. It is the domain of accessible cusp transformations. Note that an initial NFW cusp can only ever be adiabatically transformed into a cusp with slope 1/3 ≤ γ₁ ≤ 7/3. It can never be completely erased. Of course, this result is restricted to adiabatic changes. There are a variety of non-adiabatic processes that are widely believed to remove cusps, including repeated mass-loss (Read & Gilmore 2005), AGN-driven outflows (Martizzi, Teyssier & Moore 2013) and infall of satellites (El-Zant, Shlosman & Hoffman 2001).

The domain of allowed cusp slope transformations. The final cusp slope γ1 is plotted against the initial cusp slope γ0. Shaded regions cannot be reached by adiabatic changes. The upper line corresponds to γ1(γ0, −1) (the Keplerian limit) and the lower line to γ1(γ0, 2) (the harmonic limit).

Figure 3.

The domain of allowed cusp slope transformations. The final cusp slope γ₁ is plotted against the initial cusp slope γ₀. Shaded regions cannot be reached by adiabatic changes. The upper line corresponds to γ₁(γ₀, −1) (the Keplerian limit) and the lower line to γ₁(γ₀, 2) (the harmonic limit).

Open in new tab Download slide

Finally, we note that Gondolo & Silk (1999), in their study of the response of a dark matter cusp to a growing black hole, derived the limit γ₁(γ₀, −1), though only for models with 0 < γ₀ < 2. The formula is generally valid throughout the entire physical range 0 ≤ γ₀ ≤ 3, though the methods employed by Gondolo & Silk (1999) were unable to include isothermal cusp slopes and steeper.

5 CONCLUSIONS

Action–angle coordinates are a powerful tool in modern classical mechanics (see e.g. Abraham & Marsden 1978; Arnold 1989) and galactic dynamics (Binney & Tremaine 2008). Their widespread use is hampered by the fact that the actions can be difficult to find explicitly. Although numerical tools are becoming available to make this job easier (Binney 2012; Bovy 2014; Sanders & Binney 2014), there are disappointingly few potentials for which analytic results are possible.

The main result of the paper is the discovery of a simple formula for the Hamiltonian of scale-free potentials in terms of the actions. The formula recovers the well-known and exact results for the harmonic oscillator and Keplerian potentials, which have been known for over a century (see e.g. Goldstein 1980). For other power laws, the formula is exact for radial and circular orbits, but approximate for general orbits. However, numerical tests show that the errors are always very small, with mean errors across a grid of actions always <1 per cent, and maximum errors never larger than 2.5 per cent. Our formula therefore provides a useful and accurate approximate Hamiltonian for many of the commonly used galactic potentials, including the isothermal sphere and power-law cusps. It may also be improved with simple first-order corrections if need be, which reduce the mean error to always <0.6 per cent, and the maximum error <1 per cent.

We provided two applications. First, the mismatch of streams from their progenitor is described by a misalignment angle. This in turn is controlled by the Hessian of the Hamiltonian with respect to the actions (Sanders & Binney 2013). Most streams in the Milky Way halo occur at Galactocentric radii greater than 20 kpc, so that the potential is probably scale free to a good approximation. If, in addition, the potential is nearly spherical, then we have shown that the misalignment angle must be close to zero. This suggests that there are regimes in which the approximation of thin streams by orbits may be a reasonable one.

Secondly, we considered the evolution of a dark matter 1/r density cusp in a spiral galaxy. Although 1/r density cusps are known to be the endpoint of dissipationless simulations, there is ample evidence that in galaxies like the Milky Way, any such cusp must have been modified at the centre (Binney & Evans 2001). The response of a 1/r cusp to the slow build-up of a stellar bulge and disc can be computed using our formula, under the assumption of adiabatic invariance. As the overall potential becomes cored and harmonic, the dark matter density cusp becomes shallower and is nearly erased. It began as a 1/r and finished as 1/r^1/3 density cusp. Considering the more general problem of density cusps of arbitrary slope, we show that adiabatic transformations can only weaken, and never remove, such cusps. A density singularity – once present – requires some non-adiabatic process to erase it.

In this paper, we only considered scale-free potentials. Nonetheless, our methods can also be applied to spherical potentials that possess a scalelength. This requires more work because we have to devise a suitable interpolation scheme to sew the different regimes together, but it is still very effective. Evans & Williams (2014) provide a specific example, but it would be interesting to study this problem in greater generality.

Also looking to the future, it is worthwhile to consider whether similar reasoning could be applied to scale-free, axisymmetric potentials. The regularity of Poincaré surfaces of sections for axisymmetric scale-free potentials (Richstone 1982; Evans 1994) strongly hints that there are simple results to be found. This is a more complex problem, given that there are now three actions to consider rather than the two in the spherical case (J_r and J_θ + J_ϕ). Nonetheless, we can consider three special orbits: purely vertical orbits, purely radial orbits and circular orbits in the plane, for which the actions are well defined. For example, the flattened logarithmic potential is

\begin{equation} \Phi (R,z) = \frac{1}{2}v_{0}^2\log \left ((R/r_{0})^2 + (z/qr_{0})^2\right ). \end{equation}

(38)

Using the same interpolation techniques, we reach the following expression for an approximate Hamiltonian:

\begin{equation} \mathcal {H} = v_{0}^2\log \bigg (\frac{\sqrt{2\pi }(J_{R}+J_{z}/q)+\sqrt{e}\vert L_{z}\vert}{v_{0}r_{0}}\bigg ). \end{equation}

(39)

Given the fact that the actions cannot be written down as a quadratures for general orbits in this potential, this expression would have to be tested using an approximate method, such as described in Binney (2012). We look to do exactly this in future investigations.

We thank Donald Lynden-Bell and Simon Gibbons for some extremely useful conversations and insights. AW and AB acknowledge the support of STFC.

REFERENCES

Abraham

R.

Marsden

J.

,

Foundations of Mechanics

,

1978

London

Addison-Wesley

OpenURL Placeholder Text

Agnello

A.

Evans

N. W.

Romanowsky

A. J.

Brodie

J.-P.

,

MNRAS

,

2014

preprint (arXiv:1401.4461)

Arnold

V.

,

Mathematical Methods of Classical Mechanics

,

1989

Berlin

Springer-Verlag

OpenURL Placeholder Text

Binney

J.

,

MNRAS

,

2012

, vol.

426

pg.

1324

Binney

J. J.

Evans

N. W.

,

MNRAS

,

2001

, vol.

327

pg.

L27

Binney

J.

Tremaine

S.

,

Galactic Dynamics

,

2008

2nd edn

Princeton, NJ

Princeton Univ. Press

OpenURL Placeholder Text

Born

M.

,

The Mechanics of the Atom

,

1927

London

G. Bell and Sons

OpenURL Placeholder Text

Bovy

J.

,

ApJ

,

2014

preprint (arXiv:1401.2985)

El-Zant

A.

Shlosman

I.

Hoffman

Y.

,

ApJ

,

2001

, vol.

560

pg.

636

Evans

N. W.

,

MNRAS

,

1994

, vol.

267

pg.

333

Evans

N. W.

Williams

A.

,

MNRAS

,

2014

submitted

Evans

N. W.

de Zeeuw

P. T.

Lynden-Bell

D.

,

MNRAS

,

1990

, vol.

244

pg.

111

Goldstein

H.

,

Classical Mechanics

,

1980

2nd edn

London

Addison-Wesley

OpenURL Placeholder Text

Gondolo

P.

Silk

J.

,

Phys. Rev. Lett.

,

1999

, vol.

83

pg.

1719

Goodman

J.

Binney

J.

,

MNRAS

,

1984

, vol.

207

pg.

511

Grillmair

C. J.

Dionatos

O.

,

ApJ

,

2006

, vol.

643

pg.

L17

Katz

H.

McGaugh

S. S.

Sellwood

J. A.

de Blok

W. J. G.

,

MNRAS

,

2014

, vol.

439

pg.

1897

Koposov

S. E.

Rix

H.-W.

Hogg

D. W.

,

ApJ

,

2010

, vol.

712

pg.

260

Lynden-Bell

D.

,

MNRAS

,

2010

, vol.

402

pg.

1937

Martizzi

D.

Teyssier

R.

Moore

B.

,

MNRAS

,

2013

, vol.

432

pg.

1947

Navarro

J. F.

Frenk

C. S.

White

S. D. M.

,

ApJ

,

1996

, vol.

462

pg.

563

Read

J. I.

Gilmore

G.

,

MNRAS

,

2005

, vol.

356

pg.

107

Richstone

D. O.

,

ApJ

,

1982

, vol.

252

pg.

496

Sanders

J. L.

Binney

J.

,

MNRAS

,

2013

, vol.

433

pg.

1813

Sanders

J. L.

Binney

J.

,

MNRAS

,

2014

preprint (arXiv:1401.3600)

© 2014 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society

Download all slides

Views

759

Altmetric

Total Views 759

371 Pageviews

388 PDF Downloads

Since 1/1/2017

Month:	Total Views:
January 2017	1
February 2017	3
March 2017	4
August 2017	1
October 2017	6
November 2017	5
December 2017	5
January 2018	2
February 2018	5
March 2018	22
April 2018	18
May 2018	17
June 2018	17
July 2018	17
August 2018	10
September 2018	19
October 2018	14
November 2018	8
December 2018	12
January 2019	16
February 2019	11
March 2019	34
April 2019	21
May 2019	15
June 2019	12
July 2019	19
August 2019	11
September 2019	7
October 2019	9
November 2019	7
December 2019	10
January 2020	10
February 2020	3
March 2020	13
April 2020	4
May 2020	4
June 2020	13
July 2020	5
August 2020	11
September 2020	14
October 2020	1
November 2020	10
December 2020	10
January 2021	6
February 2021	6
March 2021	6
April 2021	9
May 2021	7
June 2021	4
July 2021	3
September 2021	4
October 2021	7
November 2021	3
December 2021	1
January 2022	9
February 2022	12
April 2022	3
May 2022	1
June 2022	10
July 2022	6
August 2022	3
September 2022	12
October 2022	8
November 2022	1
December 2022	2
January 2023	2
February 2023	6
March 2023	2
April 2023	2
May 2023	2
June 2023	6
July 2023	7
August 2023	9
September 2023	4
October 2023	15
November 2023	16
December 2023	12
January 2024	4
February 2024	5
March 2024	4
April 2024	8
May 2024	13
June 2024	8
July 2024	4
August 2024	6
September 2024	6
October 2024	9
November 2024	7
December 2024	5
January 2025	7
February 2025	1
March 2025	8
April 2025	2