# Galaxy halo expansions: a new biorthogonal family of potential-density pairs

Galaxy halo expansions: a new biorthogonal family of potential-density pairs Abstract Efficient expansions of the gravitational field of (dark) haloes have two main uses in the modelling of galaxies: first, they provide a compact representation of numerically constructed (or real) cosmological haloes, incorporating the effects of triaxiality, lopsidedness or other distortion. Secondly, they provide the basis functions for self-consistent field expansion algorithms used in the evolution of N-body systems. We present a new family of biorthogonal potential-density pairs constructed using the Hankel transform of the Laguerre polynomials. The lowest order density basis functions are double-power-law profiles cusped like ρ ∼ r−2+1/α at small radii with asymptotic density fall-off like ρ ∼ r−3−1/(2α). Here, α is a parameter satisfying α ≥ 1/2. The family therefore spans the range of inner density cusps found in numerical simulations, but has much shallower – and hence more realistic – outer slopes than the corresponding members of the only previously known family deduced by Zhao and exemplified by Hernquist & Ostriker. When α = 1, the lowest order density profile has an inner density cusp of ρ ∼ r−1 and an outer density slope of ρ ∼ r−3.5, similar to the famous Navarro, Frenk & White (NFW) model. For this reason, we demonstrate that our new expansion provides a more accurate representation of flattened NFW haloes than the competing Hernquist-Ostriker expansion. We utilize our new expansion by analysing a suite of numerically constructed haloes and providing the distributions of the expansion coefficients. methods: numerical, galaxies: kinematics and dynamics 1 INTRODUCTION A galaxy is often modelled as a sum of simple analytical components, such as a Hernquist (1990) bulge or a Navarro, Frenk & White (1997, NFW) or axisymmetric logarithmic halo (Binney & Tremaine 1987; Evans 1993). Although crude, such an approach produces manageable zeroth-order models from which further insights can be garnered. However, in the era of powerful numerical simulations of galaxy formation as well as accurate positional and kinematic data in our own Milky Way, such broad-brush models fail to capture the complexity of structures observed in both simulations and data. A richer and more flexible way of representing the gravitational field of galaxies (and their components) is to use basis-function or halo-expansion techniques (see e.g. Clutton-Brock 1972; Hernquist & Ostriker 1992; Zhao 1996; Lowing et al. 2011). In this framework, small higher order deviations away from the zeroth-order galaxy model are efficiently captured by terms in an orthogonal series (analogous to a Fourier series), truncated at some order nmax. Separate sets of orthogonal basis functions are used for the density and potential, so the method is sometimes known as the biorthogonal series expansion method. The two series are related through the Poisson equation, enabling the same set of coefficients to be used for both. This represents a significant improvement in flexibility, allowing axisymmetric, triaxial, lopsided, or distorted density distributions to be built up from an underlying simple model. The original motivation for the introduction of biorthogonal basis sets was the study of the stability of collisionless stellar systems. When decomposed into basis functions, the normal modes of stellar systems correspond to the eigenvectors of matrix equations (e.g. Polyachenko & Shukhman 1981; Fridman & Polyachenko 1984; Weinberg 1989; Palmer 1994; Evans & Read 1998). Stability studies of spherical galaxies often use spherical Bessel functions as the radial basis functions and spherical harmonics as the angular basis functions. These are natural choices as they are the eigenfunctions of the Laplacian in spherical coordinates. However, spherical Bessel functions form a discrete basis only if the radial range is finite, so that the galaxy model has to be truncated or of finite extent. For a semi-infinite radial range, Saha (1991) used a transformation to map the radial range to a finite one, and then used the brute force Gram–Schmidt method to build orthogonal basis functions. Subsequently, basis-function techniques were used to provide algorithms to evolve collisionless stellar systems, partly as a check on the results of linear stability theory. For example, Allen, Palmer & Papaloizou (1990) used a biorthogonal spherical Bessel basis-function expansion to study the radial and circular orbit instability in spherical galaxy models. As instabilities in stellar systems often arise from nearly resonant orbits, so accurate modelling of the precession of individual orbits for several orbital periods is very important. Hernquist & Ostriker (1992) clearly articulated the advantages of basis-function expansions for which the zeroth-order model actually resembles a galaxy. They used the clever method pioneered by Clutton-Brock (1972) to transform the (radial part of the) Laplacian into an equation whose eigenfunctions are already known. They tellingly remarked that ‘for reasons not immediately obvious to us, Clutton-Brock's approach has been virtually ignored in the literature’. In the Hernquist & Ostriker (1992) basis-function expansion, the zeroth-order model is a Hernquist (1990) density profile, which resembles an elliptical galaxy or dark matter halo with ρ ∼ 1/r at small radii and ρ ∼ 1/r4 at large radii. Direct simulation of N-body gravitational systems requires O(N2) operations to calculate the forces on the particles. However, stellar systems are typically collisionless, meaning two-body relaxation effects between bodies are small and the particles chiefly move under the influence of the smooth, mean gravitational field. It is this fact that enables the use of efficient basis-function expansion techniques for calculating the force (see e.g. Dehnen 2014). When forces are calculated from the truncated series, we automatically avoid spurious relaxation effects that may otherwise have been introduced due to undersampling of the number of particles in the system. The coefficients are only calculated once (per time-step), and the series of functions is evaluated once per particle to find the force, so such methods are O(N) in the number of particles. More recently, basis-function techniques have enjoyed a resurgence, inspired by a very different motivation. Lowing et al. (2011) pointed out the value of describing the evolution of direct and computationally expensive simulations of galaxy formation using basis functions. Snapshots of the original simulation with redshift are decomposed with basis functions. The entire simulation can then be replayed many times at will. For example, new objects can be inserted into the simulations and their behaviour studied as if they had been present originally. This technique has already been exploited in studies of the evolution of tidal streams in the Galaxy (Ngan et al. 2015, 2016). However, the gamut of applications is far broader, including the evolution of accreting subhaloes, the dynamics of satellite galaxies and globular clusters in dark haloes, and the response of galactic discs to the changing dark-halo potential. Despite its significant advantages and diverse applications, an important limitation of the basis-function technique is the relative paucity of available sets of ‘analytical’ basis sets – i.e. expressible in a closed form, and preferably obeying some form of recurrence relation so that computation of the basis functions is efficient. Any such basis sets are useful in their own right, but they are particularly interesting if the zeroth-order potential-density pair is close to that of galaxies, allowing for efficient representation using only a few terms. In fact, both Lowing et al. (2011) and Ngan et al. (2015) decomposed numerical halo simulations into the Hernquist & Ostriker (1992) basis-function expansion, as the alternatives are very limited. The purpose of this paper is to broaden the application of the basis-function technique by introducing new sets of biorthogonal basis functions and discussing their applications. We begin by reviewing the presently available expansions and presenting a general approach to the construction of biorthogonal expansions in Section 2 before presenting new specific cases in Section 3. Readers eager to apply the new basis expansion can find the appropriate expressions at the beginning of Section 4. We continue by expanding a series of N-body cosmological haloes that are imitations of the Milky Way dark halo using our new expansion method. We sum up and point out further avenues of development, both on the theoretical and applicable side, in Section 5. 2 THEORY OF BIORTHOGONAL EXPANSIONS 2.1 History The story really begins with the two pioneering papers by Clutton-Brock. Seeking solutions to the cylindrical Poisson equation, Clutton-Brock (1972) introduced the Bessel functions as a set of biorthogonal functions with a continuous eigenvalue on the interval (0, ∞). From these functions, a biorthogonal set with discrete index n was constructed by integrating a particular set of orthogonal functions gn(k) against this continuous eigenvalue (i.e. performing a Hankel transform). This procedure produced a basis set suitable for infinitesimally thin disc galaxies. For Clutton-Brock's choice of gn(k), the lowest order potential-density pair is that of the Kuzmin (1956) or Toomre (1963) disc, though he did not explicitly evaluate the integral and instead left the result as a recurrence relation. Subsequently, Aoki & Iye (1978) reproduced the results of Clutton-Brock (1972), but they explicitly evaluated the Hankel transform to give results written in a closed form. For spherical systems, Clutton-Brock (1973) took a different approach, constructing a biorthogonal basis set for spherical galaxies directly out of the Gegenbauer polynomials together with a particular change of variable. The zeroth-order density profile reproduced the Plummer (1911) model, but this idea was then extended to cover the Hernquist (1990) model by Hernquist & Ostriker (1992) before finally being generalized by Zhao (1996) to a family of basis sets for the hypervirial models (Evans & An 2005). The extension of the Clutton-Brock (1972) approach has only received limited attention in the literature, but we shall demonstrate that it provides a powerful alternative route into constructing further basis sets. Polyachenko & Shukhman (1981) considered the Bessel functions over a finite interval (0, a), where a is the radial extent of the galaxy. This gives a biorthogonal basis set with discrete indices, although the expansion suffers from the issue that the Bessel functions by themselves do not resemble any particular well-known galactic profile. Unlike Clutton-Brock, they did not consider the case of a continuous eigenvalue on the semi-infinite interval (0, ∞). They did however identify the additional degree of freedom, which allows for an adjustment in the asymptotic behaviour of the density basis functions and which we also exploit later. More recently, Rahmati & Jalali (2009) explicitly extended the derivation of Clutton-Brock (1972) and Aoki & Iye (1978) to the spherical case, making an analogous choice of gn(k), and hence deriving a basis set whose lowest order density is the perfect sphere (de Zeeuw 1985). Our work is a generalization of this result, analogous to the way in which Zhao (1996) is a generalization of Clutton-Brock (1973). In fact, Rahmati & Jalali (2009)'s basis set corresponds to the lowest member of our new family of biorthogonal pairs. 2.2 Orthogonal bessel solutions to the Poisson equation We seek solutions in series to the Poisson equation ∇2Φ = 4πGρ (where we work with G = 1 throughout). We represent the potential and density of a self-gravitating system as a sum over basis functions   $$\Phi ({\boldsymbol r}) = \sum _{nlm} \, C_{nlm} \, \Phi _{nlm}({\boldsymbol r}) \qquad \qquad \rho ({\boldsymbol r}) = \sum _{nlm} \, C_{nlm} \, \rho _{nlm}({\boldsymbol r}),$$ (1)where the basis functions are labelled by the tuple (n, l, m). The condition of biorthogonality is written as   $$-\int {\rm d}^{\,3} {\boldsymbol r} \, \Phi _{nlm} \, \rho ^*_{n^\prime l^\prime m^\prime } = N_{nlm} \, \delta _{nlm}^{n^\prime l^\prime m^\prime }, \qquad \qquad \nabla ^2\Phi _{nlm} = 4\pi \rho _{nlm}.$$ (2)We factor out the spherical harmonics Ylm(θ, ϕ) as   $$\Phi _{nlm}({\boldsymbol r}) \equiv \Phi _{nl}(r) Y_{lm}(\theta , \phi ),\qquad \qquad \rho _{nlm}({\boldsymbol r}) \equiv \rho _{nl}(r) Y_{lm}(\theta , \phi ),$$ (3)leaving the equation satisfied by the radial basis functions as   $$\left(\nabla ^2 - \frac{l(l+1)}{r^2}\right)\Phi _{nl} = 4\pi \rho _{nl}.$$ (4) On a finite domain with spherical symmetry, it is natural to choose as a basis the eigenfunctions of the Laplacian, $$f_{klm}({\boldsymbol r}) = j_l(kr)\,Y_{lm}(\theta , \phi )$$, where jl(kr) = r−1/2Jl + 1/2(kr) is a spherical Bessel function. The Bessel function of order μ is defined by the equation   \begin{eqnarray} -\frac{{\rm d}}{{\rm d}z} \left( z^{2\mu +1} \frac{{\rm d}}{{\rm d}z} \left( y_k(z) \right) \right) = k^2 z^{2\mu +1} y_k(z), \ \ \ \ y_k(z) \equiv z^{-\mu } J_\mu (kz), \end{eqnarray} (5)and the eigenfunction fklm obeys the equation ∇2fklm = k2 fklm. However, on an infinite domain, we cannot use these eigenfunctions directly for a series expansion as the radial eigenvalue k becomes continuous: every element of (0, ∞) is a valid eigenvalue k. We therefore convert the basis to a discrete one; first choosing a set of functions gn(k) that are orthogonal on the interval (0, ∞) with respect to the measure k dk and then integrating against fklm (essentially performing a Hankel transform). This produces a biorthogonal basis set, which is suitable for expanding both the potential Φ and the density ρ. Given that we are using different sets of functions to represent the potential and the density, we can exploit an additional freedom first noticed by Polyachenko & Shukhman (1981) and introduce the parameter α in the r-dependence of the Bessel functions. Accordingly, we define the potential and density basis functions as   \begin{eqnarray} \Phi _{nl}(r) \equiv \frac{-A_{nlm}}{\sqrt{r}} \int _0^\infty {\rm d}k \, g_n(k) \, J_{\mu } (kr^{1/(2\alpha )}), \qquad \rho _{nl}(r) \equiv \frac{A_{nlm}}{16\pi \alpha ^2} r^{1/\alpha - 5/2} \int _0^\infty {\rm d}k \, k^2 \, g_n(k) \, J_{\mu } (kr^{1/(2\alpha )}), \end{eqnarray} (6)where μ ≡ α (1 + 2l), and Anlm are constants that will be adjusted to give the desired overall normalization. These functions are biorthogonal in the sense of (2) if the auxiliary functions gn(k) are suitably chosen. Specifically, they must satisfy an orthogonality condition which can be derived thus,   \begin{eqnarray} -\int _0^\infty {\rm d}r \, r^2 \, \Phi _{nl} \, \rho _{n^\prime l} &=& \frac{A_{nlm}\,A_{n^\prime lm}}{8\pi \alpha } \, \int _0^\infty {\rm d}k \, g_n(k) \int _0^\infty {\rm d}q \, q^2 g_{n^\prime }(q) \int _0^\infty {\rm d}z \, z \, J_{\mu }(k z) J_{\mu }(q z) \qquad (z \equiv r^{1/(2\alpha )}) \nonumber\\ &=& \frac{A_{nlm}\,A_{n^\prime lm}}{8\pi \alpha } \, \int _0^\infty {\rm d}k \, k \, g_n(k) \, g_{n^\prime }(k) = \frac{A_{nlm}^2}{8\pi \alpha } \, I_n \, \delta _{n n^\prime }. \end{eqnarray} (7)Here, In is the normalization constant for the functions gn(k) are orthogonal with respect to the weight function k. Let us also write Jlm as the normalization constant for the spherical harmonics Ylm. The overall normalization constant is therefore   $$N_{nlm} = \frac{I_n \, J_{lm} \, A_{nlm}^2}{8\pi \alpha },$$ (8)so that the functions are biorthonormal (i.e. Nnlm = 1) if we set   $$A_{nlm} = \sqrt{\frac{8\pi \alpha }{I_n\,J_{lm}}},$$ (9)in equation (6). The coefficients Cnlm can be calculated by the following integrals over ρ and Φ. When the density is formed from a cloud of point particles as in a numerical simulation, the integral reduces to a sum, which (given this choice of normalization) is written as   \begin{eqnarray} \rho ({\boldsymbol r}) = \sum _i\,m_i \delta ^3({\boldsymbol r} - {\boldsymbol r}_{i}), \qquad C_{nlm} = -\int {\rm d}^{\,3}{\boldsymbol r} \, \Phi _{nlm}({\boldsymbol r})\,\rho ({\boldsymbol r}) = \sum _i m_i\,\Phi _{nlm}({\boldsymbol r}_{i}). \end{eqnarray} (10)One of the advantages of biorthogonal basis sets is that the coefficients in both the density and the potential are the same and can be evaluated very easily from an N-body realization. Note that the use of a discrete sum instead of a continuous integration implies that there is always some numerical (Poisson) noise introduced by discreteness effects. 2.3 Physical boundary conditions The approach presented above yields general sets of biorthogonal functions, but there are further considerations when wishing to construct basis sets for physical systems. Consider the following argument about the required asymptotic behaviour of the radial basis functions. For a given radial order n and angular order l, we define the mass enclosed by a shell with inner radius a and outer radius b as   \begin{eqnarray} M_{nl}(a, b) = 4\pi \int _a^b {\rm d}r \, r^2 \, \rho _{nl} = \int _a^b {\rm d}r \, \left[\frac{{\rm d}}{{\rm d}r}\left(r^2\frac{{\rm d}}{{\rm d}r}\Phi _{nl}\right) - l(l+1)\Phi _{nl}\right]. \end{eqnarray} (11)In order for the density and potential to describe a physical system, we impose the constraint that the mass enclosed by a thin shell goes to zero at the origin and at infinity, and that the potential remains finite everywhere. That is, we require that   \begin{eqnarray} \lim _{a\rightarrow 0} M_{nl}(0,a) = 0, \qquad \lim _{a\rightarrow \infty } M_{nl}(a, \infty ) = 0, \qquad \Phi _{nl}(r) < \infty . \end{eqnarray} (12)These translate into equations that must be satisfied by the derivative of Mnl, which is proportional to ρnl, and so examining (11) we see that the asymptotic solutions for the potential must satisfy: Φnl ∼ rl as r → 0; and Φnl ∼ r−1−l as r → ∞. 3 CONSTRUCTION OF THE NEW BASIS SET 3.1 Zeroth-order models With our general framework in place, we now turn to the job of producing expansions whose zeroth-order model resembles realistic galactic components. This is a lengthy business, so we begin by describing some of the zeroth-order models in our sequence of basis-function expansions. Of course, the zeroth-order models are all spherically symmetric, but the higher order terms in the basis-function expansion describe the effects of flattening or lopsidedness or distortion. Hence, the method can build very general density distributions around the zeroth-order model. The models are labelled by a parameter α; however, those with 0 < α < 1/2 have vanishing density at the origin, and so we consider them no further. Models of interest in galactic dynamics have α ≥ 1/2. The behaviour of the density at small and large radii is   \begin{eqnarray} \lim _{r\rightarrow 0} \rho _{000} (r) \sim r^{-2 + 1/\alpha }, \qquad \lim _{r\rightarrow \infty } \rho _{000} (r) \ \sim r^{-3 -1/(2\alpha )}, \end{eqnarray} (13)so that the α = 1/2 model is cored, whilst the remainder are cusped. The density of all the models falls off with a logarithmic slope between −3 and −4. We proceed to list some of the zeroth-order models that are obtained when α is integer or half-integer, as in these cases the potential reduces to elementary functions; in general, real values of α are permitted, but the potential must be evaluated using the incomplete beta function [see equation (18)]. When α = 1/2, the zeroth-order model is the perfect sphere of de Zeeuw (1985), which has density and potential   \begin{eqnarray} \rho _{000} = \frac{1}{\sqrt{2} \, \pi ^{3/2}} \, \frac{1}{\left(1 + r^2\right)^{2}}, \qquad \Phi _{000} = -\sqrt{2\over \pi } \frac{\arctan {r}}{r}. \end{eqnarray} (14)This is the spherical limit of the perfect ellipsoid, which provides triaxial densities that are close to the luminosity profiles of elliptical galaxies. There is a large literature on these models as the potential is separable or Stäckel, and so the orbits and action-angles can be found as quadratures (e.g. Binney & Tremaine 1987; Bertin 2014). Notice that the density has a harmonic core at the centre but falls like ρ ∼ r−4 at large radii. When α = 1, we obtain the super-NFW model (e.g. Evans & An 2006; An & Zhao 2013; Lilley, Evans & Sanders 2017). This has a central density cusp with ρ ∝ r−1 similar to the Hernquist (1990) or NFW models. Asymptotically, the density falls like ρ ∝ r−3.5 at large radii, which is somewhat faster than the NFW model, but has the distinct advantage of finite total mass. There is evidence to suggest that the NFW may be too shallow at large radii based on the recent work on the splashback radius by Diemer et al. (2017). The potential-density pair of the super-NFW model is   \begin{eqnarray} \rho _{000} = \frac{3}{16\pi } \, \frac{1}{r(1 + r)^{5/2}}, \qquad \Phi _{000} = \frac{-1}{1 + r + \sqrt{1 + r}}. \end{eqnarray} (15)This is ideal for modelling the profiles of dark haloes found in cosmological simulations. This potential-density pair may also be used as a simple model of a galaxy or star cluster. Lilley et al. (2017) show that isotropic and anisotropic distribution functions and solutions of the Jeans equations can be obtained, whilst the projected density is close to de Vaucouleurs and Sérsic profiles. There are of course other finite-mass models which have such desirable properties. A number of models with ρ ∝ r−1 at the centre and with an asymptotic fall off with logarithmic gradient between −3 and −4 are already available in the literature (e.g. An & Evans 2006; An & Zhao 2013). In general, these models are not the zeroth-order terms of any biorthorgonal expansion, so their axisymmetric and triaxial siblings are much more cumbersome to construct. More strongly cusped models can also be obtained as zeroth-order models. When α = 3/2, we have   \begin{eqnarray} \rho _{000} = \frac{2\sqrt{2}}{9 \pi ^{3/2}} \, \frac{1}{r^{4/3}\left(1 + r^{2/3}\right)^{3}}, \qquad \Phi _{000} = \sqrt{2\over \pi }\Bigl [ {1\over r^{2/3}(1 + r^{2/3})} - {{\arctan } (r^{1/3})\over r}\Bigr ]. \end{eqnarray} (16)The central density cusp is now ρ ∼ r−4/3, similar to the inner regions of the massive cosmological clusters studied by Diemand et al. (2005). When α = 2, we have   \begin{eqnarray} \rho _{000} = \frac{15}{64 \pi } \, \frac{1}{r^{3/2}(1 + r^{1/2})^{7/2}}, \qquad \Phi _{000} = -{1 + 2\sqrt{1+\sqrt{r}} \over (1 + \sqrt{1+\sqrt{r}})^2 (1+\sqrt{r})^{3/2}}. \end{eqnarray} (17)This has an inner density profile ρ ∼ r−3/2, similar to the very steepest cusps found in cosmological simulations (Moore et al. 1998). Explicitly, the lowest order basis functions are   \begin{eqnarray} \rho _{000} = {2^{\alpha +1}\Gamma (3/2+\alpha )\over \sqrt{\pi }} {1 \over r^{2 - 1/\alpha }(1+r^{1/\alpha })^{\alpha + 3/2}}\qquad \Phi _{000} = -{2^{\alpha -1} \Gamma (\alpha + 3/2) \over \sqrt{\pi }} {\mathcal {B}_\chi (\alpha ,1/2) \over r}, \qquad \chi \equiv \frac{r^{1/\alpha }}{1 + r^{1/\alpha }}, \end{eqnarray} (18)where $$\mathcal {B}_x(a,b)$$ is the incomplete beta function. The full suite of models has a range of inner density profiles suitable for representing dark haloes with cores and weak or strong cusps, and α can be tuned to match the behaviour of a given halo (see Section 4.3). 3.2 Choice of auxiliary function Now that we have given the spherically symmetric lowest order models, we show how to construct the associated higher order basis functions. We introduce a two-parameter family of auxiliary functions gn(k). We derive an expression for the potential basis functions but find that the requirement of physicality reduces this set to a single-parameter family. This reduction of complexity allows us to express the potential and density basis functions in a more succinct fashion in Sections 3.3 and 3.4. Inspired by Rahmati & Jalali (2009),1 and initially following their calculations, we now choose   \begin{eqnarray} g_n(k) \equiv L^{(\eta )}_n(2k) \, k^{(\eta -1)/2} \, \exp (-k), \end{eqnarray} (19)where $$L^{(\eta )}_n(x)$$ are the associated Laguerre polynomials (Olver et al. 2016, Section 18.3), and η is (for now) a free parameter. The Laguerre polynomials are orthogonal on (0, ∞) with respect to a weight function ω(x) ≡ xη exp (−x),   \begin{eqnarray} \int _0^\infty {\rm d}x \, L^{(\eta )}_n(x) \, L^{(\eta )}_{n^\prime }(x) \omega (x) = \delta _{n n^\prime } {\Gamma (n+\eta +1)\over \Gamma (n+1)}. \end{eqnarray} (20)Our gn(k) consist of a Laguerre polynomial multiplied by a factor of $$\sqrt{\omega (2k)/k}$$, hence ensuring the orthogonality on (0, ∞) with respect to k dk. We can obtain an expression for Φnl using the integral equation 6.621(1) of Gradshteyn & Ryzhik (2014),   \begin{eqnarray} \int _0^\infty {\rm d}u \, \exp (-us) u^\nu J_\mu (u) = \Gamma (\nu + \mu + 1) \, (1 + s^2)^{-(1 + \nu )/2} \, P^{(-\mu )}_\nu \left(\frac{s}{\sqrt{1 + s^2}}\right). \end{eqnarray} (21)Setting u = kr1/(2α) and s = r−1/(2α), and using the explicit polynomial representation of the Laguerre polynomials   $$L^{(\eta )}_n(x) = \sum _{j=0}^n (-1)^j \frac{1}{j!} \left({{n+\eta }\atop n-j}\right) x^j,$$ (22)we find   $$\Phi _{nl} = \sum _{j=0}^n A_{n j\mu \eta } \left[ r \left(1 + r^{1/\alpha }\right)^{(\eta +1)/2 + j}\right]^{-1/2} \, P_{\frac{\eta -1}{2} + j}^{\left(-\mu \right)}\left(\frac{1}{\sqrt{1 + r^{1/\alpha }}}\right), \qquad A_{n j\mu \eta } \equiv \frac{(-2)^j}{j!} \left({{n+\eta }\atop n-j}\right) \Gamma \left(j + \frac{\eta +1}{2} + \mu \right).$$ (23)We use the results of Section 2.3 on asymptotic limits to fix the parameter η. The asymptotic behaviour of the associated Legendre function (Olver et al. 2016, equation 14.8(i)) as z → 1 from above is $$P^{(\mu )}_\nu (z) \sim (1 - z)^{-\mu /2}$$. Therefore, as r → 0, every term in Φnl goes as   $$\Phi _{nl} \sim r^{-1/2} \, \left(1 - \frac{1}{\sqrt{1 + r^{1/\alpha }}}\right)^{\mu /2} \propto r^{\mu /(2\alpha ) - 1/2} = r^l,$$ (24)so the first limit implied by (12) is already satisfied. On the other hand, as r → ∞, $$1/\sqrt{1 + r^{1/\alpha }} \sim r^{-1/(2\alpha )}$$, so the associated Legendre function goes to a constant (Olver et al. 2016, equation 14.5.1) and we are left with the prefactor   $$\Phi _{nl} \sim \left[ r \left(r^{1/\alpha }\right)^{(\eta +1)/2 + j}\right]^{-1/2} \rightarrow r^{-1/2 - (\eta + 1)/(4\alpha )},$$ (25)where we have used the fact that the j = 0 term dominates the sum as r → ∞. By (12), we require this limit to be r−l−1, which we can achieve by setting η = 2μ − 1, so the gn-normalization constant is   $$I_n = \frac{\Gamma (n+2\mu )}{2^{2\mu }\,n!}.$$ (26)Now that η is no longer a free parameter, we can obtain the basis functions in their most convenient form. 3.3 Recurrence relation for the potential basis functions The sum (23) is unsatisfactory for several reasons: (i) it is numerically unstable for high n, (ii) n and j are coupled in such a way that calculating n basis functions requires O(n2) operations, and (iii) associated Legendre functions of non-integer or negative degree or order are rarely implemented numerically. Therefore, departing from Rahmati & Jalali (2009), we seek a superior expression, which can be obtained using the recurrence relation for the Laguerre polynomials,   $$n \, L_n^{(\alpha )}(x) = (n + \alpha ) \, L_{n-1}^{(\alpha )}(x) - x \, L_{n-1}^{(\alpha +1)}(x).$$ (27)Using equation (19), we can immediately write down   $$n \, g_n(k) = (n + 2\mu - 1) \, g_{n-1}(k) - 2 \exp (-k) \, k^\mu \, L_{n-1}^{(2\mu )}(2k).$$ (28)So, we obtain the following recurrence relation for the basis functions   $$n \, \Phi _{nl} = (n + 2\mu - 1) \, \Phi _{n-1,l} + {2A_{nlm}\over \sqrt{r}} \int _0^\infty k^\mu \, L_{n-1}^{(2\mu )}(2k) \, \exp (-k) \, J_\mu (kz) \, {\rm d}k, \qquad (z \equiv r^{1/(2\alpha )}).$$ (29)Evaluating the latter integral requires some work. First, we note that a generating function for the Laguerre polynomials is   \begin{eqnarray} \sum _{n=0}^\infty \, t^n \, L_n^{(\lambda )}(k) = \frac{\exp \left(-tk/(1-t) \right)}{\left(1-t\right)^{\lambda +1}}, \end{eqnarray} (30)so that, using the following Hankel transform (Gradshteyn & Ryzhik 2014, equation 6.623(1))   \begin{eqnarray} \int _0^\infty \, x^\nu \, \mathrm{e}^{-ax} \, J_\nu (xy) \, {\rm d}x = \frac{2^\nu \, \Gamma (\nu + 1/2)}{\sqrt{\pi }} \, \frac{y^\nu }{\left(a^2 + y^2\right)^{\nu + 1/2}}, \end{eqnarray} (31)as well as the generating function for the Gegenbauer polynomials $$C_n^{(\lambda )}(\xi )$$  \begin{eqnarray} \sum _{n=0}^\infty \, t^n \, C_n^{(\lambda )}(\xi ) = \frac{1}{\left(1 - 2\xi t + t^2\right)^\lambda } = \frac{\left(1 + z^2\right)^\lambda }{\left[(1+t)^2 + (1-t)^2z^2\right]^\lambda }, \qquad \left(\xi \equiv \frac{z^2 - 1}{z^2 + 1} \equiv \frac{r^{1/\alpha } - 1}{r^{1/\alpha } + 1}\right), \end{eqnarray} (32)we can find the Hankel transform of the remaining term in the recurrence relation,   \begin{eqnarray} \sum _{n=0}^\infty \, t^n \, \int _0^\infty k^\mu \, L_n^{(2\mu )}(2k) \, \exp (-k) \, J_\mu (kz) \, {\rm d}k & =& \int _0^\infty k^\mu \frac{\exp (-k(1+t)/(1-t))}{\left(1-t\right)^{2\mu + 1}}J_\mu (kz)\,{\rm d}k = \frac{2^\mu \, \Gamma (\mu + 1/2)}{\sqrt{\pi }} \, \frac{z^\mu }{\left[(1+t)^2 + (1-t)^2 z^2\right]^{\mu + 1/2}} \nonumber \\ &=& \sum _{n=0}^\infty t^n \, \frac{2^\mu \, \Gamma (\mu + 1/2)}{\sqrt{\pi }} \, \frac{z^\mu }{\left(1+z^2\right)^{\mu + 1/2}} \, C_n^{(\mu + 1/2)}(\xi ). \end{eqnarray} (33)So the recurrence relation becomes   \begin{eqnarray} n \, \Phi _{nl} = (n + 2\mu - 1) \, \Phi _{n-1,l} + A_{nlm} \frac{2^{\mu +1} \, \Gamma (\mu +1/2)}{\sqrt{\pi }} \, \frac{r^l \, C_{n-1}^{(\mu +1/2)}\left(\xi \right)}{\left(1 + r^{1/\alpha }\right)^{\mu +1/2}}. \end{eqnarray} (34)This means that, in theory, the (n + 1)-th basis function can be trivially calculated from the n-th by simply adding on a single extra term. Evaluating n basis functions therefore requires O(n) steps (although see Section 4.1). A similar approach is taken in Aoki & Iye (1978) in the disc setting. So far, we have only stated without proof the (n = 0, l = 0) case. We now obtain an explicit expression for all the potential basis functions, from which the zeroth order (n = 0, l ≥ 0) can be extracted. We begin by writing our auxiliary function as   \begin{eqnarray} g_n(k) = k^{\mu - 1} \, \exp (-k) \, L_n^{(2\mu - 1)}(2k) = k^{\mu - 1} \, \exp (-k) \left( L_n^{(2\mu - 1)}(0) - 2k\sum _{j=0}^{n-1} \frac{\left({{n+2\mu - 1}\atop n-1-j}\right)}{(n-j)\left({{n}\atop j}\right)} L_j^{(2\mu )}(2k)\right), \end{eqnarray} (35)where again μ = α (1 + 2l). Aside from the constant term, the Hankel transforms can be done using equation (33). To evaluate the Hankel transform of the constant term, we use equation (21) to find [denoting by $$\mathcal {B}_x(a,b)$$ the incomplete beta function]   \begin{eqnarray} \int _0^\infty \, k^{\mu -1} \, \exp (-k) \, J_\mu (kz) \, {\rm d}k = \frac{\Gamma (2\mu )}{\left(1+z^2\right)^{\mu /2}} \, P^{(-\mu )}_{\mu -1}\left(\frac{1}{\sqrt{1+z^2}}\right) = \frac{\Gamma (2\mu )}{\Gamma (\mu )\,2^\mu \,z^\mu } \, \mathcal {B}_\chi (\mu ,1/2), \qquad \left(\chi \equiv \frac{r^{1/\alpha }}{1 + r^{1/\alpha }}\right). \end{eqnarray} (36)So, reassembling the terms, we obtain an expression for the potential as   \begin{eqnarray} \Phi _{nl}(r) = -A_{nlm} \, \frac{2^\mu \, \Gamma (\mu + 1/2) \, \Gamma (n + 2\mu )}{\sqrt{\pi } \, n!} \left[\frac{\mathcal {B}_\chi (\mu ,1/2)}{2\,\Gamma (2\mu )\,r^{1+l}} - \frac{2r^l}{\left(1 + r^{1/\alpha }\right)^{\mu +1/2}}\sum _{j=0}^{n-1}\frac{j! \, C_j^{(\mu +1/2)}(\xi )}{\Gamma (j + 2\mu + 1)}\right], \qquad \left( \chi \equiv \frac{1+\xi }{2}\right), \end{eqnarray} (37)with μ ≡ α (1 + 2l). The only special functions required to evaluate the potential are the incomplete beta function and the Gegenbauer polynomials, which are standard library functions in many numerical software packages (both are included in the GNU Scientific Library, for example). Using the results above, we see that it is easy to write down and compute a basis set for all real values of α. 3.4 The density basis functions A similar calculation can be performed to find expressions for the density basis functions, though it is simpler as the generating function can be applied immediately with no further manipulation of gn(k). First, we note the following Hankel transform (Gradshteyn & Ryzhik 2014, equation 6.623(2))   \begin{equation*} \int _0^\infty x^{\nu +1} \, \exp (-ax) \, J_\nu (xy) \, {\rm d}x = \frac{2^{\nu +1} \, \Gamma (\nu + 3/2)}{\sqrt{\pi }} \, \frac{a \, y^\nu }{(a^2 + y^2)^{\nu + 3/2}}. \end{equation*} Then, using (30) and (32) as above, we have   \begin{eqnarray} \sum _{n=0}^\infty t^n \, \rho _{nl} \propto r^{1/\alpha - 5/2} \, \sum _{n=0}^\infty \, t^n\int _0^\infty k^{\mu + 1} \, L_n^{(2\mu -1)}(2k) \, \mathrm{e}^{-k} \, J_\mu (kz) \, {\rm d}k = \frac{r^{1/\alpha - 5/2}}{\left(1-t\right)^{2\mu }} \, \int _0^\infty k^{\mu +1} \, \mathrm{e}^{-k(1+t)/(1-t)} \, J_\mu (kz) \, {\rm d}k \end{eqnarray} (38)  \begin{eqnarray} \qquad = \frac{2^{\mu +1} \, \Gamma (\mu + 3/2)}{\sqrt{\pi }} \, r^{1/\alpha - 2 + l} \, \frac{(1+t)(1-t)^2}{\left[(1+t)^2 + (1-t)^2z^2\right]^{\mu + 3/2}} \end{eqnarray} (39)  \begin{eqnarray} \qquad = \frac{2^{\mu +1} \, \Gamma (\mu + 3/2)}{\sqrt{\pi }} \, \frac{r^{1/\alpha -2+l}}{(1+r^{1/\alpha })^{3/2+\mu }} \, \sum _{n=0}^\infty t^n \, \left[C_n^{(\mu + 3/2)}(\xi ) - C_{n-1}^{(\mu + 3/2)}(\xi ) - C_{n-2}^{(\mu + 3/2)}(\xi ) + C_{n-3}^{(\mu + 3/2)}(\xi )\right]. \end{eqnarray} (40)The last line can be simplified by adding together two Gegenbauer polynomial recursion relations (Gradshteyn & Ryzhik 2014, equations 8.933(2) and 8.933(3)), resulting in   \begin{eqnarray} \rho _{nl}(r) = \frac{A_{nlm}}{16\pi \alpha ^2} \, \frac{2^{\mu +1}\Gamma (\mu + 1/2)}{\sqrt{\pi }}\frac{r^{1/\alpha -2+l}}{(1+r^{1/\alpha })^{3/2+\mu }}\left[\left(n+\mu +1/2\right)\,C_n^{(\mu +1/2)}(\xi )-\left(n+\mu -1/2\right)\,C_{n-1}^{(\mu + 1/2)}(\xi )\right]. \end{eqnarray} (41)Just as for the potential basis functions, evaluating n density basis functions requires only O(n) steps. It is now straightforward to verify that the basis functions satisfy equation (4). If the spherical harmonics are normalized to Jlm = 4π and we set Anlm = 1, the overall normalization constant for the basis functions written in equations (37) and (41) is   \begin{eqnarray} N_{nl} = \frac{\Gamma (2\mu + n)}{2^{2\mu } \, n! \, 8\pi \alpha }. \end{eqnarray} (42)For numerical purposes, it may be desirable to redefine Anlm to incorporate more of the n- and l-dependent prefactors in equations (37) and (41) (see Section 4.1). 3.5 Comparison with the Zhao and Hernquist-Ostriker basis sets Prior to our work, only a single family of biorthogonal potential-density basis functions was known (although Rahmati & Jalali 2009, provided an alternate single set of solutions with no free parameters which our results now encompass). This basis set was introduced by Zhao (1996). It is written in terms of a parameter α that affects both the inner and outer slopes of the double power law.   \begin{eqnarray} \Phi _{nl}^{\mathrm{Zhao}} = \frac{r^l}{\left(1 + r^{1/\alpha }\right)^{\mu }} \, C^{(\mu + 1/2)}_n\left(\xi \right), \ \ \ \ \rho _{nl}^{\mathrm{Zhao}} = \frac{r^{l - 2 + 1/\alpha }}{\left(1 + r^{1/\alpha }\right)^{\mu + 2}} \, C^{(\mu + 1/2)}_n\left(\xi \right), \end{eqnarray} (43)where μ and ξ are defined as in equation (37). The specific case when α = 1/2 was discovered by Clutton-Brock (1973) and has the Plummer (1911) model as its lowest order term. The case when α = 1 was first proposed by Hernquist & Ostriker (1992) and is the most widely used of all the basis-function expansions, with the Hernquist (1990) model as its lowest order term. We note immediately the close similarity in form to our previous definitions in (37) and (41), where the basis functions are also expressed in terms of the Gegenbauer polynomials. Comparing the lowest order density of the Zhao basis set with that of our new basis set,   \begin{equation*} \rho _{00} \propto \frac{1}{r^{2 - 1/\alpha }\left(1 + r^{1/\alpha }\right)^{\alpha + 3/2}}, \qquad \rho _{00}^{\mathrm{Zhao}} \propto \frac{1}{r^{2 - 1/\alpha} (1 + r^{1/\alpha })^{\alpha + 2}}. \end{equation*} Evidently, the only difference is the shallower outer slope in the former case. This is significant as popular models for dark matter haloes tend to have outer slopes closer to r−3. For example, the generalized NFW profile (Navarro et al. 2004) has ρ ∝ r−γ(1 + r)γ−3, with values for the inner slope γ ranging from 0.7 to 1.5, and the outer slope fixed to r−3. In both our and Zhao's expansion families, when the inner slope is fixed (γ ≡ 2 − 1/α), the asymptotic outer slope is then constrained. To give some examples, if we set γ = 0.7, then Zhao's outer slope is r−4.3 whereas ours is r−3.65; and if γ = 1.5 then Zhao's is r−3.5 and ours is r−3.25. We therefore expect that our basis functions will be more efficient than Zhao's for describing typical dark matter haloes, such as those resimulated by Lowing et al. (2011). Before we pass on to numerical implementation of our new family, we remark that our and Zhao's biorthonormal expansions do not by any means exhaust all the possibilities. For example, we have identified a further family of basis expansions for which the zeroth-order members have density profiles that follow Einasto-like laws (cf. Einasto 1965; Merritt et al. 2006). This will be the subject of a later paper. 4 NUMERICAL PERFORMANCE OF THE BASIS-FUNCTION EXPANSION We turn to the application of our new basis expansion to the representation of cosmological haloes. Before inspecting both analytic and numerical haloes in Sections 4.2 and 4.3, we present a formulation of the basis expansion that is computationally friendly. 4.1 Numerical implementation It is very important to take advantage of a number of recursion relations for purposes of speed. It is also more efficient to factor out the constant parts of the potential-density pairs such that only the spatially dependent pieces are calculated for each particle. We therefore modify the definition of the basis functions as   \begin{eqnarray} {\hat{\rho }}_{nlm}({\boldsymbol r}) = \hat{\rho }_{nl}(r)Y_{lm}(\theta , \phi ) \qquad {\hat{\rho }}_{nl} = \frac{r^{1/\alpha -2+l}}{(1+r^{1/\alpha })^{\mu +3/2}}\left (\left(n+\mu +\frac{1}{2}\right)C_n^{(\mu +1/2)}(\xi )-\left(n+\mu -\frac{1}{2}\right)C_{n-1}^{(\mu +1/2)}(\xi )\right), \end{eqnarray} (44)  \begin{eqnarray} {\hat{\Phi }}_{nlm}({\boldsymbol r}) = \hat{\Phi }_{nl}(r)Y_{lm}(\theta , \phi )\qquad {\hat{\Phi }}_{nl} = \frac{\mathcal {B}_\chi (\mu ,1/2)}{2\,r^{1+l}} - \frac{2 r^l}{(z^2+1)^{\mu + 1/2}} \, \sum _{j=0}^{n-1} \frac{j!\Gamma (2\mu )}{\Gamma (2\mu +j+1)} \, C_j^{(\mu + 1/2)}(\xi ), \end{eqnarray} (45)so that Poisson's equation becomes $$\nabla ^2 {\hat{\Phi }}_{nlm}=4\pi K_{nl} {\hat{\rho }}_{nlm}$$ with   \begin{eqnarray} K_{nl} = -\frac{n!\Gamma (2\mu )}{8\pi \alpha ^2\Gamma (2\mu +n)}, \end{eqnarray} (46)and the normalization constant   \begin{eqnarray} \hat{N}_{nlm} = \int \mathrm{d}^3{\boldsymbol r}\, \hat{\Phi }^{}_{nlm}({\boldsymbol r})\hat{\rho }^*_{nlm}({\boldsymbol r})=J_{lm}\frac{\alpha \sqrt{\pi }\Gamma (\mu )}{2^{1+2\mu }\Gamma (\frac{1}{2}+\mu )}, \qquad J_{lm}=\int \mathrm{d}\phi \, \mathrm{d}\theta \sin \theta \,Y_{lm}(\theta , \phi ) Y_{lm}^*(\theta , \phi ). \end{eqnarray} (47)The Gegenbauer polynomials can be constructed recursively using the relation (Gradshteyn & Ryzhik 2014, equation 8.933(1))   \begin{eqnarray} n C_n^{(w)}(\xi ) = 2(w+n-1)\xi C_{n-1}^{(w)}(\xi )-(2w+n-2)C_{n-2}^{(w)}(\xi ), \end{eqnarray} (48)where $$C_0^{(w)}(\xi )=1$$ and $$C_1^{(w)}(\xi )=2w\xi$$. To construct the potential function $${\hat{\Phi }}_{nl}$$, we define   \begin{eqnarray} A_n=\frac{2n!\Gamma (2\mu )}{\Gamma (2\mu +n+1)}, \end{eqnarray} (49)which satisfies the recurrence relation   \begin{eqnarray} \frac{A_{n+1}}{A_{n}}=\frac{n+1}{n+1+2\mu }, \ \ \ A_0=1/\mu , \end{eqnarray} (50)such that the ladder of potential functions at fixed l is given by   \begin{eqnarray} \hat{\Phi }_{nl}=\hat{\Phi }_{(n-1)l}-\frac{r^l}{(1+z^2)^{\mu +1/2}}A_{n-1}C_{n-1}^{(\mu +1/2)}(\xi ). \end{eqnarray} (51)A naive implementation of this recursion relation, wherein one builds up the higher n terms by starting from Φ0l, results in large errors for high n. In the left-hand panel of Fig. 1, we show the logarithm of the relative difference between the potential computed using this naive approach and an arbitrary precision calculation from Mathematica. The error grows with increasing n and l, as the computation requires taking a small difference between large numbers. To see this, we note that a valid series expansion of the incomplete beta function is [using Olver et al. (2016, 8.17.8) and Fields & Wimp (1961, equation 2.5)]   \begin{eqnarray} \mathcal {B}_\chi (\mu ,1/2) = \frac{\chi ^\mu }{\sqrt{1+z^2}} \sum _{j=0}^\infty A_j C_j^{(\mu +1/2)}(\xi ). \end{eqnarray} (52)A comparison with (51) shows that the terms in $$\hat{\Phi }_{nl}$$ tend to zero as n → ∞, resulting in the aforementioned catastrophic cancellation. This expression, however, provides us with an alternative method of computing $$\hat{\Phi }_{nl}$$. At some large order N, e.g. N = 2nmax, we assume that AN ≈ 0. Then we can write   \begin{eqnarray} \hat{\Phi }_{nl} \approx \frac{r^l}{(1+z^2)^{\mu +1/2}}\sum _{j=n}^{N}A_jC_j^{(\mu +1/2)}(\xi ), \end{eqnarray} (53)which requires us to calculate $$C^{(\mu +1/2)}_{n}(\xi )$$ and An recursively up to order N. We then recursively construct the potential basis functions $$\hat{\Phi }_{nl}$$ downwards using equation (51) where now all $$\hat{\Phi }_{nl}$$ are accurate to the magnitude of AN. This procedure2 results in the reduced errors shown in the right-hand panel of Fig. 1. We only perform this procedure for l > 4 as for l ≤ 4 the naive implementation is satisfactory and the decay of $$\hat{\Phi }_{nl}$$ is weak for small l. Figure 1. View largeDownload slide Relative error in potential basis-function calculation using two methods: on the left, the result of using a forward recursion scheme and on the right using a combination of forward and backward recursion. The relative error is the difference between the potential computed by the basis expansion and an arbitrary precision calculation from Mathematica computed at r/rs = 1.3 using α = 1.2. Note that the floor on the error in the right-hand panel is set by the precision of our Mathematica calculation and can be lowered if so desired. We only perform the ‘forwards+backwards’ procedure for l > 4. Figure 1. View largeDownload slide Relative error in potential basis-function calculation using two methods: on the left, the result of using a forward recursion scheme and on the right using a combination of forward and backward recursion. The relative error is the difference between the potential computed by the basis expansion and an arbitrary precision calculation from Mathematica computed at r/rs = 1.3 using α = 1.2. Note that the floor on the error in the right-hand panel is set by the precision of our Mathematica calculation and can be lowered if so desired. We only perform the ‘forwards+backwards’ procedure for l > 4. Finally, note that the $$\hat{N}_{nl}$$ are independent of n and the Knl satisfy the properties   \begin{eqnarray} \frac{K_{(n+1)l}}{K_{nl}}=\frac{n+1}{n+2\mu }, \ \ \ K_{0l}=-\frac{1}{8\pi \alpha ^2}. \end{eqnarray} (54)With these definitions, a set of coefficients $$\hat{C}_{nlm}$$ is computed from a cloud of particles, and the potential and density reconstructed as   \begin{eqnarray} \hat{C}_{nlm} = \frac{1}{\hat{N}_{nlm}K_{nl}}\sum _i m_i \hat{\Phi }_{nl}(r_i)Y_{lm}(\theta _i, \phi _i), \qquad \Phi = \sum _{nlm}\hat{C}_{nlm}\hat{\Phi }_{nlm}, \qquad \rho = \sum _{nlm}K_{nl}\hat{C}_{nlm}\hat{\rho }_{nlm}. \end{eqnarray} (55) 4.2 Analytical haloes 4.2.1 Spherical NFW models We first consider the reconstruction of spherical NFW haloes using the radial terms of the expansion. We compare the α = 1 member of our family of expansions with the Hernquist & Ostriker (1992, hereafter HO) expansion, which is the α = 1 member of the Zhao (1996) family. Both basis sets have lowest order densities with 1/r cusps as r → 0. In Fig. 2, we show the expansion of a spherical NFW potential and density, using only radial terms (n ≥ 0, l = 0). With an equal number of terms, our α = 1 basis set performs better than the corresponding Hernquist-Ostriker set due to its closer approximation to the NFW profile in the asymptotic fall-off of the lowest order density, ρ ∼ r−7/2. The corresponding behaviour of the Hernquist-Ostriker basis set is ρ ∼ r−4. Note that the r-axis of Fig. 2 is logarithmically scaled and measured in units of the scale length, so we in fact get several hundred additional scale lengths of accurate behaviour. The improved convergence at large radii also reduces the amplitudes of the oscillations at smaller radii. The expansion of the potential is always more accurate than that of the density because the oscillations are integrated over, and therefore effectively smoothed; the accuracy of the radial acceleration expansion lies between that of the potential and the density. According to Fig. 3, the coefficients used in the spherical NFW expansion empirically follow power laws for both basis sets. In our expansion, Cn00 ∼ n−2 whereas in the Hernquist-Ostriker expansion Cn00 ∼ n−3. Figure 2. View largeDownload slide Reconstruction of a spherical NFW potential (left) and density (right) with the new super-NFW basis set (blue) and the Hernquist-Ostriker basis set (red). The distance is given in units of the NFW scale length. Both expansions use radial terms up to n = 20 and no angular terms (l = 0). Both expansions oscillate around the true value, which is represented by the horizontal grey line. Figure 2. View largeDownload slide Reconstruction of a spherical NFW potential (left) and density (right) with the new super-NFW basis set (blue) and the Hernquist-Ostriker basis set (red). The distance is given in units of the NFW scale length. Both expansions use radial terms up to n = 20 and no angular terms (l = 0). Both expansions oscillate around the true value, which is represented by the horizontal grey line. Figure 3. View largeDownload slide The run of the radial expansion coefficients with n for a spherical NFW model using the new basis set (blue) and the Hernquist-Ostriker basis set (red). Also plotted as dashed lines are the n−2 and n−3 curves, suggesting that the coefficients fall off asymptotically like n−2 in our case and like n−3 for the Hernquist-Ostriker case. Figure 3. View largeDownload slide The run of the radial expansion coefficients with n for a spherical NFW model using the new basis set (blue) and the Hernquist-Ostriker basis set (red). Also plotted as dashed lines are the n−2 and n−3 curves, suggesting that the coefficients fall off asymptotically like n−2 in our case and like n−3 for the Hernquist-Ostriker case. To make quantitative statements about the error, we follow Vasiliev (2013) and calculate the integrated squared error. This is the mass-weighted fractional density error defined by   \begin{eqnarray} ISE = \int _{r > r_\mathrm{min} \atop r < r_\mathrm{max}} \frac{\left|\rho _\mathrm{exact} - \rho _\mathrm{approx}\right|^2}{\rho _\mathrm{exact}} \, {\rm d}^{\,3}{\boldsymbol r}. \end{eqnarray} (56)Here, ρapprox is the density reconstructed using basis functions up to a radial order of nmax. Although other measures of error can be constructed, Vasiliev's suggestion is appealing as it is mass-weighted and does not bias the result towards the outer parts of the model. We use rmin = 0.01 and rmax = 100, as this is the range over which we expect the expansions to be most applicable in astrophysical problems. The run of this error measure with nmax is shown on the left in Fig. 4 for the two spherical NFW expansions under consideration. It is clear both that higher values of α give better accuracy, and that for given α and nmax our basis set is more accurate than Zhao's in the task of reproducing NFW haloes. To illustrate this difference, on the right in Fig. 4 we show how the optimum α (the minimum of each plot) varies with differing nmax, with our basis set performing better than Zhao's in each and every case. Figure 4. View largeDownload slide Left: The integrated squared error as defined in equation (56) between the exact NFW profile and the profile reconstructed with radial basis functions up to order nmax for different values of α in Zhao's basis set (red lines) and the new basis set in this paper (blue lines). Right: The integrated squared error between the exact NFW profile and the profile reconstructed with radial basis functions, against the parameter α that characterises each basis set for different numbers of radial basis functions. Again, red lines use Zhao's basis set and blue lines use the new basis set. Figure 4. View largeDownload slide Left: The integrated squared error as defined in equation (56) between the exact NFW profile and the profile reconstructed with radial basis functions up to order nmax for different values of α in Zhao's basis set (red lines) and the new basis set in this paper (blue lines). Right: The integrated squared error between the exact NFW profile and the profile reconstructed with radial basis functions, against the parameter α that characterises each basis set for different numbers of radial basis functions. Again, red lines use Zhao's basis set and blue lines use the new basis set. A detailed error analysis of basis-function techniques is carried out in Kalapotharakos, Efthymiopoulos & Voglis (2008), who claim that a significant factor in obtaining high accuracy is choosing a basis set whose lowest order density function has the correct inner slope behaviour. For a spherical NFW model, this would imply that the α = 1 expansion is best, as the zeroth-order model behaves like ρ ∼ r−1 at small radii. It is possible to choose combinations of the density basis functions so that the most diverging terms at the centre cancel, as HO already showed by building a cored-density model from the n = 0 and n = 1 monopole terms of their basis set. This delicate balance though is susceptible to numerical error. To reproduce accurately the precession of orbits that pass close to the centre, the criterion of Kalapotharakos et al. (2008) seems reasonable. However, it is at odds with the plots in Fig. 4, which show that models with α ≈ 3 or ρ ∼ r−5/3 provide the smallest integrated square error in the density, if nmax ≈ 20. A similar result was obtained by Vasiliev (2013), who used an identical integrated error measure to the present work and concluded that a somewhat higher value of α is preferable for providing a global fit to cosmological haloes, even at the expense of accuracy of the inner slope exponent near the centre. Clearly, the best choice depends on the application in hand. 4.2.2 Flattened NFW models It is important to consider flattened objects, as the haloes found in N-body simulations are generically flattened or triaxial (e.g. Jing & Suto 2002; Allgood et al. 2006). We test the performance by attempting to reconstruct a flattened density profile and including the angular terms or (l, m) terms in the series. An axisymmetric NFW density profile is $$\rho = {\bar{m}}^{-1}(1+ \bar{m})^{-2}$$, where $${\bar{m}}^2 = x^2 + y^2 + z^2/q^2$$. This means the density is stratified on similar concentric spheroids with axis ratio q. Figs 5 and 6 show the expansion of a flattened (q = 0.8) NFW density, and the corresponding acceleration due to the potential. In each case, we compare the α = 1 member of our family of expansions with the Hernquist-Ostriker expansion, which is the α = 1 member of the Zhao (1996) family. Typically, we use nmax = 20 radial basis functions and lmax = 12 angular basis functions. Figure 5. View largeDownload slide Expansion of a flattened (q = 0.8) NFW density. Viewing angle θ is shown on each plot. Both expansions use radial terms up to n = 20 and angular terms up to l = 12. The error measure is Δρ ≡ log |1 − ρapprox/ρexact| (lower is better; the dips are due to the oscillations around the true value). Figure 5. View largeDownload slide Expansion of a flattened (q = 0.8) NFW density. Viewing angle θ is shown on each plot. Both expansions use radial terms up to n = 20 and angular terms up to l = 12. The error measure is Δρ ≡ log |1 − ρapprox/ρexact| (lower is better; the dips are due to the oscillations around the true value). Figure 6. View largeDownload slide Expansion of the radial acceleration in a flattened (q = 0.8) NFW potential. Viewing angle θ is shown on each plot. Both expansions use radial terms up to n = 20 and angular terms up to l = 12. The error measure is $$\Delta _{a_r} \equiv \log {\left|1 - a_{r\,\mathrm{approx}}/a_{r\,\mathrm{exact}}\right|}$$. Figure 6. View largeDownload slide Expansion of the radial acceleration in a flattened (q = 0.8) NFW potential. Viewing angle θ is shown on each plot. Both expansions use radial terms up to n = 20 and angular terms up to l = 12. The error measure is $$\Delta _{a_r} \equiv \log {\left|1 - a_{r\,\mathrm{approx}}/a_{r\,\mathrm{exact}}\right|}$$. Each of the reconstructed quantities is plotted along three polar angles (θ). Note that the convergence is always superior nearer the equator (θ = π/2) than at the poles (θ = 0). This is a feature of any expansion involving spherical harmonics and can be remedied by introducing additional ‘l’ terms. In the flattened case, both expansions lose accuracy in the very inner and outer parts of the haloes. As the l-dependence of both our α = 1 set and the HO set is similar, we do not expect either basis set to be favoured in this regard. However, the superior behaviour in the outskirts of our basis set is maintained. Lowing et al. (2011) used the HO basis set to represent haloes, where on the order of tens of terms are used in both the angular and radial directions. They found errors of < 10 per cent are achieved over a few hundred kiloparsecs; we expect our basis set to provide improved accuracy in this regime. 4.3 Numerical haloes We now analyse a collection of ten Milky Way-like dark matter haloes, extracted from a suite of cosmological N-body zoom-in simulations. These simulations are run with the N-body part of gadget-3 which is similar to gadget-2 (Springel 2005). The zoom-in strategy follows Oñorbe et al. (2014) and all initial conditions are generated with music (Hahn & Abel 2011). Cosmological parameters are taken from the Planck Collaboration XVI et al. (2014) with h = 0.679, Ωb = 0.0481, Ω0 = 0.306, $$\Omega _\Lambda =0.694$$, σ8 = 0.827, and ns = 0.962. In order to select our haloes, we first simulate a 50 h−1 Mpc box with 5123 particles from z = 50 to 0. We use rockstar (Behroozi, Wechsler & Wu 2013) to identify haloes and select Milky Way-like haloes which have virial masses between 7.5 × 1011and2 × 1012 M⊙, no major mergers since z = 1, and no haloes with half the mass of Milky Way analogue's mass within 2 h−1 Mpc. For each selected halo, we select all particles within 10 virial radii and run an intermediate resolution zoom-in whose maximum resolution is 20483, corresponding to a particle mass of 1.8 × 106 M⊙. This intermediate step helps to reduce the contamination from low-resolution particles in our final, high-resolution zoom-in. For the final zoom-in, we take the intermediate resolution simulation and select all particles within 7.5 virial radii. We then run a zoom-in with a maximum resolution of 40963, corresponding to a particle mass of 2.23 × 105 M⊙. All of our high-resolution zoom-ins are uncontaminated within 1 h−1 Mpc of the main halo. The detailed properties of each halo are given in Table 1. Table 1. The properties of each numerical halo in our sample. The haloes are listed in the order that they are displayed in the figures. N is the number of particles, Mv is the virial mass, rv is the virial radius, rs is the scale length, riso is the distance at which the slope is isothermal (so riso = rs/(1/2 + α)α), α parametrizes the inner and outer slopes of the density basis functions, p is the y–x axis ratio, and q is the z–x axis ratio. N/106  Mv/(1012 M⊙)  rv/kpc  rs/kpc  riso/kpc  α  p  q  Figure  10.3  1.88  325  65.2  33.6  1.22  0.860  0.811  7  5.0  0.91  256  37.8  20.5  1.18  0.972  0.816  7  8.7  1.57  306  47.6  25.8  1.18  0.889  0.761  7  7.4  1.36  292  52.8  25.9  1.26  0.800  0.733  7  5.9  1.07  269  58.1  30.7  1.20  0.804  0.795  7  5.5  0.89  254  46.0  23.7  1.22  0.909  0.834  7  4.9  0.89  253  35.7  19.8  1.16  0.807  0.780  8  11.3  1.62  310  123  38.1  1.59  0.858  0.769  8  7.7  1.11  273  98.5  32.9  1.54  0.926  0.779  8  7.6  1.71  315  132.6  47.6  1.49  0.861  0.730  10  N/106  Mv/(1012 M⊙)  rv/kpc  rs/kpc  riso/kpc  α  p  q  Figure  10.3  1.88  325  65.2  33.6  1.22  0.860  0.811  7  5.0  0.91  256  37.8  20.5  1.18  0.972  0.816  7  8.7  1.57  306  47.6  25.8  1.18  0.889  0.761  7  7.4  1.36  292  52.8  25.9  1.26  0.800  0.733  7  5.9  1.07  269  58.1  30.7  1.20  0.804  0.795  7  5.5  0.89  254  46.0  23.7  1.22  0.909  0.834  7  4.9  0.89  253  35.7  19.8  1.16  0.807  0.780  8  11.3  1.62  310  123  38.1  1.59  0.858  0.769  8  7.7  1.11  273  98.5  32.9  1.54  0.926  0.779  8  7.6  1.71  315  132.6  47.6  1.49  0.861  0.730  10  View Large For each halo, we take all of the particles within 500 kpc of the main halo in z = 0 snapshot. This corresponds to between 5–12 million particles, depending on the halo mass. We wish to investigate the ability of our new basis sets to represent these numerical haloes. To this end, we must first choose the two global parameters that specify the basis set – the scale length rs and the parameter α. We need not worry at this stage about the overall normalization (related to the total mass) because this will be set automatically when performing the sum over particles via equation (10). To set rs and α, we first fit the zeroth-order density function ρ000, which is a spherically symmetric model. Because we know the virial radius rv for each halo, we perform the fitting procedure in terms of the dimensionless concentration c ≡ rv/rs, as is standard in the literature. The particles are binned logarithmically in radius, and a non-linear least-squares algorithm then adjusts α and c to minimize the difference between the logarithms of the inferred bin density and the model density. For this fitting procedure, we use a form of the zeroth-order density function that is parametrized by the mass enclosed by the virial radius, i.e. which satisfies $$\int _0^{r_{\rm v}} 4\pi \rho (r) r^2 {\rm d}r = M_{\rm v}$$. This is   \begin{eqnarray} \rho (r) = \left[\mathcal {B}_{\chi _v}(\alpha ,1/2) - \frac{c/\alpha }{\left(1 + c^{1/\alpha }\right)^{\alpha +1/2}}\right]^{-1} \, \frac{M_{\rm v} \, (\alpha + 1/2)}{4 \pi \alpha ^2 \, (r_{\rm v}/c)^3} \, \frac{(cr/r_{\rm v})^{1/\alpha - 2}}{\left(1 + (cr/r_{\rm v})^{1/\alpha }\right)^{\alpha + 3/2}}, \end{eqnarray} (57)where χv ≡ c1/α/(1 + c1/α). We can now perform a full expansion, using for each halo the basis set with the determined best values of α and rs. To compare the accuracy in the reproduction of the density, we draw contour plots in each principal plane, as shown in Figs 7, 8, and 10. The smooth underlying density distribution in each principal plane of the original numerical halo is estimated by taking particles in a slab of width 2 kpc around the plane, then creating a two-dimensional histogram with each bin having an area of approximately 3 kpc2. The density as reconstructed from the basis-function expansion is sampled along rays in each principal plane, and then interpolated on to a grid. The same 12 contours are drawn on each plot, spaced approximately logarithmically between 106 and 2 × 108  M⊙ kpc−3. Figure 7. View largeDownload slide From left to right: Density contours in the three principal planes of the original numerical halo, then density contours in the same three principal planes of the halo as reconstructed with basis functions up to n = 20 and l = 12. Distances are in kpc, and the densities are in  M⊙kpc−3. Figure 7. View largeDownload slide From left to right: Density contours in the three principal planes of the original numerical halo, then density contours in the same three principal planes of the halo as reconstructed with basis functions up to n = 20 and l = 12. Distances are in kpc, and the densities are in  M⊙kpc−3. Figure 8. View largeDownload slide Continuation of Fig. 7. The latter two haloes possess a large number of small satellites, so for these we display on the right two reconstructions: top, a reconstruction of the full halo; and bottom, one where the particles energetically bound to the 50 most massive satellites have been removed. This illustrates the issues surrounding expansion accuracy in the presence of a sub-structure. Figure 8. View largeDownload slide Continuation of Fig. 7. The latter two haloes possess a large number of small satellites, so for these we display on the right two reconstructions: top, a reconstruction of the full halo; and bottom, one where the particles energetically bound to the 50 most massive satellites have been removed. This illustrates the issues surrounding expansion accuracy in the presence of a sub-structure. The haloes displayed in Figs 7 and 8 are rotated such that the principal planes are aligned with the coordinate axes (shortest axis along the z-axis), so that the angular expansion coefficients can be compared meaningfully. Distributions of the expansion coefficients for these nine numerical haloes are shown in Fig. 9. One could produce an artificial ‘halo’ with geometry typical for this family of nine numerical haloes by drawing coefficients from the distributions shown in these plots. From Figs 7 and 8, we see that key features of the haloes are resolved correctly, including: (i) the orientation of the principal axes, (ii) the axis ratios in the three principal planes, and (iii) the run of ellipticity with radius. The size of the smallest resolvable feature is limited by the distance between the roots of the polynomial used to define the highest order function used in the reconstruction. Figure 9. View largeDownload slide We denote normalized coefficients by tildes ($$\tilde{C}_{nlm} \equiv C_{nlm}/C_{000}$$) and the standard deviation of a variable x by σ(x). Panels 1 and 2: the radial coefficients $$\tilde{C}_{n00}$$ for nine numerical haloes. Panel 3: the angular coefficients $$\tilde{C}_{0l0}$$. Panel 4: the angular coefficients $$\tilde{C}_{0ll}$$. The latter three plots are normalized by each coefficient's standard deviation to make details easier to see. The trend in the C0l0 coefficients indicates the flattening of the haloes along the z-axis, and the trend in the $$\tilde{C}_{0ll}$$ coefficients shows the elongation along the x-axis. Figure 9. View largeDownload slide We denote normalized coefficients by tildes ($$\tilde{C}_{nlm} \equiv C_{nlm}/C_{000}$$) and the standard deviation of a variable x by σ(x). Panels 1 and 2: the radial coefficients $$\tilde{C}_{n00}$$ for nine numerical haloes. Panel 3: the angular coefficients $$\tilde{C}_{0l0}$$. Panel 4: the angular coefficients $$\tilde{C}_{0ll}$$. The latter three plots are normalized by each coefficient's standard deviation to make details easier to see. The trend in the C0l0 coefficients indicates the flattening of the haloes along the z-axis, and the trend in the $$\tilde{C}_{0ll}$$ coefficients shows the elongation along the x-axis. Two haloes in Fig. 8 demonstrate one pitfall of the method. The presence of unresolved but massive subhaloes can cause a blow-up in the higher order coefficients of the expansion. This would be cancelled by yet higher order terms, but as the expansion is truncated these are not present. This problem is generic – it affects other basis sets, such as Zhao's, to a greater or lesser degree – but it can be remedied by removing the unresolved subhaloes by an automated halo-finding procedure, as demonstrated in the figure. When these subhaloes are removed, the ‘optimal’ values of α and rs (as determined by fitting equation (57)) are reduced to values more characteristic of the other haloes in the sample. The final halo, displayed in Fig. 10, provides a serious challenge. It is accompanied by a massive close-in satellite with a mass of 1.9 × 1011 M⊙, and thus has a highly aspherical geometry. We therefore examine this halo in greater detail in Fig. 10 using nmax = 20 and lmax = 12 (middle panel) and nmax = 40 and lmax = 40 (lower panel). Remarkably, the overall structure of both the halo and the large satellite is correctly resolved even with nmax = 20. This is impressive, as we might have suspected at the outset that two basis-function expansions, centred on each object, would be necessary to reproduce the merging structure. Figure 10. View largeDownload slide A halo with a prominent, massive satellite. From top to bottom: the original halo; the reconstruction with maximum order n = 20 and l = 12; the same with n = 40 and l = 40. Figure 10. View largeDownload slide A halo with a prominent, massive satellite. From top to bottom: the original halo; the reconstruction with maximum order n = 20 and l = 12; the same with n = 40 and l = 40. 5 CONCLUSIONS Biorthogonal density and potential basis functions provide useful and flexible ways of describing realistic dark matter haloes and galaxies, which may be aspherical, triaxial, or further misshapen. The coefficients of these basis-function expansions can be found easily by summing over the particles in an N-body realization and used to reconstruct both the potential and density. We have discovered a completely new family of biorthogonal potential-density pairs, parametrized in terms of α. The zeroth-order model has a density ρ ∼ r−2+1/α at small radii and ρ ∼ r−3−1/(2α) at large radii. This double-power-law profile has a central logarithmic slope between 0 and −2 and an asymptotic logarithmic slope between −3 and −4, making it perfect for representing dark haloes. The zeroth-order model with α = 1/2 has a harmonic core and is the celebrated perfect sphere of de Zeeuw (1985), whilst the model with α = 1 has a 1/r central density cusp and is the super-NFW model (Lilley et al. 2017). For each of these zeroth-order models, we provide a biorthogonal basis-function expansion in terms of standard special functions readily available in numerical libraries. This extends the zeroth-order model into the highly realistic regime of flattened, distorted, and triaxial density profiles. Previously, the only known family of biorthogonal potential-density pairs was the one outlined by Zhao (1996), of which the most widely used member is the HO expansion. The zeroth-order model has a density ρ ∼ r−2+1/α at small radii and ρ ∼ r−3−1/α at large radii. Although the central density cusp is the same, the outer density falls off rather more quickly than in our expansion, making Zhao's family less well-matched to modelling dark haloes. We have demonstrated the capabilities of our basis-function expansions by using them to recreate spherical and flattened dark haloes with analytic densities of Navarro et al. (1997) form. In particular, we showed that our method represents a significant improvement over the HO expansion, giving a more accurate reproduction of the density, potential, and radial force of NFW-like models. Additionally, we decomposed 10 simulated cosmological haloes using our basis functions, computing the coefficients as simple sums over the particles. This yielded very encouraging results with highly flattened and triaxial dark-halo density distributions well reproduced with typically 20 radial and 12 angular terms (giving a total of 20 × 122 ≈ 3000 terms). Simulated dark haloes can be lopsided, distorted or twisted, especially if there is a nearby companion exerting strong tidal forces (cf. the Milky Way and the Magellanic Clouds). Particular striking is the ability of our basis-function expansion to reproduce the density of a highly asymmetric dark halo which is in the process of accreting a companion large subhalo. Lowing et al. (2011) used the HO expansion to create potentials approximating the structure of the dark haloes in the Aquarius simulations. They found that the radial force from the expansion sometimes differed from that in the simulations when the central logarithmic slope of the numerical halo density was shallower than the slope of −1 of the zeroth-order model. They also highlighted a new application for biorthogonal basis-function expansions as a tool for data compression of large and expensive N-body simulations. By calculating sets of coefficients from numerical snapshots at different redshifts and interpolating between them, they obtained a compact and efficient summary of the halo's evolution throughout time. This is very powerful because new objects can then be inserted and the entire simulation can be re-run cheaply making the method well-suited to study a range of problems, such as the dynamics of small satellites, the stripping of globular clusters or the shape of tidal streams. Another important area of application is to the fitting of data on the Milky Way galaxy, which has increased substantially in both quality and quantity over the last few years. Models of the Milky Way galaxy are assembled from simple building blocks. Results of calculations using such ‘pieces of Lego’ are often very troubling. For example, when the position and velocity data of stars in the Sagittarius stream are fitted to such models, the conclusion is that the dark halo is triaxial with the short and long axes in the Galactic plane (Law & Majewski 2010). This configuration is unstable (see e.g. Debattista et al. 2013) and in conflict with observational data on the disc (Kuijken & Tremaine 1994). The strong suspicion is that the inflexible model of the Milky Way's potential prevented proper exploration of parameter space and artificially confined the solution for the dark matter distribution into an unrealistic straitjacket. A completely new way of representing the Galactic dark halo is needed with the advent of large-scale photometric, spectroscopic and astrometric surveys of the Galaxy. For model fitting, the dark-halo potential should be represented by distributions of the coefficients that can be used as priors in Bayesian inference from the data, rather than say a single number (the flattening) in a predetermined and unadaptable density law. Of course, the shape of a dark halo depends on the nature of the dark matter particle and on the extent of feedback processes (see e.g. Sellwood 2004; Macciò et al. 2012). The shape also depends on a host of other factors, including the mass of the halo, its environment (isolated versus group), its recent history (e.g. late infall of a large subhalo), and the presence or absence of a disc. It will be interesting to test our new basis-function expansion method on the full variety of numerically constructed haloes, and understand how the distribution of coefficients changes with the underlying physics. The main impediment to efficient exploitation of these ideas is that so few biorthogonal pairs are known. Our new discovery helps in this regard, but it has not exhausted the supply of such expansions. We will show in a later paper how to extend the methods in this paper to other cosmologically inspired dark-halo density laws. Acknowledgements JLS and EJL acknowledge the support of the STFC. We thank members of the Institute of Astronomy Streams group for discussions and comments as this work was in progress. The anonymous referee is thanked for a number of helpful suggestions. Footnotes 1 Who were in turn motivated by the analogous disc case considered in Clutton-Brock (1972). 2 It is analogous to Miller's method in numerical analysis, apparently introduced by J.C.P. Miller for the computation of Bessel functions (see e.g. Abramowitz & Stegun 1965). REFERENCES Abramowitz M., Stegun I. A., 1965, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables . Dover Publ., New York Allen A. J., Palmer P. L., Papaloizou J., 1990, MNRAS , 242, 576 https://doi.org/10.1093/mnras/242.4.576 CrossRef Search ADS   Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS , 367, 1781 https://doi.org/10.1111/j.1365-2966.2006.10094.x CrossRef Search ADS   An J. H., Evans N. W., 2006, AJ , 131, 782 https://doi.org/10.1086/499305 CrossRef Search ADS   An J., Zhao H., 2013, MNRAS , 428, 2805 https://doi.org/10.1093/mnras/sts175 CrossRef Search ADS   Aoki S., Iye M., 1978, PASJ , 30, 519 Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ , 762, 109 https://doi.org/10.1088/0004-637X/762/2/109 CrossRef Search ADS   Bertin G., 2014, Dynamics of Galaxies . Cambridge Univ. Press, Cambridge Google Scholar CrossRef Search ADS   Binney J., Tremaine S., 1987, Galactic Dynamics . Princeton Univ. Press, Princeton, NJ, p. 747 Clutton-Brock M., 1972, Ap&SS , 16, 101 https://doi.org/10.1007/BF00643095 CrossRef Search ADS   Clutton-Brock M., 1973, Ap&SS , 23, 55 https://doi.org/10.1007/BF00647652 CrossRef Search ADS   de Zeeuw T., 1985, MNRAS , 216, 273 https://doi.org/10.1093/mnras/216.2.273 CrossRef Search ADS   Debattista V. P., Roškar R., Valluri M., Quinn T., Moore B., Wadsley J., 2013, MNRAS , 434, 2971 https://doi.org/10.1093/mnras/stt1217 CrossRef Search ADS   Dehnen W., 2014, Comput. Astrophys. Cosmol. , 1, 1 https://doi.org/10.1186/s40668-014-0001-7 CrossRef Search ADS   Diemand J., Zemp M., Moore B., Stadel J., Carollo C. M., 2005, MNRAS , 364, 665 https://doi.org/10.1111/j.1365-2966.2005.09601.x CrossRef Search ADS   Diemer B., Mansfield P., Kravtsov A. V., More S., 2017, ApJ , 843, 140 https://doi.org/10.3847/1538-4357/aa79ab CrossRef Search ADS   Einasto J., 1965, Tr. Astrofizicheskogo Inst. Alma-Ata , 5, 87 Evans N. W., 1993, MNRAS , 260, 191 https://doi.org/10.1093/mnras/260.1.191 CrossRef Search ADS   Evans N. W., An J. H., 2005, MNRAS , 360, 492 https://doi.org/10.1111/j.1365-2966.2005.09078.x CrossRef Search ADS   Evans N. W., An J. H., 2006, Phys. Rev. D , 73, 023524 https://doi.org/10.1103/PhysRevD.73.023524 CrossRef Search ADS   Evans N. W., Read J. C. A., 1998, MNRAS , 300, 106 https://doi.org/10.1046/j.1365-8711.1998.01864.x CrossRef Search ADS   Fields J. L., Wimp J., 1961, Math. Comput. , 15, 390 https://doi.org/10.1090/S0025-5718-1961-0125992-3 CrossRef Search ADS   Fridman A. M., Polyachenko V. L., 1984, Physics of Gravitating Systems I: Equilibrium and Stability . Springer Verlag, Berlin. Gradshteyn I., Ryzhik I., 2014, Table of Integrals, Series and Products , 8th edn. Academic Press, New York Hahn O., Abel T., 2011, MNRAS , 415, 2101 https://doi.org/10.1111/j.1365-2966.2011.18820.x CrossRef Search ADS   Hernquist L., 1990, ApJ , 356, 359 https://doi.org/10.1086/168845 CrossRef Search ADS   Hernquist L., Ostriker J. P., 1992, ApJ , 386, 375 (HO) https://doi.org/10.1086/171025 CrossRef Search ADS   Jing Y. P., Suto Y., 2002, ApJ , 574, 538 https://doi.org/10.1086/341065 CrossRef Search ADS   Kalapotharakos C., Efthymiopoulos C., Voglis N., 2008, MNRAS , 383, 971 https://doi.org/10.1111/j.1365-2966.2007.12417.x CrossRef Search ADS   Kuijken K., Tremaine S., 1994, ApJ , 421, 178 https://doi.org/10.1086/173635 CrossRef Search ADS   Kuzmin G. G., 1956, Publ. Tartu Astrofizica Obs. , 33, 75 Law D. R., Majewski S. R., 2010, ApJ , 714, 229 https://doi.org/10.1088/0004-637X/714/1/229 CrossRef Search ADS   Lilley E. J., Evans N. W., Sanders J. L., 2017, MNRAS , submitted Lowing B., Jenkins A., Eke V., Frenk C., 2011, MNRAS , 416, 2697 https://doi.org/10.1111/j.1365-2966.2011.19222.x CrossRef Search ADS   Macciò A. V., Paduroiu S., Anderhalden D., Schneider A., Moore B., 2012, MNRAS , 424, 1105 https://doi.org/10.1111/j.1365-2966.2012.21284.x CrossRef Search ADS   Merritt D., Graham A. W., Moore B., Diemand J., Terzić B., 2006, AJ , 132, 2685 https://doi.org/10.1086/508988 CrossRef Search ADS   Moore B., Governato F., Quinn T., Stadel J., Lake G., 1998, ApJ , 499, L5 https://doi.org/10.1086/311333 CrossRef Search ADS   Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ , 490, 493 https://doi.org/10.1086/304888 CrossRef Search ADS   Navarro J. F. et al.  , 2004, MNRAS , 349, 1039 https://doi.org/10.1111/j.1365-2966.2004.07586.x CrossRef Search ADS   Ngan W., Bozek B., Carlberg R. G., Wyse R. F. G., Szalay A. S., Madau P., 2015, ApJ , 803, 75 https://doi.org/10.1088/0004-637X/803/2/75 CrossRef Search ADS   Ngan W., Carlberg R. G., Bozek B., Wyse R. F. G., Szalay A. S., Madau P., 2016, ApJ , 818, 194 https://doi.org/10.3847/0004-637X/818/2/194 CrossRef Search ADS   Olver F. W. J., Olde Daalhuis A. B., Lozier D. W., Schneider B. I., Boisvert R. F., Clark C. W., Miller B. R., Saunders B. V., 2016, NIST Digital Library of Mathematical Functions . Available at: http://dlmf.nist.gov/ Oñorbe J., Garrison-Kimmel S., Maller A. H., Bullock J. S., Rocha M., Hahn O., 2014, MNRAS , 437, 1894 https://doi.org/10.1093/mnras/stt2020 CrossRef Search ADS   Palmer P. L., 1994, Astrophysics and Space Science Library, Vol. 185, Stability of Collisionless Stellar Systems: Mechanisms for the Dynamical Structure of Galaxies . Springer, Berlin Planck Collaboration XVI et al.  , 2014, A&A , 571, A16 https://doi.org/10.1051/0004-6361/201321591 CrossRef Search ADS   Plummer H. C., 1911, MNRAS , 71, 460 https://doi.org/10.1093/mnras/71.5.460 CrossRef Search ADS   Polyachenko V. L., Shukhman I. G., 1981, Sov. Astron. , 25, 533 Rahmati A., Jalali M. A., 2009, MNRAS , 393, 1459 https://doi.org/10.1111/j.1365-2966.2008.14226.x CrossRef Search ADS   Saha P., 1991, MNRAS , 248, 494 https://doi.org/10.1093/mnras/248.3.494 CrossRef Search ADS   Sellwood J. A., 2004, in Ryder S., Pisano D., Walker M., Freeman K., eds, IAU Symp. 220, Dark Matter in Galaxies . Kluwer, Dordrecht, p. 27 Springel V., 2005, MNRAS , 364, 1105 https://doi.org/10.1111/j.1365-2966.2005.09655.x CrossRef Search ADS   Toomre A., 1963, ApJ , 138, 385 https://doi.org/10.1086/147653 CrossRef Search ADS   Vasiliev E., 2013, MNRAS , 434, 3174 https://doi.org/10.1093/mnras/stt1235 CrossRef Search ADS   Weinberg M. D., 1989, MNRAS , 239, 549 https://doi.org/10.1093/mnras/239.2.549 CrossRef Search ADS   Zhao H., 1996, MNRAS , 278, 488 https://doi.org/10.1093/mnras/278.2.488 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# Galaxy halo expansions: a new biorthogonal family of potential-density pairs

, Volume 476 (2) – May 1, 2018
18 pages

/lp/ou_press/galaxy-halo-expansions-a-new-biorthogonal-family-of-potential-density-l4B8eYI1IN
Publisher
Oxford University Press
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty296
Publisher site
See Article on Publisher Site

### Abstract

Abstract Efficient expansions of the gravitational field of (dark) haloes have two main uses in the modelling of galaxies: first, they provide a compact representation of numerically constructed (or real) cosmological haloes, incorporating the effects of triaxiality, lopsidedness or other distortion. Secondly, they provide the basis functions for self-consistent field expansion algorithms used in the evolution of N-body systems. We present a new family of biorthogonal potential-density pairs constructed using the Hankel transform of the Laguerre polynomials. The lowest order density basis functions are double-power-law profiles cusped like ρ ∼ r−2+1/α at small radii with asymptotic density fall-off like ρ ∼ r−3−1/(2α). Here, α is a parameter satisfying α ≥ 1/2. The family therefore spans the range of inner density cusps found in numerical simulations, but has much shallower – and hence more realistic – outer slopes than the corresponding members of the only previously known family deduced by Zhao and exemplified by Hernquist & Ostriker. When α = 1, the lowest order density profile has an inner density cusp of ρ ∼ r−1 and an outer density slope of ρ ∼ r−3.5, similar to the famous Navarro, Frenk & White (NFW) model. For this reason, we demonstrate that our new expansion provides a more accurate representation of flattened NFW haloes than the competing Hernquist-Ostriker expansion. We utilize our new expansion by analysing a suite of numerically constructed haloes and providing the distributions of the expansion coefficients. methods: numerical, galaxies: kinematics and dynamics 1 INTRODUCTION A galaxy is often modelled as a sum of simple analytical components, such as a Hernquist (1990) bulge or a Navarro, Frenk & White (1997, NFW) or axisymmetric logarithmic halo (Binney & Tremaine 1987; Evans 1993). Although crude, such an approach produces manageable zeroth-order models from which further insights can be garnered. However, in the era of powerful numerical simulations of galaxy formation as well as accurate positional and kinematic data in our own Milky Way, such broad-brush models fail to capture the complexity of structures observed in both simulations and data. A richer and more flexible way of representing the gravitational field of galaxies (and their components) is to use basis-function or halo-expansion techniques (see e.g. Clutton-Brock 1972; Hernquist & Ostriker 1992; Zhao 1996; Lowing et al. 2011). In this framework, small higher order deviations away from the zeroth-order galaxy model are efficiently captured by terms in an orthogonal series (analogous to a Fourier series), truncated at some order nmax. Separate sets of orthogonal basis functions are used for the density and potential, so the method is sometimes known as the biorthogonal series expansion method. The two series are related through the Poisson equation, enabling the same set of coefficients to be used for both. This represents a significant improvement in flexibility, allowing axisymmetric, triaxial, lopsided, or distorted density distributions to be built up from an underlying simple model. The original motivation for the introduction of biorthogonal basis sets was the study of the stability of collisionless stellar systems. When decomposed into basis functions, the normal modes of stellar systems correspond to the eigenvectors of matrix equations (e.g. Polyachenko & Shukhman 1981; Fridman & Polyachenko 1984; Weinberg 1989; Palmer 1994; Evans & Read 1998). Stability studies of spherical galaxies often use spherical Bessel functions as the radial basis functions and spherical harmonics as the angular basis functions. These are natural choices as they are the eigenfunctions of the Laplacian in spherical coordinates. However, spherical Bessel functions form a discrete basis only if the radial range is finite, so that the galaxy model has to be truncated or of finite extent. For a semi-infinite radial range, Saha (1991) used a transformation to map the radial range to a finite one, and then used the brute force Gram–Schmidt method to build orthogonal basis functions. Subsequently, basis-function techniques were used to provide algorithms to evolve collisionless stellar systems, partly as a check on the results of linear stability theory. For example, Allen, Palmer & Papaloizou (1990) used a biorthogonal spherical Bessel basis-function expansion to study the radial and circular orbit instability in spherical galaxy models. As instabilities in stellar systems often arise from nearly resonant orbits, so accurate modelling of the precession of individual orbits for several orbital periods is very important. Hernquist & Ostriker (1992) clearly articulated the advantages of basis-function expansions for which the zeroth-order model actually resembles a galaxy. They used the clever method pioneered by Clutton-Brock (1972) to transform the (radial part of the) Laplacian into an equation whose eigenfunctions are already known. They tellingly remarked that ‘for reasons not immediately obvious to us, Clutton-Brock's approach has been virtually ignored in the literature’. In the Hernquist & Ostriker (1992) basis-function expansion, the zeroth-order model is a Hernquist (1990) density profile, which resembles an elliptical galaxy or dark matter halo with ρ ∼ 1/r at small radii and ρ ∼ 1/r4 at large radii. Direct simulation of N-body gravitational systems requires O(N2) operations to calculate the forces on the particles. However, stellar systems are typically collisionless, meaning two-body relaxation effects between bodies are small and the particles chiefly move under the influence of the smooth, mean gravitational field. It is this fact that enables the use of efficient basis-function expansion techniques for calculating the force (see e.g. Dehnen 2014). When forces are calculated from the truncated series, we automatically avoid spurious relaxation effects that may otherwise have been introduced due to undersampling of the number of particles in the system. The coefficients are only calculated once (per time-step), and the series of functions is evaluated once per particle to find the force, so such methods are O(N) in the number of particles. More recently, basis-function techniques have enjoyed a resurgence, inspired by a very different motivation. Lowing et al. (2011) pointed out the value of describing the evolution of direct and computationally expensive simulations of galaxy formation using basis functions. Snapshots of the original simulation with redshift are decomposed with basis functions. The entire simulation can then be replayed many times at will. For example, new objects can be inserted into the simulations and their behaviour studied as if they had been present originally. This technique has already been exploited in studies of the evolution of tidal streams in the Galaxy (Ngan et al. 2015, 2016). However, the gamut of applications is far broader, including the evolution of accreting subhaloes, the dynamics of satellite galaxies and globular clusters in dark haloes, and the response of galactic discs to the changing dark-halo potential. Despite its significant advantages and diverse applications, an important limitation of the basis-function technique is the relative paucity of available sets of ‘analytical’ basis sets – i.e. expressible in a closed form, and preferably obeying some form of recurrence relation so that computation of the basis functions is efficient. Any such basis sets are useful in their own right, but they are particularly interesting if the zeroth-order potential-density pair is close to that of galaxies, allowing for efficient representation using only a few terms. In fact, both Lowing et al. (2011) and Ngan et al. (2015) decomposed numerical halo simulations into the Hernquist & Ostriker (1992) basis-function expansion, as the alternatives are very limited. The purpose of this paper is to broaden the application of the basis-function technique by introducing new sets of biorthogonal basis functions and discussing their applications. We begin by reviewing the presently available expansions and presenting a general approach to the construction of biorthogonal expansions in Section 2 before presenting new specific cases in Section 3. Readers eager to apply the new basis expansion can find the appropriate expressions at the beginning of Section 4. We continue by expanding a series of N-body cosmological haloes that are imitations of the Milky Way dark halo using our new expansion method. We sum up and point out further avenues of development, both on the theoretical and applicable side, in Section 5. 2 THEORY OF BIORTHOGONAL EXPANSIONS 2.1 History The story really begins with the two pioneering papers by Clutton-Brock. Seeking solutions to the cylindrical Poisson equation, Clutton-Brock (1972) introduced the Bessel functions as a set of biorthogonal functions with a continuous eigenvalue on the interval (0, ∞). From these functions, a biorthogonal set with discrete index n was constructed by integrating a particular set of orthogonal functions gn(k) against this continuous eigenvalue (i.e. performing a Hankel transform). This procedure produced a basis set suitable for infinitesimally thin disc galaxies. For Clutton-Brock's choice of gn(k), the lowest order potential-density pair is that of the Kuzmin (1956) or Toomre (1963) disc, though he did not explicitly evaluate the integral and instead left the result as a recurrence relation. Subsequently, Aoki & Iye (1978) reproduced the results of Clutton-Brock (1972), but they explicitly evaluated the Hankel transform to give results written in a closed form. For spherical systems, Clutton-Brock (1973) took a different approach, constructing a biorthogonal basis set for spherical galaxies directly out of the Gegenbauer polynomials together with a particular change of variable. The zeroth-order density profile reproduced the Plummer (1911) model, but this idea was then extended to cover the Hernquist (1990) model by Hernquist & Ostriker (1992) before finally being generalized by Zhao (1996) to a family of basis sets for the hypervirial models (Evans & An 2005). The extension of the Clutton-Brock (1972) approach has only received limited attention in the literature, but we shall demonstrate that it provides a powerful alternative route into constructing further basis sets. Polyachenko & Shukhman (1981) considered the Bessel functions over a finite interval (0, a), where a is the radial extent of the galaxy. This gives a biorthogonal basis set with discrete indices, although the expansion suffers from the issue that the Bessel functions by themselves do not resemble any particular well-known galactic profile. Unlike Clutton-Brock, they did not consider the case of a continuous eigenvalue on the semi-infinite interval (0, ∞). They did however identify the additional degree of freedom, which allows for an adjustment in the asymptotic behaviour of the density basis functions and which we also exploit later. More recently, Rahmati & Jalali (2009) explicitly extended the derivation of Clutton-Brock (1972) and Aoki & Iye (1978) to the spherical case, making an analogous choice of gn(k), and hence deriving a basis set whose lowest order density is the perfect sphere (de Zeeuw 1985). Our work is a generalization of this result, analogous to the way in which Zhao (1996) is a generalization of Clutton-Brock (1973). In fact, Rahmati & Jalali (2009)'s basis set corresponds to the lowest member of our new family of biorthogonal pairs. 2.2 Orthogonal bessel solutions to the Poisson equation We seek solutions in series to the Poisson equation ∇2Φ = 4πGρ (where we work with G = 1 throughout). We represent the potential and density of a self-gravitating system as a sum over basis functions   $$\Phi ({\boldsymbol r}) = \sum _{nlm} \, C_{nlm} \, \Phi _{nlm}({\boldsymbol r}) \qquad \qquad \rho ({\boldsymbol r}) = \sum _{nlm} \, C_{nlm} \, \rho _{nlm}({\boldsymbol r}),$$ (1)where the basis functions are labelled by the tuple (n, l, m). The condition of biorthogonality is written as   $$-\int {\rm d}^{\,3} {\boldsymbol r} \, \Phi _{nlm} \, \rho ^*_{n^\prime l^\prime m^\prime } = N_{nlm} \, \delta _{nlm}^{n^\prime l^\prime m^\prime }, \qquad \qquad \nabla ^2\Phi _{nlm} = 4\pi \rho _{nlm}.$$ (2)We factor out the spherical harmonics Ylm(θ, ϕ) as   $$\Phi _{nlm}({\boldsymbol r}) \equiv \Phi _{nl}(r) Y_{lm}(\theta , \phi ),\qquad \qquad \rho _{nlm}({\boldsymbol r}) \equiv \rho _{nl}(r) Y_{lm}(\theta , \phi ),$$ (3)leaving the equation satisfied by the radial basis functions as   $$\left(\nabla ^2 - \frac{l(l+1)}{r^2}\right)\Phi _{nl} = 4\pi \rho _{nl}.$$ (4) On a finite domain with spherical symmetry, it is natural to choose as a basis the eigenfunctions of the Laplacian, $$f_{klm}({\boldsymbol r}) = j_l(kr)\,Y_{lm}(\theta , \phi )$$, where jl(kr) = r−1/2Jl + 1/2(kr) is a spherical Bessel function. The Bessel function of order μ is defined by the equation   \begin{eqnarray} -\frac{{\rm d}}{{\rm d}z} \left( z^{2\mu +1} \frac{{\rm d}}{{\rm d}z} \left( y_k(z) \right) \right) = k^2 z^{2\mu +1} y_k(z), \ \ \ \ y_k(z) \equiv z^{-\mu } J_\mu (kz), \end{eqnarray} (5)and the eigenfunction fklm obeys the equation ∇2fklm = k2 fklm. However, on an infinite domain, we cannot use these eigenfunctions directly for a series expansion as the radial eigenvalue k becomes continuous: every element of (0, ∞) is a valid eigenvalue k. We therefore convert the basis to a discrete one; first choosing a set of functions gn(k) that are orthogonal on the interval (0, ∞) with respect to the measure k dk and then integrating against fklm (essentially performing a Hankel transform). This produces a biorthogonal basis set, which is suitable for expanding both the potential Φ and the density ρ. Given that we are using different sets of functions to represent the potential and the density, we can exploit an additional freedom first noticed by Polyachenko & Shukhman (1981) and introduce the parameter α in the r-dependence of the Bessel functions. Accordingly, we define the potential and density basis functions as   \begin{eqnarray} \Phi _{nl}(r) \equiv \frac{-A_{nlm}}{\sqrt{r}} \int _0^\infty {\rm d}k \, g_n(k) \, J_{\mu } (kr^{1/(2\alpha )}), \qquad \rho _{nl}(r) \equiv \frac{A_{nlm}}{16\pi \alpha ^2} r^{1/\alpha - 5/2} \int _0^\infty {\rm d}k \, k^2 \, g_n(k) \, J_{\mu } (kr^{1/(2\alpha )}), \end{eqnarray} (6)where μ ≡ α (1 + 2l), and Anlm are constants that will be adjusted to give the desired overall normalization. These functions are biorthogonal in the sense of (2) if the auxiliary functions gn(k) are suitably chosen. Specifically, they must satisfy an orthogonality condition which can be derived thus,   \begin{eqnarray} -\int _0^\infty {\rm d}r \, r^2 \, \Phi _{nl} \, \rho _{n^\prime l} &=& \frac{A_{nlm}\,A_{n^\prime lm}}{8\pi \alpha } \, \int _0^\infty {\rm d}k \, g_n(k) \int _0^\infty {\rm d}q \, q^2 g_{n^\prime }(q) \int _0^\infty {\rm d}z \, z \, J_{\mu }(k z) J_{\mu }(q z) \qquad (z \equiv r^{1/(2\alpha )}) \nonumber\\ &=& \frac{A_{nlm}\,A_{n^\prime lm}}{8\pi \alpha } \, \int _0^\infty {\rm d}k \, k \, g_n(k) \, g_{n^\prime }(k) = \frac{A_{nlm}^2}{8\pi \alpha } \, I_n \, \delta _{n n^\prime }. \end{eqnarray} (7)Here, In is the normalization constant for the functions gn(k) are orthogonal with respect to the weight function k. Let us also write Jlm as the normalization constant for the spherical harmonics Ylm. The overall normalization constant is therefore   $$N_{nlm} = \frac{I_n \, J_{lm} \, A_{nlm}^2}{8\pi \alpha },$$ (8)so that the functions are biorthonormal (i.e. Nnlm = 1) if we set   $$A_{nlm} = \sqrt{\frac{8\pi \alpha }{I_n\,J_{lm}}},$$ (9)in equation (6). The coefficients Cnlm can be calculated by the following integrals over ρ and Φ. When the density is formed from a cloud of point particles as in a numerical simulation, the integral reduces to a sum, which (given this choice of normalization) is written as   \begin{eqnarray} \rho ({\boldsymbol r}) = \sum _i\,m_i \delta ^3({\boldsymbol r} - {\boldsymbol r}_{i}), \qquad C_{nlm} = -\int {\rm d}^{\,3}{\boldsymbol r} \, \Phi _{nlm}({\boldsymbol r})\,\rho ({\boldsymbol r}) = \sum _i m_i\,\Phi _{nlm}({\boldsymbol r}_{i}). \end{eqnarray} (10)One of the advantages of biorthogonal basis sets is that the coefficients in both the density and the potential are the same and can be evaluated very easily from an N-body realization. Note that the use of a discrete sum instead of a continuous integration implies that there is always some numerical (Poisson) noise introduced by discreteness effects. 2.3 Physical boundary conditions The approach presented above yields general sets of biorthogonal functions, but there are further considerations when wishing to construct basis sets for physical systems. Consider the following argument about the required asymptotic behaviour of the radial basis functions. For a given radial order n and angular order l, we define the mass enclosed by a shell with inner radius a and outer radius b as   \begin{eqnarray} M_{nl}(a, b) = 4\pi \int _a^b {\rm d}r \, r^2 \, \rho _{nl} = \int _a^b {\rm d}r \, \left[\frac{{\rm d}}{{\rm d}r}\left(r^2\frac{{\rm d}}{{\rm d}r}\Phi _{nl}\right) - l(l+1)\Phi _{nl}\right]. \end{eqnarray} (11)In order for the density and potential to describe a physical system, we impose the constraint that the mass enclosed by a thin shell goes to zero at the origin and at infinity, and that the potential remains finite everywhere. That is, we require that   \begin{eqnarray} \lim _{a\rightarrow 0} M_{nl}(0,a) = 0, \qquad \lim _{a\rightarrow \infty } M_{nl}(a, \infty ) = 0, \qquad \Phi _{nl}(r) < \infty . \end{eqnarray} (12)These translate into equations that must be satisfied by the derivative of Mnl, which is proportional to ρnl, and so examining (11) we see that the asymptotic solutions for the potential must satisfy: Φnl ∼ rl as r → 0; and Φnl ∼ r−1−l as r → ∞. 3 CONSTRUCTION OF THE NEW BASIS SET 3.1 Zeroth-order models With our general framework in place, we now turn to the job of producing expansions whose zeroth-order model resembles realistic galactic components. This is a lengthy business, so we begin by describing some of the zeroth-order models in our sequence of basis-function expansions. Of course, the zeroth-order models are all spherically symmetric, but the higher order terms in the basis-function expansion describe the effects of flattening or lopsidedness or distortion. Hence, the method can build very general density distributions around the zeroth-order model. The models are labelled by a parameter α; however, those with 0 < α < 1/2 have vanishing density at the origin, and so we consider them no further. Models of interest in galactic dynamics have α ≥ 1/2. The behaviour of the density at small and large radii is   \begin{eqnarray} \lim _{r\rightarrow 0} \rho _{000} (r) \sim r^{-2 + 1/\alpha }, \qquad \lim _{r\rightarrow \infty } \rho _{000} (r) \ \sim r^{-3 -1/(2\alpha )}, \end{eqnarray} (13)so that the α = 1/2 model is cored, whilst the remainder are cusped. The density of all the models falls off with a logarithmic slope between −3 and −4. We proceed to list some of the zeroth-order models that are obtained when α is integer or half-integer, as in these cases the potential reduces to elementary functions; in general, real values of α are permitted, but the potential must be evaluated using the incomplete beta function [see equation (18)]. When α = 1/2, the zeroth-order model is the perfect sphere of de Zeeuw (1985), which has density and potential   \begin{eqnarray} \rho _{000} = \frac{1}{\sqrt{2} \, \pi ^{3/2}} \, \frac{1}{\left(1 + r^2\right)^{2}}, \qquad \Phi _{000} = -\sqrt{2\over \pi } \frac{\arctan {r}}{r}. \end{eqnarray} (14)This is the spherical limit of the perfect ellipsoid, which provides triaxial densities that are close to the luminosity profiles of elliptical galaxies. There is a large literature on these models as the potential is separable or Stäckel, and so the orbits and action-angles can be found as quadratures (e.g. Binney & Tremaine 1987; Bertin 2014). Notice that the density has a harmonic core at the centre but falls like ρ ∼ r−4 at large radii. When α = 1, we obtain the super-NFW model (e.g. Evans & An 2006; An & Zhao 2013; Lilley, Evans & Sanders 2017). This has a central density cusp with ρ ∝ r−1 similar to the Hernquist (1990) or NFW models. Asymptotically, the density falls like ρ ∝ r−3.5 at large radii, which is somewhat faster than the NFW model, but has the distinct advantage of finite total mass. There is evidence to suggest that the NFW may be too shallow at large radii based on the recent work on the splashback radius by Diemer et al. (2017). The potential-density pair of the super-NFW model is   \begin{eqnarray} \rho _{000} = \frac{3}{16\pi } \, \frac{1}{r(1 + r)^{5/2}}, \qquad \Phi _{000} = \frac{-1}{1 + r + \sqrt{1 + r}}. \end{eqnarray} (15)This is ideal for modelling the profiles of dark haloes found in cosmological simulations. This potential-density pair may also be used as a simple model of a galaxy or star cluster. Lilley et al. (2017) show that isotropic and anisotropic distribution functions and solutions of the Jeans equations can be obtained, whilst the projected density is close to de Vaucouleurs and Sérsic profiles. There are of course other finite-mass models which have such desirable properties. A number of models with ρ ∝ r−1 at the centre and with an asymptotic fall off with logarithmic gradient between −3 and −4 are already available in the literature (e.g. An & Evans 2006; An & Zhao 2013). In general, these models are not the zeroth-order terms of any biorthorgonal expansion, so their axisymmetric and triaxial siblings are much more cumbersome to construct. More strongly cusped models can also be obtained as zeroth-order models. When α = 3/2, we have   \begin{eqnarray} \rho _{000} = \frac{2\sqrt{2}}{9 \pi ^{3/2}} \, \frac{1}{r^{4/3}\left(1 + r^{2/3}\right)^{3}}, \qquad \Phi _{000} = \sqrt{2\over \pi }\Bigl [ {1\over r^{2/3}(1 + r^{2/3})} - {{\arctan } (r^{1/3})\over r}\Bigr ]. \end{eqnarray} (16)The central density cusp is now ρ ∼ r−4/3, similar to the inner regions of the massive cosmological clusters studied by Diemand et al. (2005). When α = 2, we have   \begin{eqnarray} \rho _{000} = \frac{15}{64 \pi } \, \frac{1}{r^{3/2}(1 + r^{1/2})^{7/2}}, \qquad \Phi _{000} = -{1 + 2\sqrt{1+\sqrt{r}} \over (1 + \sqrt{1+\sqrt{r}})^2 (1+\sqrt{r})^{3/2}}. \end{eqnarray} (17)This has an inner density profile ρ ∼ r−3/2, similar to the very steepest cusps found in cosmological simulations (Moore et al. 1998). Explicitly, the lowest order basis functions are   \begin{eqnarray} \rho _{000} = {2^{\alpha +1}\Gamma (3/2+\alpha )\over \sqrt{\pi }} {1 \over r^{2 - 1/\alpha }(1+r^{1/\alpha })^{\alpha + 3/2}}\qquad \Phi _{000} = -{2^{\alpha -1} \Gamma (\alpha + 3/2) \over \sqrt{\pi }} {\mathcal {B}_\chi (\alpha ,1/2) \over r}, \qquad \chi \equiv \frac{r^{1/\alpha }}{1 + r^{1/\alpha }}, \end{eqnarray} (18)where $$\mathcal {B}_x(a,b)$$ is the incomplete beta function. The full suite of models has a range of inner density profiles suitable for representing dark haloes with cores and weak or strong cusps, and α can be tuned to match the behaviour of a given halo (see Section 4.3). 3.2 Choice of auxiliary function Now that we have given the spherically symmetric lowest order models, we show how to construct the associated higher order basis functions. We introduce a two-parameter family of auxiliary functions gn(k). We derive an expression for the potential basis functions but find that the requirement of physicality reduces this set to a single-parameter family. This reduction of complexity allows us to express the potential and density basis functions in a more succinct fashion in Sections 3.3 and 3.4. Inspired by Rahmati & Jalali (2009),1 and initially following their calculations, we now choose   \begin{eqnarray} g_n(k) \equiv L^{(\eta )}_n(2k) \, k^{(\eta -1)/2} \, \exp (-k), \end{eqnarray} (19)where $$L^{(\eta )}_n(x)$$ are the associated Laguerre polynomials (Olver et al. 2016, Section 18.3), and η is (for now) a free parameter. The Laguerre polynomials are orthogonal on (0, ∞) with respect to a weight function ω(x) ≡ xη exp (−x),   \begin{eqnarray} \int _0^\infty {\rm d}x \, L^{(\eta )}_n(x) \, L^{(\eta )}_{n^\prime }(x) \omega (x) = \delta _{n n^\prime } {\Gamma (n+\eta +1)\over \Gamma (n+1)}. \end{eqnarray} (20)Our gn(k) consist of a Laguerre polynomial multiplied by a factor of $$\sqrt{\omega (2k)/k}$$, hence ensuring the orthogonality on (0, ∞) with respect to k dk. We can obtain an expression for Φnl using the integral equation 6.621(1) of Gradshteyn & Ryzhik (2014),   \begin{eqnarray} \int _0^\infty {\rm d}u \, \exp (-us) u^\nu J_\mu (u) = \Gamma (\nu + \mu + 1) \, (1 + s^2)^{-(1 + \nu )/2} \, P^{(-\mu )}_\nu \left(\frac{s}{\sqrt{1 + s^2}}\right). \end{eqnarray} (21)Setting u = kr1/(2α) and s = r−1/(2α), and using the explicit polynomial representation of the Laguerre polynomials   $$L^{(\eta )}_n(x) = \sum _{j=0}^n (-1)^j \frac{1}{j!} \left({{n+\eta }\atop n-j}\right) x^j,$$ (22)we find   $$\Phi _{nl} = \sum _{j=0}^n A_{n j\mu \eta } \left[ r \left(1 + r^{1/\alpha }\right)^{(\eta +1)/2 + j}\right]^{-1/2} \, P_{\frac{\eta -1}{2} + j}^{\left(-\mu \right)}\left(\frac{1}{\sqrt{1 + r^{1/\alpha }}}\right), \qquad A_{n j\mu \eta } \equiv \frac{(-2)^j}{j!} \left({{n+\eta }\atop n-j}\right) \Gamma \left(j + \frac{\eta +1}{2} + \mu \right).$$ (23)We use the results of Section 2.3 on asymptotic limits to fix the parameter η. The asymptotic behaviour of the associated Legendre function (Olver et al. 2016, equation 14.8(i)) as z → 1 from above is $$P^{(\mu )}_\nu (z) \sim (1 - z)^{-\mu /2}$$. Therefore, as r → 0, every term in Φnl goes as   $$\Phi _{nl} \sim r^{-1/2} \, \left(1 - \frac{1}{\sqrt{1 + r^{1/\alpha }}}\right)^{\mu /2} \propto r^{\mu /(2\alpha ) - 1/2} = r^l,$$ (24)so the first limit implied by (12) is already satisfied. On the other hand, as r → ∞, $$1/\sqrt{1 + r^{1/\alpha }} \sim r^{-1/(2\alpha )}$$, so the associated Legendre function goes to a constant (Olver et al. 2016, equation 14.5.1) and we are left with the prefactor   $$\Phi _{nl} \sim \left[ r \left(r^{1/\alpha }\right)^{(\eta +1)/2 + j}\right]^{-1/2} \rightarrow r^{-1/2 - (\eta + 1)/(4\alpha )},$$ (25)where we have used the fact that the j = 0 term dominates the sum as r → ∞. By (12), we require this limit to be r−l−1, which we can achieve by setting η = 2μ − 1, so the gn-normalization constant is   $$I_n = \frac{\Gamma (n+2\mu )}{2^{2\mu }\,n!}.$$ (26)Now that η is no longer a free parameter, we can obtain the basis functions in their most convenient form. 3.3 Recurrence relation for the potential basis functions The sum (23) is unsatisfactory for several reasons: (i) it is numerically unstable for high n, (ii) n and j are coupled in such a way that calculating n basis functions requires O(n2) operations, and (iii) associated Legendre functions of non-integer or negative degree or order are rarely implemented numerically. Therefore, departing from Rahmati & Jalali (2009), we seek a superior expression, which can be obtained using the recurrence relation for the Laguerre polynomials,   $$n \, L_n^{(\alpha )}(x) = (n + \alpha ) \, L_{n-1}^{(\alpha )}(x) - x \, L_{n-1}^{(\alpha +1)}(x).$$ (27)Using equation (19), we can immediately write down   $$n \, g_n(k) = (n + 2\mu - 1) \, g_{n-1}(k) - 2 \exp (-k) \, k^\mu \, L_{n-1}^{(2\mu )}(2k).$$ (28)So, we obtain the following recurrence relation for the basis functions   $$n \, \Phi _{nl} = (n + 2\mu - 1) \, \Phi _{n-1,l} + {2A_{nlm}\over \sqrt{r}} \int _0^\infty k^\mu \, L_{n-1}^{(2\mu )}(2k) \, \exp (-k) \, J_\mu (kz) \, {\rm d}k, \qquad (z \equiv r^{1/(2\alpha )}).$$ (29)Evaluating the latter integral requires some work. First, we note that a generating function for the Laguerre polynomials is   \begin{eqnarray} \sum _{n=0}^\infty \, t^n \, L_n^{(\lambda )}(k) = \frac{\exp \left(-tk/(1-t) \right)}{\left(1-t\right)^{\lambda +1}}, \end{eqnarray} (30)so that, using the following Hankel transform (Gradshteyn & Ryzhik 2014, equation 6.623(1))   \begin{eqnarray} \int _0^\infty \, x^\nu \, \mathrm{e}^{-ax} \, J_\nu (xy) \, {\rm d}x = \frac{2^\nu \, \Gamma (\nu + 1/2)}{\sqrt{\pi }} \, \frac{y^\nu }{\left(a^2 + y^2\right)^{\nu + 1/2}}, \end{eqnarray} (31)as well as the generating function for the Gegenbauer polynomials $$C_n^{(\lambda )}(\xi )$$  \begin{eqnarray} \sum _{n=0}^\infty \, t^n \, C_n^{(\lambda )}(\xi ) = \frac{1}{\left(1 - 2\xi t + t^2\right)^\lambda } = \frac{\left(1 + z^2\right)^\lambda }{\left[(1+t)^2 + (1-t)^2z^2\right]^\lambda }, \qquad \left(\xi \equiv \frac{z^2 - 1}{z^2 + 1} \equiv \frac{r^{1/\alpha } - 1}{r^{1/\alpha } + 1}\right), \end{eqnarray} (32)we can find the Hankel transform of the remaining term in the recurrence relation,   \begin{eqnarray} \sum _{n=0}^\infty \, t^n \, \int _0^\infty k^\mu \, L_n^{(2\mu )}(2k) \, \exp (-k) \, J_\mu (kz) \, {\rm d}k & =& \int _0^\infty k^\mu \frac{\exp (-k(1+t)/(1-t))}{\left(1-t\right)^{2\mu + 1}}J_\mu (kz)\,{\rm d}k = \frac{2^\mu \, \Gamma (\mu + 1/2)}{\sqrt{\pi }} \, \frac{z^\mu }{\left[(1+t)^2 + (1-t)^2 z^2\right]^{\mu + 1/2}} \nonumber \\ &=& \sum _{n=0}^\infty t^n \, \frac{2^\mu \, \Gamma (\mu + 1/2)}{\sqrt{\pi }} \, \frac{z^\mu }{\left(1+z^2\right)^{\mu + 1/2}} \, C_n^{(\mu + 1/2)}(\xi ). \end{eqnarray} (33)So the recurrence relation becomes   \begin{eqnarray} n \, \Phi _{nl} = (n + 2\mu - 1) \, \Phi _{n-1,l} + A_{nlm} \frac{2^{\mu +1} \, \Gamma (\mu +1/2)}{\sqrt{\pi }} \, \frac{r^l \, C_{n-1}^{(\mu +1/2)}\left(\xi \right)}{\left(1 + r^{1/\alpha }\right)^{\mu +1/2}}. \end{eqnarray} (34)This means that, in theory, the (n + 1)-th basis function can be trivially calculated from the n-th by simply adding on a single extra term. Evaluating n basis functions therefore requires O(n) steps (although see Section 4.1). A similar approach is taken in Aoki & Iye (1978) in the disc setting. So far, we have only stated without proof the (n = 0, l = 0) case. We now obtain an explicit expression for all the potential basis functions, from which the zeroth order (n = 0, l ≥ 0) can be extracted. We begin by writing our auxiliary function as   \begin{eqnarray} g_n(k) = k^{\mu - 1} \, \exp (-k) \, L_n^{(2\mu - 1)}(2k) = k^{\mu - 1} \, \exp (-k) \left( L_n^{(2\mu - 1)}(0) - 2k\sum _{j=0}^{n-1} \frac{\left({{n+2\mu - 1}\atop n-1-j}\right)}{(n-j)\left({{n}\atop j}\right)} L_j^{(2\mu )}(2k)\right), \end{eqnarray} (35)where again μ = α (1 + 2l). Aside from the constant term, the Hankel transforms can be done using equation (33). To evaluate the Hankel transform of the constant term, we use equation (21) to find [denoting by $$\mathcal {B}_x(a,b)$$ the incomplete beta function]   \begin{eqnarray} \int _0^\infty \, k^{\mu -1} \, \exp (-k) \, J_\mu (kz) \, {\rm d}k = \frac{\Gamma (2\mu )}{\left(1+z^2\right)^{\mu /2}} \, P^{(-\mu )}_{\mu -1}\left(\frac{1}{\sqrt{1+z^2}}\right) = \frac{\Gamma (2\mu )}{\Gamma (\mu )\,2^\mu \,z^\mu } \, \mathcal {B}_\chi (\mu ,1/2), \qquad \left(\chi \equiv \frac{r^{1/\alpha }}{1 + r^{1/\alpha }}\right). \end{eqnarray} (36)So, reassembling the terms, we obtain an expression for the potential as   \begin{eqnarray} \Phi _{nl}(r) = -A_{nlm} \, \frac{2^\mu \, \Gamma (\mu + 1/2) \, \Gamma (n + 2\mu )}{\sqrt{\pi } \, n!} \left[\frac{\mathcal {B}_\chi (\mu ,1/2)}{2\,\Gamma (2\mu )\,r^{1+l}} - \frac{2r^l}{\left(1 + r^{1/\alpha }\right)^{\mu +1/2}}\sum _{j=0}^{n-1}\frac{j! \, C_j^{(\mu +1/2)}(\xi )}{\Gamma (j + 2\mu + 1)}\right], \qquad \left( \chi \equiv \frac{1+\xi }{2}\right), \end{eqnarray} (37)with μ ≡ α (1 + 2l). The only special functions required to evaluate the potential are the incomplete beta function and the Gegenbauer polynomials, which are standard library functions in many numerical software packages (both are included in the GNU Scientific Library, for example). Using the results above, we see that it is easy to write down and compute a basis set for all real values of α. 3.4 The density basis functions A similar calculation can be performed to find expressions for the density basis functions, though it is simpler as the generating function can be applied immediately with no further manipulation of gn(k). First, we note the following Hankel transform (Gradshteyn & Ryzhik 2014, equation 6.623(2))   \begin{equation*} \int _0^\infty x^{\nu +1} \, \exp (-ax) \, J_\nu (xy) \, {\rm d}x = \frac{2^{\nu +1} \, \Gamma (\nu + 3/2)}{\sqrt{\pi }} \, \frac{a \, y^\nu }{(a^2 + y^2)^{\nu + 3/2}}. \end{equation*} Then, using (30) and (32) as above, we have   \begin{eqnarray} \sum _{n=0}^\infty t^n \, \rho _{nl} \propto r^{1/\alpha - 5/2} \, \sum _{n=0}^\infty \, t^n\int _0^\infty k^{\mu + 1} \, L_n^{(2\mu -1)}(2k) \, \mathrm{e}^{-k} \, J_\mu (kz) \, {\rm d}k = \frac{r^{1/\alpha - 5/2}}{\left(1-t\right)^{2\mu }} \, \int _0^\infty k^{\mu +1} \, \mathrm{e}^{-k(1+t)/(1-t)} \, J_\mu (kz) \, {\rm d}k \end{eqnarray} (38)  \begin{eqnarray} \qquad = \frac{2^{\mu +1} \, \Gamma (\mu + 3/2)}{\sqrt{\pi }} \, r^{1/\alpha - 2 + l} \, \frac{(1+t)(1-t)^2}{\left[(1+t)^2 + (1-t)^2z^2\right]^{\mu + 3/2}} \end{eqnarray} (39)  \begin{eqnarray} \qquad = \frac{2^{\mu +1} \, \Gamma (\mu + 3/2)}{\sqrt{\pi }} \, \frac{r^{1/\alpha -2+l}}{(1+r^{1/\alpha })^{3/2+\mu }} \, \sum _{n=0}^\infty t^n \, \left[C_n^{(\mu + 3/2)}(\xi ) - C_{n-1}^{(\mu + 3/2)}(\xi ) - C_{n-2}^{(\mu + 3/2)}(\xi ) + C_{n-3}^{(\mu + 3/2)}(\xi )\right]. \end{eqnarray} (40)The last line can be simplified by adding together two Gegenbauer polynomial recursion relations (Gradshteyn & Ryzhik 2014, equations 8.933(2) and 8.933(3)), resulting in   \begin{eqnarray} \rho _{nl}(r) = \frac{A_{nlm}}{16\pi \alpha ^2} \, \frac{2^{\mu +1}\Gamma (\mu + 1/2)}{\sqrt{\pi }}\frac{r^{1/\alpha -2+l}}{(1+r^{1/\alpha })^{3/2+\mu }}\left[\left(n+\mu +1/2\right)\,C_n^{(\mu +1/2)}(\xi )-\left(n+\mu -1/2\right)\,C_{n-1}^{(\mu + 1/2)}(\xi )\right]. \end{eqnarray} (41)Just as for the potential basis functions, evaluating n density basis functions requires only O(n) steps. It is now straightforward to verify that the basis functions satisfy equation (4). If the spherical harmonics are normalized to Jlm = 4π and we set Anlm = 1, the overall normalization constant for the basis functions written in equations (37) and (41) is   \begin{eqnarray} N_{nl} = \frac{\Gamma (2\mu + n)}{2^{2\mu } \, n! \, 8\pi \alpha }. \end{eqnarray} (42)For numerical purposes, it may be desirable to redefine Anlm to incorporate more of the n- and l-dependent prefactors in equations (37) and (41) (see Section 4.1). 3.5 Comparison with the Zhao and Hernquist-Ostriker basis sets Prior to our work, only a single family of biorthogonal potential-density basis functions was known (although Rahmati & Jalali 2009, provided an alternate single set of solutions with no free parameters which our results now encompass). This basis set was introduced by Zhao (1996). It is written in terms of a parameter α that affects both the inner and outer slopes of the double power law.   \begin{eqnarray} \Phi _{nl}^{\mathrm{Zhao}} = \frac{r^l}{\left(1 + r^{1/\alpha }\right)^{\mu }} \, C^{(\mu + 1/2)}_n\left(\xi \right), \ \ \ \ \rho _{nl}^{\mathrm{Zhao}} = \frac{r^{l - 2 + 1/\alpha }}{\left(1 + r^{1/\alpha }\right)^{\mu + 2}} \, C^{(\mu + 1/2)}_n\left(\xi \right), \end{eqnarray} (43)where μ and ξ are defined as in equation (37). The specific case when α = 1/2 was discovered by Clutton-Brock (1973) and has the Plummer (1911) model as its lowest order term. The case when α = 1 was first proposed by Hernquist & Ostriker (1992) and is the most widely used of all the basis-function expansions, with the Hernquist (1990) model as its lowest order term. We note immediately the close similarity in form to our previous definitions in (37) and (41), where the basis functions are also expressed in terms of the Gegenbauer polynomials. Comparing the lowest order density of the Zhao basis set with that of our new basis set,   \begin{equation*} \rho _{00} \propto \frac{1}{r^{2 - 1/\alpha }\left(1 + r^{1/\alpha }\right)^{\alpha + 3/2}}, \qquad \rho _{00}^{\mathrm{Zhao}} \propto \frac{1}{r^{2 - 1/\alpha} (1 + r^{1/\alpha })^{\alpha + 2}}. \end{equation*} Evidently, the only difference is the shallower outer slope in the former case. This is significant as popular models for dark matter haloes tend to have outer slopes closer to r−3. For example, the generalized NFW profile (Navarro et al. 2004) has ρ ∝ r−γ(1 + r)γ−3, with values for the inner slope γ ranging from 0.7 to 1.5, and the outer slope fixed to r−3. In both our and Zhao's expansion families, when the inner slope is fixed (γ ≡ 2 − 1/α), the asymptotic outer slope is then constrained. To give some examples, if we set γ = 0.7, then Zhao's outer slope is r−4.3 whereas ours is r−3.65; and if γ = 1.5 then Zhao's is r−3.5 and ours is r−3.25. We therefore expect that our basis functions will be more efficient than Zhao's for describing typical dark matter haloes, such as those resimulated by Lowing et al. (2011). Before we pass on to numerical implementation of our new family, we remark that our and Zhao's biorthonormal expansions do not by any means exhaust all the possibilities. For example, we have identified a further family of basis expansions for which the zeroth-order members have density profiles that follow Einasto-like laws (cf. Einasto 1965; Merritt et al. 2006). This will be the subject of a later paper. 4 NUMERICAL PERFORMANCE OF THE BASIS-FUNCTION EXPANSION We turn to the application of our new basis expansion to the representation of cosmological haloes. Before inspecting both analytic and numerical haloes in Sections 4.2 and 4.3, we present a formulation of the basis expansion that is computationally friendly. 4.1 Numerical implementation It is very important to take advantage of a number of recursion relations for purposes of speed. It is also more efficient to factor out the constant parts of the potential-density pairs such that only the spatially dependent pieces are calculated for each particle. We therefore modify the definition of the basis functions as   \begin{eqnarray} {\hat{\rho }}_{nlm}({\boldsymbol r}) = \hat{\rho }_{nl}(r)Y_{lm}(\theta , \phi ) \qquad {\hat{\rho }}_{nl} = \frac{r^{1/\alpha -2+l}}{(1+r^{1/\alpha })^{\mu +3/2}}\left (\left(n+\mu +\frac{1}{2}\right)C_n^{(\mu +1/2)}(\xi )-\left(n+\mu -\frac{1}{2}\right)C_{n-1}^{(\mu +1/2)}(\xi )\right), \end{eqnarray} (44)  \begin{eqnarray} {\hat{\Phi }}_{nlm}({\boldsymbol r}) = \hat{\Phi }_{nl}(r)Y_{lm}(\theta , \phi )\qquad {\hat{\Phi }}_{nl} = \frac{\mathcal {B}_\chi (\mu ,1/2)}{2\,r^{1+l}} - \frac{2 r^l}{(z^2+1)^{\mu + 1/2}} \, \sum _{j=0}^{n-1} \frac{j!\Gamma (2\mu )}{\Gamma (2\mu +j+1)} \, C_j^{(\mu + 1/2)}(\xi ), \end{eqnarray} (45)so that Poisson's equation becomes $$\nabla ^2 {\hat{\Phi }}_{nlm}=4\pi K_{nl} {\hat{\rho }}_{nlm}$$ with   \begin{eqnarray} K_{nl} = -\frac{n!\Gamma (2\mu )}{8\pi \alpha ^2\Gamma (2\mu +n)}, \end{eqnarray} (46)and the normalization constant   \begin{eqnarray} \hat{N}_{nlm} = \int \mathrm{d}^3{\boldsymbol r}\, \hat{\Phi }^{}_{nlm}({\boldsymbol r})\hat{\rho }^*_{nlm}({\boldsymbol r})=J_{lm}\frac{\alpha \sqrt{\pi }\Gamma (\mu )}{2^{1+2\mu }\Gamma (\frac{1}{2}+\mu )}, \qquad J_{lm}=\int \mathrm{d}\phi \, \mathrm{d}\theta \sin \theta \,Y_{lm}(\theta , \phi ) Y_{lm}^*(\theta , \phi ). \end{eqnarray} (47)The Gegenbauer polynomials can be constructed recursively using the relation (Gradshteyn & Ryzhik 2014, equation 8.933(1))   \begin{eqnarray} n C_n^{(w)}(\xi ) = 2(w+n-1)\xi C_{n-1}^{(w)}(\xi )-(2w+n-2)C_{n-2}^{(w)}(\xi ), \end{eqnarray} (48)where $$C_0^{(w)}(\xi )=1$$ and $$C_1^{(w)}(\xi )=2w\xi$$. To construct the potential function $${\hat{\Phi }}_{nl}$$, we define   \begin{eqnarray} A_n=\frac{2n!\Gamma (2\mu )}{\Gamma (2\mu +n+1)}, \end{eqnarray} (49)which satisfies the recurrence relation   \begin{eqnarray} \frac{A_{n+1}}{A_{n}}=\frac{n+1}{n+1+2\mu }, \ \ \ A_0=1/\mu , \end{eqnarray} (50)such that the ladder of potential functions at fixed l is given by   \begin{eqnarray} \hat{\Phi }_{nl}=\hat{\Phi }_{(n-1)l}-\frac{r^l}{(1+z^2)^{\mu +1/2}}A_{n-1}C_{n-1}^{(\mu +1/2)}(\xi ). \end{eqnarray} (51)A naive implementation of this recursion relation, wherein one builds up the higher n terms by starting from Φ0l, results in large errors for high n. In the left-hand panel of Fig. 1, we show the logarithm of the relative difference between the potential computed using this naive approach and an arbitrary precision calculation from Mathematica. The error grows with increasing n and l, as the computation requires taking a small difference between large numbers. To see this, we note that a valid series expansion of the incomplete beta function is [using Olver et al. (2016, 8.17.8) and Fields & Wimp (1961, equation 2.5)]   \begin{eqnarray} \mathcal {B}_\chi (\mu ,1/2) = \frac{\chi ^\mu }{\sqrt{1+z^2}} \sum _{j=0}^\infty A_j C_j^{(\mu +1/2)}(\xi ). \end{eqnarray} (52)A comparison with (51) shows that the terms in $$\hat{\Phi }_{nl}$$ tend to zero as n → ∞, resulting in the aforementioned catastrophic cancellation. This expression, however, provides us with an alternative method of computing $$\hat{\Phi }_{nl}$$. At some large order N, e.g. N = 2nmax, we assume that AN ≈ 0. Then we can write   \begin{eqnarray} \hat{\Phi }_{nl} \approx \frac{r^l}{(1+z^2)^{\mu +1/2}}\sum _{j=n}^{N}A_jC_j^{(\mu +1/2)}(\xi ), \end{eqnarray} (53)which requires us to calculate $$C^{(\mu +1/2)}_{n}(\xi )$$ and An recursively up to order N. We then recursively construct the potential basis functions $$\hat{\Phi }_{nl}$$ downwards using equation (51) where now all $$\hat{\Phi }_{nl}$$ are accurate to the magnitude of AN. This procedure2 results in the reduced errors shown in the right-hand panel of Fig. 1. We only perform this procedure for l > 4 as for l ≤ 4 the naive implementation is satisfactory and the decay of $$\hat{\Phi }_{nl}$$ is weak for small l. Figure 1. View largeDownload slide Relative error in potential basis-function calculation using two methods: on the left, the result of using a forward recursion scheme and on the right using a combination of forward and backward recursion. The relative error is the difference between the potential computed by the basis expansion and an arbitrary precision calculation from Mathematica computed at r/rs = 1.3 using α = 1.2. Note that the floor on the error in the right-hand panel is set by the precision of our Mathematica calculation and can be lowered if so desired. We only perform the ‘forwards+backwards’ procedure for l > 4. Figure 1. View largeDownload slide Relative error in potential basis-function calculation using two methods: on the left, the result of using a forward recursion scheme and on the right using a combination of forward and backward recursion. The relative error is the difference between the potential computed by the basis expansion and an arbitrary precision calculation from Mathematica computed at r/rs = 1.3 using α = 1.2. Note that the floor on the error in the right-hand panel is set by the precision of our Mathematica calculation and can be lowered if so desired. We only perform the ‘forwards+backwards’ procedure for l > 4. Finally, note that the $$\hat{N}_{nl}$$ are independent of n and the Knl satisfy the properties   \begin{eqnarray} \frac{K_{(n+1)l}}{K_{nl}}=\frac{n+1}{n+2\mu }, \ \ \ K_{0l}=-\frac{1}{8\pi \alpha ^2}. \end{eqnarray} (54)With these definitions, a set of coefficients $$\hat{C}_{nlm}$$ is computed from a cloud of particles, and the potential and density reconstructed as   \begin{eqnarray} \hat{C}_{nlm} = \frac{1}{\hat{N}_{nlm}K_{nl}}\sum _i m_i \hat{\Phi }_{nl}(r_i)Y_{lm}(\theta _i, \phi _i), \qquad \Phi = \sum _{nlm}\hat{C}_{nlm}\hat{\Phi }_{nlm}, \qquad \rho = \sum _{nlm}K_{nl}\hat{C}_{nlm}\hat{\rho }_{nlm}. \end{eqnarray} (55) 4.2 Analytical haloes 4.2.1 Spherical NFW models We first consider the reconstruction of spherical NFW haloes using the radial terms of the expansion. We compare the α = 1 member of our family of expansions with the Hernquist & Ostriker (1992, hereafter HO) expansion, which is the α = 1 member of the Zhao (1996) family. Both basis sets have lowest order densities with 1/r cusps as r → 0. In Fig. 2, we show the expansion of a spherical NFW potential and density, using only radial terms (n ≥ 0, l = 0). With an equal number of terms, our α = 1 basis set performs better than the corresponding Hernquist-Ostriker set due to its closer approximation to the NFW profile in the asymptotic fall-off of the lowest order density, ρ ∼ r−7/2. The corresponding behaviour of the Hernquist-Ostriker basis set is ρ ∼ r−4. Note that the r-axis of Fig. 2 is logarithmically scaled and measured in units of the scale length, so we in fact get several hundred additional scale lengths of accurate behaviour. The improved convergence at large radii also reduces the amplitudes of the oscillations at smaller radii. The expansion of the potential is always more accurate than that of the density because the oscillations are integrated over, and therefore effectively smoothed; the accuracy of the radial acceleration expansion lies between that of the potential and the density. According to Fig. 3, the coefficients used in the spherical NFW expansion empirically follow power laws for both basis sets. In our expansion, Cn00 ∼ n−2 whereas in the Hernquist-Ostriker expansion Cn00 ∼ n−3. Figure 2. View largeDownload slide Reconstruction of a spherical NFW potential (left) and density (right) with the new super-NFW basis set (blue) and the Hernquist-Ostriker basis set (red). The distance is given in units of the NFW scale length. Both expansions use radial terms up to n = 20 and no angular terms (l = 0). Both expansions oscillate around the true value, which is represented by the horizontal grey line. Figure 2. View largeDownload slide Reconstruction of a spherical NFW potential (left) and density (right) with the new super-NFW basis set (blue) and the Hernquist-Ostriker basis set (red). The distance is given in units of the NFW scale length. Both expansions use radial terms up to n = 20 and no angular terms (l = 0). Both expansions oscillate around the true value, which is represented by the horizontal grey line. Figure 3. View largeDownload slide The run of the radial expansion coefficients with n for a spherical NFW model using the new basis set (blue) and the Hernquist-Ostriker basis set (red). Also plotted as dashed lines are the n−2 and n−3 curves, suggesting that the coefficients fall off asymptotically like n−2 in our case and like n−3 for the Hernquist-Ostriker case. Figure 3. View largeDownload slide The run of the radial expansion coefficients with n for a spherical NFW model using the new basis set (blue) and the Hernquist-Ostriker basis set (red). Also plotted as dashed lines are the n−2 and n−3 curves, suggesting that the coefficients fall off asymptotically like n−2 in our case and like n−3 for the Hernquist-Ostriker case. To make quantitative statements about the error, we follow Vasiliev (2013) and calculate the integrated squared error. This is the mass-weighted fractional density error defined by   \begin{eqnarray} ISE = \int _{r > r_\mathrm{min} \atop r < r_\mathrm{max}} \frac{\left|\rho _\mathrm{exact} - \rho _\mathrm{approx}\right|^2}{\rho _\mathrm{exact}} \, {\rm d}^{\,3}{\boldsymbol r}. \end{eqnarray} (56)Here, ρapprox is the density reconstructed using basis functions up to a radial order of nmax. Although other measures of error can be constructed, Vasiliev's suggestion is appealing as it is mass-weighted and does not bias the result towards the outer parts of the model. We use rmin = 0.01 and rmax = 100, as this is the range over which we expect the expansions to be most applicable in astrophysical problems. The run of this error measure with nmax is shown on the left in Fig. 4 for the two spherical NFW expansions under consideration. It is clear both that higher values of α give better accuracy, and that for given α and nmax our basis set is more accurate than Zhao's in the task of reproducing NFW haloes. To illustrate this difference, on the right in Fig. 4 we show how the optimum α (the minimum of each plot) varies with differing nmax, with our basis set performing better than Zhao's in each and every case. Figure 4. View largeDownload slide Left: The integrated squared error as defined in equation (56) between the exact NFW profile and the profile reconstructed with radial basis functions up to order nmax for different values of α in Zhao's basis set (red lines) and the new basis set in this paper (blue lines). Right: The integrated squared error between the exact NFW profile and the profile reconstructed with radial basis functions, against the parameter α that characterises each basis set for different numbers of radial basis functions. Again, red lines use Zhao's basis set and blue lines use the new basis set. Figure 4. View largeDownload slide Left: The integrated squared error as defined in equation (56) between the exact NFW profile and the profile reconstructed with radial basis functions up to order nmax for different values of α in Zhao's basis set (red lines) and the new basis set in this paper (blue lines). Right: The integrated squared error between the exact NFW profile and the profile reconstructed with radial basis functions, against the parameter α that characterises each basis set for different numbers of radial basis functions. Again, red lines use Zhao's basis set and blue lines use the new basis set. A detailed error analysis of basis-function techniques is carried out in Kalapotharakos, Efthymiopoulos & Voglis (2008), who claim that a significant factor in obtaining high accuracy is choosing a basis set whose lowest order density function has the correct inner slope behaviour. For a spherical NFW model, this would imply that the α = 1 expansion is best, as the zeroth-order model behaves like ρ ∼ r−1 at small radii. It is possible to choose combinations of the density basis functions so that the most diverging terms at the centre cancel, as HO already showed by building a cored-density model from the n = 0 and n = 1 monopole terms of their basis set. This delicate balance though is susceptible to numerical error. To reproduce accurately the precession of orbits that pass close to the centre, the criterion of Kalapotharakos et al. (2008) seems reasonable. However, it is at odds with the plots in Fig. 4, which show that models with α ≈ 3 or ρ ∼ r−5/3 provide the smallest integrated square error in the density, if nmax ≈ 20. A similar result was obtained by Vasiliev (2013), who used an identical integrated error measure to the present work and concluded that a somewhat higher value of α is preferable for providing a global fit to cosmological haloes, even at the expense of accuracy of the inner slope exponent near the centre. Clearly, the best choice depends on the application in hand. 4.2.2 Flattened NFW models It is important to consider flattened objects, as the haloes found in N-body simulations are generically flattened or triaxial (e.g. Jing & Suto 2002; Allgood et al. 2006). We test the performance by attempting to reconstruct a flattened density profile and including the angular terms or (l, m) terms in the series. An axisymmetric NFW density profile is $$\rho = {\bar{m}}^{-1}(1+ \bar{m})^{-2}$$, where $${\bar{m}}^2 = x^2 + y^2 + z^2/q^2$$. This means the density is stratified on similar concentric spheroids with axis ratio q. Figs 5 and 6 show the expansion of a flattened (q = 0.8) NFW density, and the corresponding acceleration due to the potential. In each case, we compare the α = 1 member of our family of expansions with the Hernquist-Ostriker expansion, which is the α = 1 member of the Zhao (1996) family. Typically, we use nmax = 20 radial basis functions and lmax = 12 angular basis functions. Figure 5. View largeDownload slide Expansion of a flattened (q = 0.8) NFW density. Viewing angle θ is shown on each plot. Both expansions use radial terms up to n = 20 and angular terms up to l = 12. The error measure is Δρ ≡ log |1 − ρapprox/ρexact| (lower is better; the dips are due to the oscillations around the true value). Figure 5. View largeDownload slide Expansion of a flattened (q = 0.8) NFW density. Viewing angle θ is shown on each plot. Both expansions use radial terms up to n = 20 and angular terms up to l = 12. The error measure is Δρ ≡ log |1 − ρapprox/ρexact| (lower is better; the dips are due to the oscillations around the true value). Figure 6. View largeDownload slide Expansion of the radial acceleration in a flattened (q = 0.8) NFW potential. Viewing angle θ is shown on each plot. Both expansions use radial terms up to n = 20 and angular terms up to l = 12. The error measure is $$\Delta _{a_r} \equiv \log {\left|1 - a_{r\,\mathrm{approx}}/a_{r\,\mathrm{exact}}\right|}$$. Figure 6. View largeDownload slide Expansion of the radial acceleration in a flattened (q = 0.8) NFW potential. Viewing angle θ is shown on each plot. Both expansions use radial terms up to n = 20 and angular terms up to l = 12. The error measure is $$\Delta _{a_r} \equiv \log {\left|1 - a_{r\,\mathrm{approx}}/a_{r\,\mathrm{exact}}\right|}$$. Each of the reconstructed quantities is plotted along three polar angles (θ). Note that the convergence is always superior nearer the equator (θ = π/2) than at the poles (θ = 0). This is a feature of any expansion involving spherical harmonics and can be remedied by introducing additional ‘l’ terms. In the flattened case, both expansions lose accuracy in the very inner and outer parts of the haloes. As the l-dependence of both our α = 1 set and the HO set is similar, we do not expect either basis set to be favoured in this regard. However, the superior behaviour in the outskirts of our basis set is maintained. Lowing et al. (2011) used the HO basis set to represent haloes, where on the order of tens of terms are used in both the angular and radial directions. They found errors of < 10 per cent are achieved over a few hundred kiloparsecs; we expect our basis set to provide improved accuracy in this regime. 4.3 Numerical haloes We now analyse a collection of ten Milky Way-like dark matter haloes, extracted from a suite of cosmological N-body zoom-in simulations. These simulations are run with the N-body part of gadget-3 which is similar to gadget-2 (Springel 2005). The zoom-in strategy follows Oñorbe et al. (2014) and all initial conditions are generated with music (Hahn & Abel 2011). Cosmological parameters are taken from the Planck Collaboration XVI et al. (2014) with h = 0.679, Ωb = 0.0481, Ω0 = 0.306, $$\Omega _\Lambda =0.694$$, σ8 = 0.827, and ns = 0.962. In order to select our haloes, we first simulate a 50 h−1 Mpc box with 5123 particles from z = 50 to 0. We use rockstar (Behroozi, Wechsler & Wu 2013) to identify haloes and select Milky Way-like haloes which have virial masses between 7.5 × 1011and2 × 1012 M⊙, no major mergers since z = 1, and no haloes with half the mass of Milky Way analogue's mass within 2 h−1 Mpc. For each selected halo, we select all particles within 10 virial radii and run an intermediate resolution zoom-in whose maximum resolution is 20483, corresponding to a particle mass of 1.8 × 106 M⊙. This intermediate step helps to reduce the contamination from low-resolution particles in our final, high-resolution zoom-in. For the final zoom-in, we take the intermediate resolution simulation and select all particles within 7.5 virial radii. We then run a zoom-in with a maximum resolution of 40963, corresponding to a particle mass of 2.23 × 105 M⊙. All of our high-resolution zoom-ins are uncontaminated within 1 h−1 Mpc of the main halo. The detailed properties of each halo are given in Table 1. Table 1. The properties of each numerical halo in our sample. The haloes are listed in the order that they are displayed in the figures. N is the number of particles, Mv is the virial mass, rv is the virial radius, rs is the scale length, riso is the distance at which the slope is isothermal (so riso = rs/(1/2 + α)α), α parametrizes the inner and outer slopes of the density basis functions, p is the y–x axis ratio, and q is the z–x axis ratio. N/106  Mv/(1012 M⊙)  rv/kpc  rs/kpc  riso/kpc  α  p  q  Figure  10.3  1.88  325  65.2  33.6  1.22  0.860  0.811  7  5.0  0.91  256  37.8  20.5  1.18  0.972  0.816  7  8.7  1.57  306  47.6  25.8  1.18  0.889  0.761  7  7.4  1.36  292  52.8  25.9  1.26  0.800  0.733  7  5.9  1.07  269  58.1  30.7  1.20  0.804  0.795  7  5.5  0.89  254  46.0  23.7  1.22  0.909  0.834  7  4.9  0.89  253  35.7  19.8  1.16  0.807  0.780  8  11.3  1.62  310  123  38.1  1.59  0.858  0.769  8  7.7  1.11  273  98.5  32.9  1.54  0.926  0.779  8  7.6  1.71  315  132.6  47.6  1.49  0.861  0.730  10  N/106  Mv/(1012 M⊙)  rv/kpc  rs/kpc  riso/kpc  α  p  q  Figure  10.3  1.88  325  65.2  33.6  1.22  0.860  0.811  7  5.0  0.91  256  37.8  20.5  1.18  0.972  0.816  7  8.7  1.57  306  47.6  25.8  1.18  0.889  0.761  7  7.4  1.36  292  52.8  25.9  1.26  0.800  0.733  7  5.9  1.07  269  58.1  30.7  1.20  0.804  0.795  7  5.5  0.89  254  46.0  23.7  1.22  0.909  0.834  7  4.9  0.89  253  35.7  19.8  1.16  0.807  0.780  8  11.3  1.62  310  123  38.1  1.59  0.858  0.769  8  7.7  1.11  273  98.5  32.9  1.54  0.926  0.779  8  7.6  1.71  315  132.6  47.6  1.49  0.861  0.730  10  View Large For each halo, we take all of the particles within 500 kpc of the main halo in z = 0 snapshot. This corresponds to between 5–12 million particles, depending on the halo mass. We wish to investigate the ability of our new basis sets to represent these numerical haloes. To this end, we must first choose the two global parameters that specify the basis set – the scale length rs and the parameter α. We need not worry at this stage about the overall normalization (related to the total mass) because this will be set automatically when performing the sum over particles via equation (10). To set rs and α, we first fit the zeroth-order density function ρ000, which is a spherically symmetric model. Because we know the virial radius rv for each halo, we perform the fitting procedure in terms of the dimensionless concentration c ≡ rv/rs, as is standard in the literature. The particles are binned logarithmically in radius, and a non-linear least-squares algorithm then adjusts α and c to minimize the difference between the logarithms of the inferred bin density and the model density. For this fitting procedure, we use a form of the zeroth-order density function that is parametrized by the mass enclosed by the virial radius, i.e. which satisfies $$\int _0^{r_{\rm v}} 4\pi \rho (r) r^2 {\rm d}r = M_{\rm v}$$. This is   \begin{eqnarray} \rho (r) = \left[\mathcal {B}_{\chi _v}(\alpha ,1/2) - \frac{c/\alpha }{\left(1 + c^{1/\alpha }\right)^{\alpha +1/2}}\right]^{-1} \, \frac{M_{\rm v} \, (\alpha + 1/2)}{4 \pi \alpha ^2 \, (r_{\rm v}/c)^3} \, \frac{(cr/r_{\rm v})^{1/\alpha - 2}}{\left(1 + (cr/r_{\rm v})^{1/\alpha }\right)^{\alpha + 3/2}}, \end{eqnarray} (57)where χv ≡ c1/α/(1 + c1/α). We can now perform a full expansion, using for each halo the basis set with the determined best values of α and rs. To compare the accuracy in the reproduction of the density, we draw contour plots in each principal plane, as shown in Figs 7, 8, and 10. The smooth underlying density distribution in each principal plane of the original numerical halo is estimated by taking particles in a slab of width 2 kpc around the plane, then creating a two-dimensional histogram with each bin having an area of approximately 3 kpc2. The density as reconstructed from the basis-function expansion is sampled along rays in each principal plane, and then interpolated on to a grid. The same 12 contours are drawn on each plot, spaced approximately logarithmically between 106 and 2 × 108  M⊙ kpc−3. Figure 7. View largeDownload slide From left to right: Density contours in the three principal planes of the original numerical halo, then density contours in the same three principal planes of the halo as reconstructed with basis functions up to n = 20 and l = 12. Distances are in kpc, and the densities are in  M⊙kpc−3. Figure 7. View largeDownload slide From left to right: Density contours in the three principal planes of the original numerical halo, then density contours in the same three principal planes of the halo as reconstructed with basis functions up to n = 20 and l = 12. Distances are in kpc, and the densities are in  M⊙kpc−3. Figure 8. View largeDownload slide Continuation of Fig. 7. The latter two haloes possess a large number of small satellites, so for these we display on the right two reconstructions: top, a reconstruction of the full halo; and bottom, one where the particles energetically bound to the 50 most massive satellites have been removed. This illustrates the issues surrounding expansion accuracy in the presence of a sub-structure. Figure 8. View largeDownload slide Continuation of Fig. 7. The latter two haloes possess a large number of small satellites, so for these we display on the right two reconstructions: top, a reconstruction of the full halo; and bottom, one where the particles energetically bound to the 50 most massive satellites have been removed. This illustrates the issues surrounding expansion accuracy in the presence of a sub-structure. The haloes displayed in Figs 7 and 8 are rotated such that the principal planes are aligned with the coordinate axes (shortest axis along the z-axis), so that the angular expansion coefficients can be compared meaningfully. Distributions of the expansion coefficients for these nine numerical haloes are shown in Fig. 9. One could produce an artificial ‘halo’ with geometry typical for this family of nine numerical haloes by drawing coefficients from the distributions shown in these plots. From Figs 7 and 8, we see that key features of the haloes are resolved correctly, including: (i) the orientation of the principal axes, (ii) the axis ratios in the three principal planes, and (iii) the run of ellipticity with radius. The size of the smallest resolvable feature is limited by the distance between the roots of the polynomial used to define the highest order function used in the reconstruction. Figure 9. View largeDownload slide We denote normalized coefficients by tildes ($$\tilde{C}_{nlm} \equiv C_{nlm}/C_{000}$$) and the standard deviation of a variable x by σ(x). Panels 1 and 2: the radial coefficients $$\tilde{C}_{n00}$$ for nine numerical haloes. Panel 3: the angular coefficients $$\tilde{C}_{0l0}$$. Panel 4: the angular coefficients $$\tilde{C}_{0ll}$$. The latter three plots are normalized by each coefficient's standard deviation to make details easier to see. The trend in the C0l0 coefficients indicates the flattening of the haloes along the z-axis, and the trend in the $$\tilde{C}_{0ll}$$ coefficients shows the elongation along the x-axis. Figure 9. View largeDownload slide We denote normalized coefficients by tildes ($$\tilde{C}_{nlm} \equiv C_{nlm}/C_{000}$$) and the standard deviation of a variable x by σ(x). Panels 1 and 2: the radial coefficients $$\tilde{C}_{n00}$$ for nine numerical haloes. Panel 3: the angular coefficients $$\tilde{C}_{0l0}$$. Panel 4: the angular coefficients $$\tilde{C}_{0ll}$$. The latter three plots are normalized by each coefficient's standard deviation to make details easier to see. The trend in the C0l0 coefficients indicates the flattening of the haloes along the z-axis, and the trend in the $$\tilde{C}_{0ll}$$ coefficients shows the elongation along the x-axis. Two haloes in Fig. 8 demonstrate one pitfall of the method. The presence of unresolved but massive subhaloes can cause a blow-up in the higher order coefficients of the expansion. This would be cancelled by yet higher order terms, but as the expansion is truncated these are not present. This problem is generic – it affects other basis sets, such as Zhao's, to a greater or lesser degree – but it can be remedied by removing the unresolved subhaloes by an automated halo-finding procedure, as demonstrated in the figure. When these subhaloes are removed, the ‘optimal’ values of α and rs (as determined by fitting equation (57)) are reduced to values more characteristic of the other haloes in the sample. The final halo, displayed in Fig. 10, provides a serious challenge. It is accompanied by a massive close-in satellite with a mass of 1.9 × 1011 M⊙, and thus has a highly aspherical geometry. We therefore examine this halo in greater detail in Fig. 10 using nmax = 20 and lmax = 12 (middle panel) and nmax = 40 and lmax = 40 (lower panel). Remarkably, the overall structure of both the halo and the large satellite is correctly resolved even with nmax = 20. This is impressive, as we might have suspected at the outset that two basis-function expansions, centred on each object, would be necessary to reproduce the merging structure. Figure 10. View largeDownload slide A halo with a prominent, massive satellite. From top to bottom: the original halo; the reconstruction with maximum order n = 20 and l = 12; the same with n = 40 and l = 40. Figure 10. View largeDownload slide A halo with a prominent, massive satellite. From top to bottom: the original halo; the reconstruction with maximum order n = 20 and l = 12; the same with n = 40 and l = 40. 5 CONCLUSIONS Biorthogonal density and potential basis functions provide useful and flexible ways of describing realistic dark matter haloes and galaxies, which may be aspherical, triaxial, or further misshapen. The coefficients of these basis-function expansions can be found easily by summing over the particles in an N-body realization and used to reconstruct both the potential and density. We have discovered a completely new family of biorthogonal potential-density pairs, parametrized in terms of α. The zeroth-order model has a density ρ ∼ r−2+1/α at small radii and ρ ∼ r−3−1/(2α) at large radii. This double-power-law profile has a central logarithmic slope between 0 and −2 and an asymptotic logarithmic slope between −3 and −4, making it perfect for representing dark haloes. The zeroth-order model with α = 1/2 has a harmonic core and is the celebrated perfect sphere of de Zeeuw (1985), whilst the model with α = 1 has a 1/r central density cusp and is the super-NFW model (Lilley et al. 2017). For each of these zeroth-order models, we provide a biorthogonal basis-function expansion in terms of standard special functions readily available in numerical libraries. This extends the zeroth-order model into the highly realistic regime of flattened, distorted, and triaxial density profiles. Previously, the only known family of biorthogonal potential-density pairs was the one outlined by Zhao (1996), of which the most widely used member is the HO expansion. The zeroth-order model has a density ρ ∼ r−2+1/α at small radii and ρ ∼ r−3−1/α at large radii. Although the central density cusp is the same, the outer density falls off rather more quickly than in our expansion, making Zhao's family less well-matched to modelling dark haloes. We have demonstrated the capabilities of our basis-function expansions by using them to recreate spherical and flattened dark haloes with analytic densities of Navarro et al. (1997) form. In particular, we showed that our method represents a significant improvement over the HO expansion, giving a more accurate reproduction of the density, potential, and radial force of NFW-like models. Additionally, we decomposed 10 simulated cosmological haloes using our basis functions, computing the coefficients as simple sums over the particles. This yielded very encouraging results with highly flattened and triaxial dark-halo density distributions well reproduced with typically 20 radial and 12 angular terms (giving a total of 20 × 122 ≈ 3000 terms). Simulated dark haloes can be lopsided, distorted or twisted, especially if there is a nearby companion exerting strong tidal forces (cf. the Milky Way and the Magellanic Clouds). Particular striking is the ability of our basis-function expansion to reproduce the density of a highly asymmetric dark halo which is in the process of accreting a companion large subhalo. Lowing et al. (2011) used the HO expansion to create potentials approximating the structure of the dark haloes in the Aquarius simulations. They found that the radial force from the expansion sometimes differed from that in the simulations when the central logarithmic slope of the numerical halo density was shallower than the slope of −1 of the zeroth-order model. They also highlighted a new application for biorthogonal basis-function expansions as a tool for data compression of large and expensive N-body simulations. By calculating sets of coefficients from numerical snapshots at different redshifts and interpolating between them, they obtained a compact and efficient summary of the halo's evolution throughout time. This is very powerful because new objects can then be inserted and the entire simulation can be re-run cheaply making the method well-suited to study a range of problems, such as the dynamics of small satellites, the stripping of globular clusters or the shape of tidal streams. Another important area of application is to the fitting of data on the Milky Way galaxy, which has increased substantially in both quality and quantity over the last few years. Models of the Milky Way galaxy are assembled from simple building blocks. Results of calculations using such ‘pieces of Lego’ are often very troubling. For example, when the position and velocity data of stars in the Sagittarius stream are fitted to such models, the conclusion is that the dark halo is triaxial with the short and long axes in the Galactic plane (Law & Majewski 2010). This configuration is unstable (see e.g. Debattista et al. 2013) and in conflict with observational data on the disc (Kuijken & Tremaine 1994). The strong suspicion is that the inflexible model of the Milky Way's potential prevented proper exploration of parameter space and artificially confined the solution for the dark matter distribution into an unrealistic straitjacket. A completely new way of representing the Galactic dark halo is needed with the advent of large-scale photometric, spectroscopic and astrometric surveys of the Galaxy. For model fitting, the dark-halo potential should be represented by distributions of the coefficients that can be used as priors in Bayesian inference from the data, rather than say a single number (the flattening) in a predetermined and unadaptable density law. Of course, the shape of a dark halo depends on the nature of the dark matter particle and on the extent of feedback processes (see e.g. Sellwood 2004; Macciò et al. 2012). The shape also depends on a host of other factors, including the mass of the halo, its environment (isolated versus group), its recent history (e.g. late infall of a large subhalo), and the presence or absence of a disc. It will be interesting to test our new basis-function expansion method on the full variety of numerically constructed haloes, and understand how the distribution of coefficients changes with the underlying physics. The main impediment to efficient exploitation of these ideas is that so few biorthogonal pairs are known. Our new discovery helps in this regard, but it has not exhausted the supply of such expansions. We will show in a later paper how to extend the methods in this paper to other cosmologically inspired dark-halo density laws. Acknowledgements JLS and EJL acknowledge the support of the STFC. We thank members of the Institute of Astronomy Streams group for discussions and comments as this work was in progress. The anonymous referee is thanked for a number of helpful suggestions. Footnotes 1 Who were in turn motivated by the analogous disc case considered in Clutton-Brock (1972). 2 It is analogous to Miller's method in numerical analysis, apparently introduced by J.C.P. Miller for the computation of Bessel functions (see e.g. Abramowitz & Stegun 1965). REFERENCES Abramowitz M., Stegun I. A., 1965, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables . Dover Publ., New York Allen A. J., Palmer P. L., Papaloizou J., 1990, MNRAS , 242, 576 https://doi.org/10.1093/mnras/242.4.576 CrossRef Search ADS   Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS , 367, 1781 https://doi.org/10.1111/j.1365-2966.2006.10094.x CrossRef Search ADS   An J. H., Evans N. W., 2006, AJ , 131, 782 https://doi.org/10.1086/499305 CrossRef Search ADS   An J., Zhao H., 2013, MNRAS , 428, 2805 https://doi.org/10.1093/mnras/sts175 CrossRef Search ADS   Aoki S., Iye M., 1978, PASJ , 30, 519 Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ , 762, 109 https://doi.org/10.1088/0004-637X/762/2/109 CrossRef Search ADS   Bertin G., 2014, Dynamics of Galaxies . Cambridge Univ. Press, Cambridge Google Scholar CrossRef Search ADS   Binney J., Tremaine S., 1987, Galactic Dynamics . Princeton Univ. Press, Princeton, NJ, p. 747 Clutton-Brock M., 1972, Ap&SS , 16, 101 https://doi.org/10.1007/BF00643095 CrossRef Search ADS   Clutton-Brock M., 1973, Ap&SS , 23, 55 https://doi.org/10.1007/BF00647652 CrossRef Search ADS   de Zeeuw T., 1985, MNRAS , 216, 273 https://doi.org/10.1093/mnras/216.2.273 CrossRef Search ADS   Debattista V. P., Roškar R., Valluri M., Quinn T., Moore B., Wadsley J., 2013, MNRAS , 434, 2971 https://doi.org/10.1093/mnras/stt1217 CrossRef Search ADS   Dehnen W., 2014, Comput. Astrophys. Cosmol. , 1, 1 https://doi.org/10.1186/s40668-014-0001-7 CrossRef Search ADS   Diemand J., Zemp M., Moore B., Stadel J., Carollo C. M., 2005, MNRAS , 364, 665 https://doi.org/10.1111/j.1365-2966.2005.09601.x CrossRef Search ADS   Diemer B., Mansfield P., Kravtsov A. V., More S., 2017, ApJ , 843, 140 https://doi.org/10.3847/1538-4357/aa79ab CrossRef Search ADS   Einasto J., 1965, Tr. Astrofizicheskogo Inst. Alma-Ata , 5, 87 Evans N. W., 1993, MNRAS , 260, 191 https://doi.org/10.1093/mnras/260.1.191 CrossRef Search ADS   Evans N. W., An J. H., 2005, MNRAS , 360, 492 https://doi.org/10.1111/j.1365-2966.2005.09078.x CrossRef Search ADS   Evans N. W., An J. H., 2006, Phys. Rev. D , 73, 023524 https://doi.org/10.1103/PhysRevD.73.023524 CrossRef Search ADS   Evans N. W., Read J. C. A., 1998, MNRAS , 300, 106 https://doi.org/10.1046/j.1365-8711.1998.01864.x CrossRef Search ADS   Fields J. L., Wimp J., 1961, Math. Comput. , 15, 390 https://doi.org/10.1090/S0025-5718-1961-0125992-3 CrossRef Search ADS   Fridman A. M., Polyachenko V. L., 1984, Physics of Gravitating Systems I: Equilibrium and Stability . Springer Verlag, Berlin. Gradshteyn I., Ryzhik I., 2014, Table of Integrals, Series and Products , 8th edn. Academic Press, New York Hahn O., Abel T., 2011, MNRAS , 415, 2101 https://doi.org/10.1111/j.1365-2966.2011.18820.x CrossRef Search ADS   Hernquist L., 1990, ApJ , 356, 359 https://doi.org/10.1086/168845 CrossRef Search ADS   Hernquist L., Ostriker J. P., 1992, ApJ , 386, 375 (HO) https://doi.org/10.1086/171025 CrossRef Search ADS   Jing Y. P., Suto Y., 2002, ApJ , 574, 538 https://doi.org/10.1086/341065 CrossRef Search ADS   Kalapotharakos C., Efthymiopoulos C., Voglis N., 2008, MNRAS , 383, 971 https://doi.org/10.1111/j.1365-2966.2007.12417.x CrossRef Search ADS   Kuijken K., Tremaine S., 1994, ApJ , 421, 178 https://doi.org/10.1086/173635 CrossRef Search ADS   Kuzmin G. G., 1956, Publ. Tartu Astrofizica Obs. , 33, 75 Law D. R., Majewski S. R., 2010, ApJ , 714, 229 https://doi.org/10.1088/0004-637X/714/1/229 CrossRef Search ADS   Lilley E. J., Evans N. W., Sanders J. L., 2017, MNRAS , submitted Lowing B., Jenkins A., Eke V., Frenk C., 2011, MNRAS , 416, 2697 https://doi.org/10.1111/j.1365-2966.2011.19222.x CrossRef Search ADS   Macciò A. V., Paduroiu S., Anderhalden D., Schneider A., Moore B., 2012, MNRAS , 424, 1105 https://doi.org/10.1111/j.1365-2966.2012.21284.x CrossRef Search ADS   Merritt D., Graham A. W., Moore B., Diemand J., Terzić B., 2006, AJ , 132, 2685 https://doi.org/10.1086/508988 CrossRef Search ADS   Moore B., Governato F., Quinn T., Stadel J., Lake G., 1998, ApJ , 499, L5 https://doi.org/10.1086/311333 CrossRef Search ADS   Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ , 490, 493 https://doi.org/10.1086/304888 CrossRef Search ADS   Navarro J. F. et al.  , 2004, MNRAS , 349, 1039 https://doi.org/10.1111/j.1365-2966.2004.07586.x CrossRef Search ADS   Ngan W., Bozek B., Carlberg R. G., Wyse R. F. G., Szalay A. S., Madau P., 2015, ApJ , 803, 75 https://doi.org/10.1088/0004-637X/803/2/75 CrossRef Search ADS   Ngan W., Carlberg R. G., Bozek B., Wyse R. F. G., Szalay A. S., Madau P., 2016, ApJ , 818, 194 https://doi.org/10.3847/0004-637X/818/2/194 CrossRef Search ADS   Olver F. W. J., Olde Daalhuis A. B., Lozier D. W., Schneider B. I., Boisvert R. F., Clark C. W., Miller B. R., Saunders B. V., 2016, NIST Digital Library of Mathematical Functions . Available at: http://dlmf.nist.gov/ Oñorbe J., Garrison-Kimmel S., Maller A. H., Bullock J. S., Rocha M., Hahn O., 2014, MNRAS , 437, 1894 https://doi.org/10.1093/mnras/stt2020 CrossRef Search ADS   Palmer P. L., 1994, Astrophysics and Space Science Library, Vol. 185, Stability of Collisionless Stellar Systems: Mechanisms for the Dynamical Structure of Galaxies . Springer, Berlin Planck Collaboration XVI et al.  , 2014, A&A , 571, A16 https://doi.org/10.1051/0004-6361/201321591 CrossRef Search ADS   Plummer H. C., 1911, MNRAS , 71, 460 https://doi.org/10.1093/mnras/71.5.460 CrossRef Search ADS   Polyachenko V. L., Shukhman I. G., 1981, Sov. Astron. , 25, 533 Rahmati A., Jalali M. A., 2009, MNRAS , 393, 1459 https://doi.org/10.1111/j.1365-2966.2008.14226.x CrossRef Search ADS   Saha P., 1991, MNRAS , 248, 494 https://doi.org/10.1093/mnras/248.3.494 CrossRef Search ADS   Sellwood J. A., 2004, in Ryder S., Pisano D., Walker M., Freeman K., eds, IAU Symp. 220, Dark Matter in Galaxies . Kluwer, Dordrecht, p. 27 Springel V., 2005, MNRAS , 364, 1105 https://doi.org/10.1111/j.1365-2966.2005.09655.x CrossRef Search ADS   Toomre A., 1963, ApJ , 138, 385 https://doi.org/10.1086/147653 CrossRef Search ADS   Vasiliev E., 2013, MNRAS , 434, 3174 https://doi.org/10.1093/mnras/stt1235 CrossRef Search ADS   Weinberg M. D., 1989, MNRAS , 239, 549 https://doi.org/10.1093/mnras/239.2.549 CrossRef Search ADS   Zhao H., 1996, MNRAS , 278, 488 https://doi.org/10.1093/mnras/278.2.488 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations