How gravitational softening affects galaxy stability I. Linear mode analysis of disc galaxies

TL;DR

This study examines how different gravitational softening kernels influence the accuracy of linear mode analysis in disc galaxy simulations, highlighting the importance of minimizing gravity bias for precise eigenmode frequency estimation.

Contribution

It introduces and tests softening kernels with reduced gravity bias, demonstrating improved accuracy in galaxy stability mode analysis compared to traditional kernels like Plummer.

Findings

Softening kernels with minimal gravity bias yield more accurate mode frequencies.

Plummer softening performs poorly in reproducing unsoftened mode properties.

Super-Newtonian forces near the softening length improve mode analysis accuracy.

Abstract

Linear perturbation is used to investigate the effect of gravitational softening on the retrieved two-armed spiral eigenmodes of razor-thin stellar discs. We explore four softening kernels with different degrees of gravity bias, and with/without compact support (compact in the sense that they yield exactly Newtonian forces outside the softening kernel). These kernels are applied to two disc galaxy models with well-known unsoftened unstable modes. We illustrate quantitatively the importance of a vanishing linear gravity bias to yield accurate frequency estimates of the unstable modes. As such, Plummer softening, while very popular amongst simulators, performs poorly in our tests. The best results, with excellent agreement between the softened and unsoftened mode properties, are obtained with softening kernels that have a reduced gravity bias, obtained by compensating for the sub-Newtonian forces at small interparticle distances with slightly super-Newtonian forces at radii near the softening length. We present examples of such kernels that, moreover, are analytically simple and computationally cheap. Finally, these results light the way to the construction of softening methods with even smaller gravity bias, although at the price of increasingly complex kernels.

Tables

Table 1. Table 1: Characteristics of the various softening methods used in this study. Here, ϕ italic-ϕ \phi , ϱ italic-ϱ \varrho , and σ 𝜎 \sigma take argument u = r / ε 𝑢 𝑟 𝜀 u=r/\varepsilon and are the dimensionless kernels for, respectively, potential, spatial and surface density, defined in equations ( 4 ) and ( 12 ). The coefficients a 0 subscript 𝑎 0 a_{0} and a 2 subscript 𝑎 2 a_{2} determine the gravitational biases (see equation  8 ). Here, ε 0 subscript 𝜀 0 \varepsilon_{0} quantifies the scaling of the kernels with non-zero a 0 subscript 𝑎 0 a_{0} to a common level of gravity bias and ε F subscript 𝜀 𝐹 \varepsilon_{F} quantifies the scaling of the kernels to a common maximum inter-particle force. A dagger indicates softening methods with compact support, where the formulæ for ϕ italic-ϕ \phi and ϱ italic-ϱ \varrho only apply at r < ε 𝑟 𝜀 r<\varepsilon or t ≡ 1 − u 2 ∈ [ 0 , 1 ] 𝑡 1 superscript 𝑢 2 0 1 t\equiv 1-u^{2}\in[0,1] . For these methods, no sensible razor-thin surface-density kernel σ ​ ( u ) 𝜎 𝑢 \sigma(u) can be provided.

name		$\varphi \left(u\right)$	$\varrho \left(u\right)$	$\sigma \left(u\right)$	$a_0$	$a_2$	${\epsilon }_0/\epsilon$	${\epsilon }_F/\epsilon$
	Newton	$\frac{1}{u}$	${\delta }_{3D}\left(u\right)$	${\delta }_{2D}\left(u\right)$	0	0	0	0
P0	Plummer	$\frac{1}{\sqrt{1+u^2}}$	$\frac{3}{4\pi }\frac{1}{{\left(1+u^2\right)}^{5/2}}$	$\frac{1}{2\pi }\frac{1}{{\left(1+u^2\right)}^{3/2}}$	$2\pi$	$\infty$	1	1
Q2	$\begin{array}{c}\text{2}\text{<sup>nd</sup>}\text{ modified}\\ \text{Kuz’min}\end{array}$	$\frac{3+\frac{5}{2}{u}^2+u^4}{{\left(1+u^2\right)}^{5/2}}$	$\frac{15}{8\pi }\frac{4-3u^2}{{\left(1+u^2\right)}^{9/2}}$	$\frac{3}{4\pi }\frac{4-u^2}{{\left(1+u^2\right)}^{7/2}}$	0	$-\frac{\pi }{3}$	—	2.568
†F3	Ferrers $n=3$	$1+\frac{1}{2}t+\frac{3}{8}{t}^2+\frac{5}{16}{t}^3+\frac{35}{128}{t}^4$	$\frac{315}{64\pi }{t}^3$	ill-behaved	$\frac{63\pi }{128}$	$\frac{7\pi }{1024}$	$\frac{63}{256}$	2.309
†L2	$\begin{array}{c}\text{2D modified}\\ \text{Ferrers }n=2\end{array}$	$1+\frac{1}{2}t+\frac{3}{8}{t}^2+\frac{5}{2}{t}^3$	$\frac{105}{8\pi }\left(1-2u^2\right)t$	ill-behaved	$0$	$-\frac{7\pi }{384}$	—	3.711

Equations

$$\frac{\partial F(\mathbf{x},\mathbf{v},t)}{\partial t}+\mathbf{v}\cdot \frac{\partial F(\mathbf{x},\mathbf{v},t)}{\partial \mathbf{x}}+\frac{\partial V(\mathbf{x},t)}{\partial \mathbf{x}}\cdot \frac{\partial F(\mathbf{x},\mathbf{v},t)}{\partial \mathbf{v}}=0.$$

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}$$

$$V(\mathbf{x},t)$$

$$\approx G\sum _i\frac{m_i}{|\mathbf{x}-{\mathbf{x}}_i|}$$

$$\frac{1}{r}\to \psi (r)=\frac{1}{\epsilon }\varphi \left(\frac{r}{\epsilon }\right),$$

$$V(\mathbf{x})\approx \hat{V}(\mathbf{x})=\sum _i\frac{G{m}_i}{\epsilon }\varphi \left(\frac{|\mathbf{x}-{\mathbf{x}}_i|}{\epsilon }\right).$$

$${\mathrm{var}}_{\mathbf{x}}(\hat{V})=⟨{\left[⟨\hat{V}(\mathbf{x})⟩-\hat{V}(\mathbf{x})\right]}^2⟩,$$

$${\mathrm{bias}}_{\mathbf{x}}(\hat{V})=⟨\hat{V}(\mathbf{x})⟩-V(\mathbf{x}).$$

$${\mathrm{bias}}_{\mathbf{x}}(\hat{V})=a_0\epsilon G\Sigma (\mathbf{x})+a_2{\epsilon }^{3}G{\nabla }^2\Sigma (\mathbf{x})+\mathcal{O}({\epsilon }^5),$$

$$a_n=\frac{2\pi }{2^n([n/2]!{)}^2}{\int }_0^{\infty }\left[1-u\varphi (u)\right]u^n\mathrm{d}u.$$

$${\mathrm{bias}}_{\mathbf{x}}(\hat{V})$$

$${\mathrm{bias}}_{\mathbf{x}}(\hat{\mathbf{a}})$$

$$\rho (r)=\frac{m}{{\epsilon }^3}\varrho \left(\frac{r}{\epsilon }\right)\text{and}\Sigma (r)=\frac{m}{{\epsilon }^2}\sigma \left(\frac{r}{\epsilon }\right),$$

$$\varrho (u)$$

$$\sigma (u)$$

$$\varphi (u)=\frac{K_n(u^2)}{{\left(1+u^2\right)}^{n+1/2}},$$

$$a_0$$

$$=\underset{x\to \infty }{\lim }\left\{(1-c_1)x+2c_1-c_0+\mathcal{O}\left(\frac{1}{x}\right)\right\}$$

$$\varphi (u)=\frac{2+u^2}{(1+u^2{)}^{3/2}}$$

$$\varrho (u)=\frac{3}{4\pi }\frac{4-u^2}{(1+u^2{)}^{7/2}}\text{ and }\sigma (u)=\frac{3}{2\pi }\frac{1}{(1+u^2{)}^{5/2}}.$$

$$a_0$$

$$=\underset{x\to \infty }{\lim }\left\{(1-c_2)x+\frac{1}{3}(8c_2-2c_1-c_0)+\mathcal{O}\left(\frac{1}{x}\right)\right\}$$

$$a_2$$

$$=\underset{x\to \infty }{\lim }\{\frac{1}{3}(1-c_2)x^3+\left(\frac{5c_2}{2}-c_1\right)x$$

$$-\frac{2}{3}(8c_2-4c_1+c_0)+\mathcal{O}\left(\frac{1}{x}\right)\}.$$

$$\epsilon \to a\epsilon \text{and}\varphi (q)\to a\varphi (aq)$$

$${\epsilon }_F=\epsilon /\sqrt{-{\varphi }_{\max }^{\prime }},$$

$${\epsilon }_0\equiv \frac{a_0}{2\pi }\epsilon ={\int }_0^{\infty }\left[1-r\psi (r)\right]\mathrm{d}r,$$

$${\epsilon }_2={\left|\frac{3a_2}{\pi }\right|}^{1/3}\epsilon .$$

$$V_{\mathrm{pert}}(r,\theta ,t)=V_{\mathrm{pert}}(r){\mathrm{e}}^{\mathrm{i}(m\theta -\omega t)}$$

$${\Sigma }_{\mathrm{resp}}(r,\theta ,t)=\int f_{\mathrm{resp}}(r,\theta ,v_r,v_{\theta },t)\mathrm{d}{v}_r\mathrm{d}{v}_{\theta }$$

$$V_{\mathrm{resp}}(\mathbf{x})=G\int {\Sigma }_{\mathrm{resp}}({\mathbf{x}}^{\prime })\psi (|\mathbf{x}-{\mathbf{x}}^{\prime }|){\mathrm{d}}^2{\mathbf{x}}^{\prime },$$

$$V_{\mathrm{pert}}(r,\theta ,t)\equiv V_{\mathrm{resp}}(r,\theta ,t)$$

Full text

We explore four softening kernels with different degrees of gravity bias, and with/without compact support (compact in the sense that they yield exactly Newtonian forces outside the softening kernel). These kernels are applied to two disc galaxy models with well-known unsoftened unstable modes. We illustrate quantitatively the importance of a vanishing linear gravity bias to yield accurate frequency estimates of the unstable modes. As such, Plummer softening, while very popular amongst simulators, performs poorly in our tests.

The best results, with excellent agreement between the softened and unsoftened mode properties, are obtained with softening kernels that have a reduced gravity bias, obtained by compensating for the sub-Newtonian forces at small interparticle distances with slightly super-Newtonian forces at radii near the softening length. We present examples of such kernels that, moreover, are analytically simple and computationally cheap. Finally, these results light the way to the construction of softening methods with even smaller gravity bias, although at the price of increasingly complex kernels.

††pagerange: How gravitational softening affects galaxy stability I. Linear mode analysis of disc galaxies–A.2††pubyear: 2009

1 Introduction

The evolution of a collisionless stellar system is determined by the collisionless Boltzmann equation (CBE)

[TABLE]

Here, $F(\boldsymbol{x},\boldsymbol{v},t)$ is the distribution function (DF), which gives the stellar phase-space density at location $(\boldsymbol{x},\boldsymbol{v})$ and time $t$. The motion of the collisionless fluid is controlled by the (positive) binding potential $V(\boldsymbol{x},t)$. The direct numerical integration of the CBE in six-dimensional phase space is in general impossible because under the CBE the DF develops ever finer structures owing to phase mixing or chaotic mixing. However, numerical schemes which smooth out such fine structure (whereby violating the CBE) are possible but taxing (Yoshikawa, Yoshida & Umemura, 2013; Schaller et al., 2014; Colombi et al., 2015).

A much more popular method for modelling collisionless stellar dynamics is an $N$-body simulation: a Monte-Carlo approach, which integrates the CBE via the method of characteristics. The DF is represented by a collection of phase-space points, called “particles”, and each particle is evolved through phase space along its characteristic curve by solving the first-order differential equations

[TABLE]

In the case of gravitational forces, the binding potential can be written as

[TABLE]

with $\rho$ the stellar mass density and $m_{i}$ and $\boldsymbol{x}_{i}$ the mass and position, respectively, of the $i^{\text{th}}$ particle. Using the approximation (3b) to the gravitational field, diverging accelerations may occur in close encounters (‘collisions’) between particles. Such collisions are an artefact of the much smaller $N$ in the simulation than in the simulated system. This problem is generally solved by “softening” gravity, when the $1/r$ potential in equations (3) is replaced by a non-diverging form,

[TABLE]

where $\psi(r)$ is the softening Green’s function, $\varepsilon$ the softening length, and $\phi$ the dimensionless softening kernel. For suitable functions $\phi$ this modifies the inter-particle interactions such that accelerations remain bounded and strongly deflecting encounters are avoided.

Unfortunately, this force modification results in a bias of gravity and hence changes the character of the physical problem being addressed by the $N$-body simulation. In practice, a balance must be found between too much softening, causing force bias, and too little softening, allowing strong encounters that render the $N$-body dynamics intractable and collisional (Merritt, 1996). Dehnen (2001) has derived asymptotic relations in the context of spherically symmetric systems, which can be used to inform the choice of the softening parameters (kernel and length) such as to minimise the resulting mean-square gravity error.

However, it remains unclear what the optimal choice of these parameters is in terms of accurately modelling the dynamics, rather than merely minimising the gravity error. The goal of this series of papers is to investigate this question by considering stellar dynamical problems which invoke non-trivial dynamics but are still simple enough that accurate solutions of the CBE are available. Specifically, we consider unstable two- and three-dimensional systems, whose eigenmodes can be accurately obtained from linear perturbation theory. In this first paper, we focus on the two-armed (multiplicity $m=2$) spiral-shaped eigenmodes of two-dimensional razor-thin disc galaxies in the limit of $N\to\infty$ (which can be achieved with computationally cheap linear mode analysis without $N$-body simulations).

Research into the origin and longevity and/or transience of spiral structure in disc galaxies has by and large relied on $N$-body simulations (Hohl, 1971; Toomre, 1977; Sellwood & Lin, 1989; Sellwood & Kahn, 1991; Sellwood, 2011; D’Onghia, Vogelsberger & Hernquist, 2013; Sellwood & Carlberg, 2014) and on (semi-)analytical mode analysis using the first-order CBE as a starting point (Zang, 1976; Kalnajs, 1977; Fridman & Polyachenko, 1984; Palmer, 1994; Vauterin & Dejonghe, 1996; Pichon & Cannon, 1997; Evans & Read, 1998; Jalali & Hunter, 2005; Jalali, 2007; Binney & Tremaine, 2008; Polyachenko & Just, 2015). Here, we use linear theory to investigate how gravitational softening affects the growth of small-amplitude eigenmodes in $N$-body simulations.

We introduce the various softening techniques explored by us in section 2. Our implementation of gravitational softening in linear mode theory is layed out in section 3. The properties of the unsoftened axially symmetric disc models are discussed in section 4. Our results are presented in section 5 and their implications are discussed in section 6.

2 Gravitational softening in two spatial dimensions

The usual motivation for softening in two-dimensional $N$-body simulations is to account for the finite thickness of the stellar disc, which is neglected in the razor-thin limit (Miller, 1971). In this interpretation, the softening length $\varepsilon$ is no longer a numerical but a physical parameter and the modification of gravity and, consequently, of the dynamics is deliberate, because real galaxies are not razor-thin (Romeo, 1992, 1997).

However, as we are interested in the errors introduced by the softening-induced modification of gravity, we cannot adopt this interpretation but must consider the razor-thin disc as the desired physical model, whose modes one may attempt to recover with an $N$-body simulation.

2.1 Gravity bias

When inserting the softening kernel (4) into equation (3b), we obtain the softened potential

[TABLE]

Here, $\hat{V}(\boldsymbol{x})$ can be interpreted as an estimate, based on the masses and positions of the simulation particles, for the true potential $V(\boldsymbol{x})$. The mean-square error made by this estimate can be decomposed into a variance,

[TABLE]

and a bias,

[TABLE]

Here, $\langle\cdot\rangle$ denotes the ensemble average of a quantity over all $N$-body realisations, in the limit $N\to\infty$, of the underlying smooth density distribution. The variance measures the mean amplitude of the random fluctuations of a softened $N$-body potential around its ensemble mean. In other words, it measures the graininess of the $N$-body potential: the error made by not softening enough. The bias measures the deviation of the ensemble mean of the softened $N$-body potential from the underlying smooth potential. This is the error made by softening too much.

For the situation of three-dimensional $N$-body simulations, Dehnen (2001) has derived analytical asymptotic relations for these quantities. An adaptation of his derivations to an $N$-body simulation of a razor-thin disc with surface density $\Sigma(\boldsymbol{x})$ gives

[TABLE]

for the bias on the potential $V$ (see Appendix A for a derivation). The coefficients $a_{n}$ depend only on the functional form of the softening kernel:

[TABLE]

This is different for three-dimensional systems, where this bias asymptotes as $\varepsilon^{2}$ at lowest order. Here, in two dimensions, the gravity biases are proportional to $\varepsilon$. This is a direct consequence of the reduced number of dimensions. Thus, in two-dimensional $N$-body simulations, the gravity bias is in general significantly stronger than in three-dimensional $N$-body simulations111 Except in three-dimensional simulations of disc galaxies with scale height $h\apprle\varepsilon$, since in three dimensions

$\displaystyle\mathrm{bias}_{\boldsymbol{x}}(\hat{V})$ $\displaystyle\approx-a_{0}^{\mathrm{3D}}\varepsilon^{2}G\rho(\boldsymbol{x})$

(10a)

$\displaystyle\mathrm{bias}_{\boldsymbol{x}}(\hat{\boldsymbol{a}})$ $\displaystyle\approx-a_{0}^{\mathrm{3D}}\varepsilon^{2}G\boldsymbol{\nabla}\rho(\boldsymbol{x})\propto\varepsilon^{2}/h$

(10b)

(Dehnen, 2001, eqs. 10) with $\rho(\boldsymbol{x})$ the spatial density..

Note that the integrand in equation (9) is always well-behaved in the limit $u\to 0$ but can be problematic for large $u$. If $1-u\phi(u)\propto u^{-p}$ for large $u$ then only coefficients $a_{n}$ with $n<p-1$ are finite; the rest come out infinite. In this case, the Taylor series (8) does not converge, but the terms at $n<p-1$ still provide a useful approximation, only the remainder grows faster than $\varepsilon^{2\lceil(p-1)/2\rceil}$.

2.2 Softening kernels

The functional forms and other properties of the softening kernels used in this study are listed in Table 1. Their interparticle interaction potentials and forces are plotted in Figure 1. In particular, we list the dimensionless interparticle interaction potential $\phi$ and the corresponding three-dimensional and two-dimensional dimensionless density kernels, denoted by $\varrho$ and $\sigma$, respectively, assigned to each point particle.

In three dimensions, gravitational softening is equivalent to estimating the spatial mass density as a superposition of spheres with density distribution $\rho(r)$ placed at the particle positions. In two dimensions, one can consider softening as a way of smoothing the overall mass distribution as a superposition of razor-thin discs with density distribution $\Sigma(r)$ at the particle positions. For a softening kernel $\phi$, these spherical or razor-thin surface density distributions are given by

[TABLE]

where the dimensionless density and surface-density kernels are given by

[TABLE]

respectively (Binney & Tremaine, 2008; Kalnajs, 1999).

Note that the softened force is that between a softened particle (with density distribution $\rho(r)$ or $\Sigma(r)$ given by eqn. (11)) and a point particle and not the force between two softened particles (Barnes, 2012). This remark applies to all softening techniques.

2.2.1 Plummer softening P0: infinite support, $a_{0}\neq 0$

In three dimensions, this popular softening method corresponds to estimating the spatial mass density as a superposition of Plummer spheres (Plummer, 1911) with scale radius $\varepsilon$ at the particle positions. In two dimensions, Plummer softening amounts to smoothing the overall mass distribution as a superposition of razor-thin Kuz’min discs (Kuz’min, 1956).

Plummer softening modifies the gravitational interaction at all interparticle separations and asymptotes to the Newtonian interaction only for infinitely large particle separations. As a result, $a_{0}^{\mathrm{3D}}=\infty$ and the gravity bias (10) in this case grows faster than $\propto\varepsilon^{2}$. In two dimensions, this method’s non-zero $a_{0}$ indicates that the gravity bias increases linearly with softening length.

2.2.2 Modified Kuz’min softening Q2: infinite support, $a_{0}=0$

One can modify the Plummer kernel such that $a_{0}$, $a_{2}$, …, $a_{k}=0$, for any chosen even $k$, in order to significantly reduce the gravity bias. For instance, we introduce the class of modified Kuz’min potentials, with an interaction potential given by

[TABLE]

here $K_{n}(u^{2})=\sum_{i=0}^{n}c_{i}u^{2i}$ is a polynomial of degree $n$ in $u^{2}$ with coefficients $c_{i}$. We refer to the $n^{\rm th}$ member of this class as Qn. Clearly, the choice $n=0$ yields the Plummer kernel, or: Q0=P0.

The extra degree of freedom that comes with the choice $n=1$, can be exploited to make $a_{0}=0$. For instance, for this Q1 kernel, one finds that

[TABLE]

and the choice $c_{0}=2$, $c_{1}=1$ makes $a_{0}=0$. Thus, the ‘Q1’ method is defined by the interparticle potential

[TABLE]

and by the corresponding 3D and 2D density distributions

[TABLE]

In order to achieve $a_{0}=0$, this softening method compensates with slightly super-Newtonian forces at $r\apprge 1.2\varepsilon$ for the substantially sub-Newtonian accelerations at small separations. Unfortunately, while we now have $a_{0}=0$, the Q1 softening kernel still has a diverging second-order coefficient $a_{2}$, such that the gravity bias grows faster than $\varepsilon^{3}$.

The extra degree of freedom provided by coefficient $c_{2}$ of kernel Q2 can be used to provide $a_{2}$ with a finite value. Indeed, for this kernel, we find that

[TABLE]

and

[TABLE]

Demanding $a_{0}$ to be zero and $a_{2}$ to be finite, leads to $c_{2}=1$, $c_{1}=5/2$, and $c_{0}=3$. All properties of this Q2 kernel are listed in Table 1.

Of course, this game can be continued to obtain $a_{2}=0$, then finite $a_{4}$, then $a_{4}=0$, etc. However, the functional form of the resulting interaction potentials becomes increasingly complex.

2.2.3 Ferrers $n=3$ softening F3: compact support, $a_{0}\neq 0$

An increasingly popular method is cubic spline softening, which in three-dimensional $N$-body simulations corresponds to replacing each particle by a cubic-spline smoothing kernel as is widely used in Smoothed Particle Hydrodynamics (SPH) codes (see e.g. Monaghan, 1992) and was introduced as a gravitational softening kernel for $N$-body/SPH codes by Hernquist & Katz (1989). Its main advantages in this context are its exactly Newtonian behaviour beyond the softening length and its dual use as hydrodynamics smoother and gravity softener. However, its interparticle potential and 3D density distribution are numerically rather unattractive due to their complex, piecewise continuous functional forms.

Here and in the remainder, we will use the phrase “compact support” to indicate that a kernel yields exactly Newtonian forces outside the softening kernel. In three dimensions, it immediately follows from Newton’s first theorem that the corresponding density distribution is zero outside the kernel, i.e. $\varrho=0$ at $u>1$. In two dimensions, this is not the case. In fact, all softening kernels with exact Newtonian gravity at separations $r>\varepsilon$ have poorly behaved corresponding razor-thin disc profiles $\sigma(u)$ (12b), with infinite spatial extent and negative values.

We here opt for the so-called Ferrers softening methods, labeled ‘Fn’, whose interaction potentials are polynomials of degree $n+1$ in the variable $t=1-u^{2}$ inside the softening length and that behave exactly Newtonian at separations $r>\varepsilon$. In three dimensions, they correspond to replacing each particle with a Ferrers (1877) sphere of order $n$. For $n=0$ this is just a homogeneous sphere. Higher-order models have spherical densities that are simple polynomials, with $n$ continuous derivatives, in the variable $t=1-u^{2}$.

For this paper, we investigate member $n=3$ from this family, with properties listed in Table 1, as an example of a softening method with compact support but with $a_{0}\neq 0$.

2.2.4 Modified Ferrers softening L2: compact support, $a_{0}=0$

It is possible to modify the Ferrers softening methods to obtain $a_{0}=0$ in order to reduce their gravitational bias while retaining the attractive property of having compact support. Here, we test the method labeled ’L2’, or ‘2D modified Ferrers $n=2$’, in Table 1. Its name derives from the fact that it’s based on the F2 kernel, which has as an interaction potential $\phi(u)=1+\frac{1}{2}t+\frac{3}{8}t^{2}+\frac{5}{16}t^{3}$, but where the coefficient of the last term is tuned to make $a_{0}=0$. This leads to the L2 interaction potential $\phi(u)=1+\frac{1}{2}t+\frac{3}{8}t^{2}+\frac{5}{2}t^{3}$.

This kernel achieves its desirable properties by having super-Newtonian accelerations for a limited range of separations close to and inside of $r=\varepsilon$.

2.3 Softening scale

The softening length and kernel are only defined up to a re-scaling: the softened potential (5) is invariant under the transformation

[TABLE]

with scaling factor $a$. This implies that the parameter $\varepsilon$ has no natural scale by itself and comparing different kernels at the same $\varepsilon$ is meaningless. Therefore, some other measure is required for such a comparison. One such measure valid for all softening kernels is the force-scaling as

[TABLE]

where $-\phi^{\prime}_{\max}$ denotes the maximum value of the derivative of $\phi$ (Springel, 2010). The ratios $\varepsilon_{F}/\varepsilon$, scaled to unity for the Plummer kernel, are listed in Table 1 for the kernels considered in this study.

For kernels with $a_{0}\neq 0$, another natural measure of the softening length is

[TABLE]

which measures the actual scale of the bias irrespective of any re-scaling. With this definition, $\varepsilon_{0}=\varepsilon$ for Plummer softening. For other softening methods used in this study, $\varepsilon_{0}$ is given in Table 1.

Likewise, softening techniques with zero $a_{0}$ but non-zero $a_{2}$ can be scaled to a common level of gravity bias via the softening length transformation

[TABLE]

With this definition, $\varepsilon=\varepsilon_{2}$ for Q2 softening.

3 Softened gravity in stability analysis

3.1 Linear mode theory

We use pyStab, a Python/C++ computer code, to analyse the stability of a razor-thin stellar disc embedded in an axially symmetric gravitational potential. The details of the mathematical formalism behind this code and of its implementation can be found in Vauterin & Dejonghe (1996), Dury et al. (2008), and De Rijcke & Voulis (2016) so we will not repeat these here. An axially symmetric disc galaxy model is characterized by a distribution function $F_{0}(E,J)$, with $E$ the specific binding energy and $J$ the specific angular momentum of a stellar orbit, and a mean gravitational potential $V_{0}(r)$. In the remainder, we will refer to this unperturbed axially symmetric state as the “base state” of the system. Note that this base state is only a correct solution of the CBE when employing the Newtonian gravitational interaction (but see below).

For any given base state, pyStab can retrieve those complex frequencies $\omega$ for which a spiral-shaped perturbation of the form

[TABLE]

constitutes an eigenmode. Here, $(r,\theta)$ are polar coordinates in the stellar disc, $m$ is the multiplicity of the spiral pattern, $\Omega_{p}=\Re\{\omega\}/m$ its pattern speed, and $\Im\{\omega\}$ its growth rate. A general perturbing potential can always be expanded in such modes and, owing to the linear approximation, these can be studied independently from each other.

In essence, pyStab solves the first-order CBE to find the response distribution function $f_{\rm resp}(r,\theta,v_{r},v_{\theta},t)$ produced by a given perturbation $V_{\rm pert}(r,\theta,t)$. This response distribution function generates the response density

[TABLE]

which in turn gives rise to the response softened gravitational potential

[TABLE]

where the integral runs over the whole surface of the stellar disc and $\psi(r)$, defined in equation (4), is the softened Green’s function for gravitational interactions, replacing the Newtonian $1/r$. Eigenmodes are then identified by the fact that

[TABLE]

and pyStab employs a matrix method (Kalnajs, 1977) to find them. The perturbing potential $V_{\rm pert}$ is expanded in a basis of potentials, $V_{\ell}$. The response to each basis potential, denoted by $V_{\ell,\rm resp}$, can likewise be expanded in this basis as

[TABLE]

If the perturbation is an eigenmode, then the ${\mathcal{C}}$ matrix can be shown to possess a unity eigenvalue (Vauterin & Dejonghe, 1996; Dury et al., 2008; De Rijcke & Voulis, 2016). This feature is exploited by pyStab to identify the eigenmodes.

The formalism contains a number of technical parameters, such as the number of orbits on which phase space is sampled (here we use $n_{\rm orbit}(n_{\rm orbit}+1)/2$ orbits with $n_{\rm orbit}=600$ in the allowed triangle of turning point – or pericentre/apocentre – space), the number $n_{\rm Fourier}$ of Fourier components in which the periodic part of the perturbing potential is expanded (here we use $n_{\rm Fourier}=80$), the number of potential-density pairs (PDPs) that is used for the expansion of the radial part of the perturbing potential and density (we use 44 PDPs), and the shape and extent of the PDP density basis functions. As in De Rijcke & Voulis (2016), we use PDP densities of the form

[TABLE]

where the average radii $r_{\ell}$ cover the relevant part of the stellar disc and are evenly spaced on a logarithmic scale so the resolution is highest in the inner regions of the disc. The widths $\sigma_{\ell}$ are automatically chosen such that consecutive basis functions are sufficiently unresolved to represent any smooth function. The position of these PDP density basis functions can be tuned to achieve a high spatial resolution there where the eigenmodes live. The corresponding PDP potentials are obtained via

[TABLE]

3.2 Introducing gravitational softening

Since we want to validate our approach by comparing particular results with published work based on numerical simulations, we mimic the strategies employed by simulators when setting up and performing $N$-body simulations of disc galaxies aimed at mode analysis. Usually, an initial condition is generated by sampling stellar particles from the distribution function $F_{0}(E,J)$ evaluated using the Newtonian gravitational potential $V_{0}(r)$, independent of the gravitational softening that is employed later on when evolving the particles through time. Moreover, the axially symmetric force field of the base state is subsequently evaluated correctly, i.e. without softening, either by directly using the analytical expression for the potential $V_{0}$ or by adding a small correction to the softened gravitational field derived from the particles. Only the non-axisymmetric force field of the growing waves is softened (Earn & Sellwood, 1995; Sellwood & Evans, 2001; Sellwood, 2012). This allows a simulator to sample particles from the correct DF evaluated in the correct potential so that at least the initial conditions of a simulation correspond to the intended base state and the particle dynamics in the axially symmetric force field is followed correctly.

Therefore, we only implement gravitational softening in the response potential $V_{\rm resp}(r,\theta,t)$, but not in the axially symmetric base state potential $V_{0}(r)$. Using this strategy, eqn. (25), and a fortiori eqn (29), is the only place where the softened gravitational interaction enters the computation of the modes. It is therefore straightforward to insert interaction potentials other than the Newtonian one into a mode analysis code. The gravity bias introduced in Section 2.1 must then be regarded as a measure for the fidelity with which the softened response potential resembles the Newtonian one. Thus, we can use linear stability theory to emulate the results expected in the large $N$ limit from $N$-body simulations of disc galaxies. In this paper, we investigate the effect on the eigenmodes in disc galaxy models from the P0, Q2, F3, and L2 softening methods listed in Table 1.

4 The base states

Below, we give the essential details of the two base states that we employ for this study. We also list the frequencies of the known eigenmodes of these base states computed for a Newtonian interparticle interaction.

4.1 The isochrone disc model

The isochrone disc is characterized by the cored density profile

[TABLE]

which self-consistently generates the gravitational potential

[TABLE]

(Henon, 1959; Kalnajs, 1976). Here, $M$ is the total mass of the stellar disc and $b$ its scale-length. As shown by Kalnajs (1976) and Earn & Sellwood (1995), a family of distribution functions that generate this potential-density pair is given by

[TABLE]

with $m_{K}$ an integer, $x=J\sqrt{2E}/GM$, and

[TABLE]

Here, $P_{m}$ is the Legendre polynomial of degree $m$ and

[TABLE]

We adopt $m_{K}=12$ for this study. The Legendre polynomial can be evaluated explicitly, allowing the integral featuring in the expression for the distribution function to be evaluated in closed form. However, the resulting expression is numerically very unstable for small $x$. Therefore, we opted to simply evaluate the integral numerically. This distribution function is only used to populate orbits with positive angular momentum.

Counter-rotating stars have been added according to the prescription given in Earn & Sellwood (1995) which, unfortunately, necessitates going back and forth between energy, angular momentum and the radial action variable, $J_{r}$:

[TABLE]

with $E^{\prime}$ the energy corresponding to a radial action $J_{r}+|J|$ and zero angular momentum. Fortunately, analytical conversion formulæ between energy, angular momentum, and radial action exist for the isochrone disc (Binney & Tremaine, 2008).

We will focus here on the bisymmetric ($m=2$) modes of this model. Pichon & Cannon (1997) provide the frequencies of three modes of this base state model, choosing units such that $G=M=b=1$, as

[TABLE]

while Jalali & Hunter (2005) find

[TABLE]

for the two main modes222 Accuracy estimates for mode frequencies derived from linear stability computations are hard to obtain since many numerical parameters come into play. Judging from the differences between published mode frequencies and from our own limited experiments with varying the values of the employed numerical parameters (described in Section 3.1) we estimate the mode frequencies to be accurate to about the percent level. .

4.2 The Mestel disc

The Mestel (1963) disc has a cusped total surface density given by

[TABLE]

which self-consistently generates a gravitational potential of the form

[TABLE]

with the surface density scale given by $\Sigma_{0}=v_{0}^{2}/2\pi Gr_{0}$. Here, $v_{0}$ is the value of the disc’s constant circular velocity. A central hole is cut out of this disc model by multiplying its distribution function (Toomre, 1977)

[TABLE]

where $q$ is a real number, with a cut-out function of the form

[TABLE]

with $x=\left(J/r_{0}v_{0}\right)^{n}$ (obviously, this also slightly suppresses the distribution function at larger $J$-values). Outside the central cut-out region, the disc’s constant radial velocity dispersion is given by $\sigma=v_{0}/\sqrt{1+q}$.

Here, we will focus on the $q=6$, $n=4$ member of this model family and adopt units such that $G=v_{0}=r_{0}=1$. Its dominant bisymmetric mode is then expected to have a frequency

[TABLE]

as shown, e.g., by Toomre (1977); Read (1997); Evans & Read (1998); Polyachenko & Just (2015). This mode owes its existence to the inner cutout: if the angular momentum cut-off is not sufficiently steep, i.e. if $n$ is too small, there is no eigenmode. The idea is that incoming trailing wave packets are (partially) reflected from this sharp inner edge to travel back outwards as leading waves. Overreflection, or swing amplification, at the evanescent zone around the corotation resonance (Mark, 1976; Toomre, 1981) sends amplified trailing waves back inwards. Inside this resonance cavity, growing modes can occur (Evans & Read, 1998).

5 Results

5.1 The isochrone disc

5.1.1 Unsoftened gravity

For the two most rapidly growing $m=2$ modes of the $m_{K}=12$ isochrone disc, we find frequencies

[TABLE]

in good agreement with published values.

However, the third mode listed by Pichon & Cannon (1997) showed up in our analysis as only the fifth fastest growing mode, with a frequency

[TABLE]

The two interloping modes at

[TABLE]

have not been described in the literature before. We confirmed that they are robust to changes of the numerical parameters in the code (resolution in phase space, number of Fourier modes, etc.) and that they exert a zero total torque on the disc, as they should, and therefore see no reason to discard them as spurious (Polyachenko & Just, 2015).

5.1.2 Softened gravity

For this base state, published information on how the properties of the two main $m=2$ modes change with softening in a numerical simulation exists. Earn & Sellwood (1995) use a polar grid code with 120,000 particles, a fixed time step, and a polar grid of 128 azimuthal and 85 radial nodes to simulate the $m_{K}=12$ isochrone disc from quiet-start initial conditions (Sellwood, 1983) using Plummer softening with different softening lengths. The results of these simulations are shown in Fig. 2 as black data points. Both the pattern speed and the growth rate of the two dominant modes appear to be declining functions of softening length. As Earn & Sellwood (1995) note: “[…] it is clear that both parts of the eigenfrequency are strongly affected by even moderate softening.” The dependence of the mode frequencies on softening length is markedly non-linear which hampers a simple extrapolation to zero softening length.

Overplotted in Fig. 2 are the results from our linear stability analysis with pyStab, using different softening prescriptions. Clearly, our results for Plummer softening agree rather well with those presented in Earn & Sellwood (1995): both the pattern speed and growth rate are non-linearly declining functions of the softening length (scaled according to equations (19) and (20) to a common maximum inter-particle force). The drop is steepest for small softening lengths and becomes shallower for larger $\varepsilon$. Especially for larger $\varepsilon$-values, numerical simulations and linear mode analysis predict the same behaviour for frequency as a function of softening length. We tentatively attribute the deviations between the simulations and linear theory at small $\varepsilon$-values to variance, i.e. to the gravity error caused by not softening enough (cf. paragraph 2.1).

Using the other softening recipes, the pattern speed and growth rate likewise decline with increasing softening length but they do so much less dramatically and with smaller deviation from a linear dependence on softening length than when using Plummer softening. Moreover, it appears that having a zero gravity bias parameter $a_{0}$ induces a much stronger effect than having compact support. This is exemplified in this case by the Q2 method (infinite support, $a_{0}=0$) yielding results much closer to the Newtonian ones than the F3 method (compact support, $a_{0}>0$). Methods that combine compact support with having $a_{0}=0$, like the L2 method, appear vastly superior, with very little deviation between the retrieved mode frequencies and the correct, Newtonian values.

However, we refer the reader to section 6 for a discussion of how to correctly interpret this apparent success.

5.2 The Mestel disc

5.2.1 Unsoftened gravity

pyStab retrieves the dominant mode of the $q=6$, $n=4$ Mestel disc, along with a number of much slower growing modes. Since the frequencies of these minor modes are sensitive to the choice of numerical parameter values, they are most likely spurious (Polyachenko & Just, 2015). Using unsoftened gravity, we find the dominant mode to have a frequency $\omega=0.876+0.128\;\mathrm{i}$, which is in good agreement with the values reported by Toomre (1977); Read (1997); Polyachenko & Just (2015) and which were computed using different mode analysis techniques and codes.

5.2.2 Softened gravity

In Fig. 3, we show how the pattern speed (top panel) and growth rate (bottom panel) of the dominant mode of this base state change with increasing softening length (scaled according to equations (19) and (20) to a common maximum inter-particle force) using the softening recipes listed in Table 1.

Overplotted in this figure, we show the frequency estimates of Sellwood & Evans (2001) for this Mestel disc model, based on $N$-body simulations with a particle-mesh code employing 2.5 million particles, and a grid of 256 azimuthal and 200 radial nodes. The 5 simulations presented here all start from exactly the same initial conditions but are evolved using different Plummer softening lengths. The agreement with our linear mode analysis is not as good as in the case of the isochrone disc. As reported by Sellwood & Evans (2001), there is a $\sim 10$ % scatter between the measured frequencies of simulations with resampled initial conditions at a constant particle number. This may be why the simulation datapoints do not converge to the Newtonian linear-mode result for zero softening length. Moreover, particle noise may have negatively affected the frequency measurements. Still, the trend followed by these simulations is in qualitative agreement with our results: both the pattern speed and the growth rate decrease with increasing softening length.

As for the isochrone disc, softening methods with $a_{0}=0$ (like Q2) stay much closer to the Newtonian mode frequency than methods with compact support but non-zero $a_{0}$ (like F3) for a given value of the softening length $\varepsilon$. Methods that combine compact support with $a_{0}=0$ (like L2) generally outperform the others.

Again, we refer the reader to section 6 for a discussion of how to correctly interpret this apparent success.

6 Discussion

6.1 Scaling to the same level of gravity bias

As mentioned in paragraph 2.3, the softening length and kernel are only defined up to a re-scaling and we advocate the scale

[TABLE]

to bring methods with non-zero $a_{0}$ to a common gravity bias level.

In Fig. 4, we show the retrieved frequencies of the modes of the isochrone and Mestel discs as a function of $\varepsilon_{0}$ for the two softening methods with non-zero gravity bias parameter $a_{0}$ (i.e. P0 and F3). Clearly, the differences between both softening methods, which are so striking in Figures 2 and 3, now largely disappear. At a given $\varepsilon_{0}$-value, all $a_{0}\neq 0$ methods perform almost equally well. Thus, it is always possible to re-scale the softening length of one $a_{0}\neq 0$ softening technique such that it approximately matches the performance of another $a_{0}\neq 0$ method. No exact matching is possible because of the higher order terms in the expansion of the gravity bias.

As can be seen in Figure 5, softening techniques with $a_{0}=0$ but non-zero $a_{2}$ can be scaled to a common level of gravity bias using the transformation (22), allowing for higher order terms in the expansion (8) for the gravity bias.

Based on these results, it seems fair to say that softening strategies with $a_{0}=0$ generally yield more accurate (i.e. Newtonian-like) mode frequencies than strategies with $a_{0}\neq 0$ because the gravity bias of the latter grows linearly with softening length $\varepsilon$ while for the former it grows much more slowly, as $\varepsilon^{3}$. However, within each of these classes of softening techniques, there is no particular reason to favour one method over another provided they are compared at (approximately) the same level of gravity bias.

6.2 Physical interpretation

We define the two-dimensional Fourier transform $\widehat{\psi}(k)$ of the interparticle interaction potential $\psi(r)$ as

[TABLE]

Based on Poisson’s equation, using separation of variables it is straightforward to show that the radial part of the gravitational response potential, which we denote here by $V_{m}(r)$, generated by an $m$-armed spiral response density of the form

[TABLE]

can be retrieved from the relation

[TABLE]

with ${\mathcal{H}}_{m}$ the Hankel transform of order $m$ (see e.g. Binney & Tremaine, 2008). The Hankel transform of order $m$ of a function $f(r)$ is defined as

[TABLE]

with $J_{m}(x)$ a Bessel function of the first kind.

Here, we will use the notation $\widehat{\psi}_{N}(k)$ for the Fourier transform of the Newtonian interaction potential, with

[TABLE]

Likewise, for the interaction potentials listed in Table 1 we find that:

[TABLE]

For the softening techniques with compact support, F3 and L2, no simple analytical expression exists for the Fourier transform of their interaction potentials but they can easily be obtained numerically.

For a given response density $\Sigma_{m}$, the softened response potential $V_{m}^{\varepsilon}$ and the unsoftened response potential $V_{m}^{0}$ are connected as

[TABLE]

with $\varepsilon$ the softening length. We plot the $k$-dependent suppression factor $\widehat{\psi}/\widehat{\psi}_{N}$ that links the Fourier transforms of the softened and unsoftened response potentials in Figure 6. The most striking consequence of gravitational softening is the suppression of the small-scale, i.e. large wavenumber $k$, structure in the Fourier transform of the response potential.

In Appendix A, we show how this suppression factor is connected to the gravity bias. More specifically, we prove that the even coefficients in the series expansion of the suppression factor around zero $k$ are directly proportional to the even coefficients in the series expansion of the gravity bias around zero $\varepsilon$ (in case the latter exist). Clearly, desiging a softening kernel to have vanishing bias coefficients is equivalent to designing an interaction kernel whose suppression factor $\widehat{\psi}/\widehat{\psi}_{N}$ is increasingly close to unity for small wavenumbers $k$.

However, just as it is not sensible to compare the various softening strategies as we did in Figures 2 and 3, it makes little sense to compare the suppression factors as a function of $k\varepsilon_{F}$. It is more meaningful to compare the suppression factors at the same level of gravity bias, i.e. as a function of $k\varepsilon_{0}$ for the softening methods with $a_{0}\neq 0$, and as a function of $k\varepsilon_{2}$ for the softening methods with $a_{0}=0$. This comparison is shown in Figure 7 and restates our previous conclusions. As a function of the re-scaled wave number $k\varepsilon_{0}$, which places the P0 and F3 kernels on an equal gravity bias footing, the suppression factors of these two softening techniques behave remarkably similar. In fact, Plummer softening leads to less suppression for large $k\varepsilon_{0}$ than F3 softening. This agrees with Figure 4 in which Plummer softening is shown to stay closer to the correct, Newtonian result than F3 softening at an equal level of gravity bias. The Q2 kernel, in turn, leads to less suppression than the L2 kernel and, as can be seen in Figure 5, it also leads to slightly better frequency estimates.

This suppression of the response potential will likely lead to an increased stability of the model galaxy. This expectation is borne out by studying the stability of axially symmetric WKBJ waves under Plummer softening, where a Toomre $Q$-value $Q<1$ now separates growing from stationary waves (Miller, 1971). This analysis has been extended to include general softening kernels by Romeo (1994, 1997). The physical background of our results and the expected influence of softening on $m=0$ WKBJ waves are, therefore, already well understood. Here, we took this work further by studying general eigenmodes beyond the WKBJ approximation and by going to $m=2$ patterns.

7 Conclusions

We use linear perturbation theory to investigate how different recipes for gravitational softening, as employed in numerical $N$-body simulations of razor-thin disc galaxies, affect predictions for the properties of the latter’s spiral eigenmodes. We specifically focus on the frequencies, i.e. pattern speeds and growth rates, of two-armed modes in the linear regime. We warn the reader that our approach does not take into account the effects of approximate force evaluations (Barnes & Hut, 1986), finite-$N$ (Dehnen, 2001; Sellwood, 2012; Fouvry et al., 2015), stochasticity (Sellwood & Debattista, 2009), implicit softening contributed by the grid in particle-mesh codes (Romeo, 1994), etc. whose respective magnitudes, moreover, may depend on the amount of gravitational softening.

We have tested our linear mode analysis approach by comparing the behaviour of the frequencies of the dominant $m=2$ modes of an isochrone disc and of the Mestel disc as a function of Plummer softening length with those found in the $N$-body simulations reported by Earn & Sellwood (1995) and Sellwood & Evans (2001). Overall, we found reasonably good agreement between linear theory and numerical simulations, also in the softened regime.

We argue that the only meaningful way of comparing softening kernels is to scale them to the same gravity bias level. In this paper, we show how this scaling can be achieved, based on the results of Dehnen (2001). Thus, it is always possible to re-scale the softening length of one softening technique such that it matches the performance of another method with the same dependence of gravity bias on softening length.

We have shown that softening methods with a vanishing lowest-order term in the expansion of the gravity bias as a function of softening length (in two dimensions, this is a linear term; in three dimensions, this term is quadratic in the softening length) and whose gravity bias therefore grows slowly with increasing softening length (e.g. the Q2 and L2 methods discussed in this paper) provide more accurate mode frequency estimates than methods with a non-zero lowest-order term (e.g. the P0 and F3 methods). Softening methods with zero lowest-order term compensate the sub-Newtonian forces deep inside the kernel with super-Newtonian forces near $r\sim\varepsilon$.

Kernels with compact support, in the sense that they yield exactly Newtonian forces outside of the softening kernel, perhaps somewhat counter-intuitively, do not necessarily provide more accurate frequency estimates than kernels with infinite extent. For instance, when compared at a common level of gravity bias, the Plummer kernel (P0) provides more accurate frequency estimates than the F3 kernel. Likewise, the Q2 kernel outperforms the L2 kernel in this regard.

The relative merit of a softening kernel can be judged from its suppression of the small-scale, i.e. large wavenumber $k$, structure in the Fourier transform of its response potential. The stronger this suppression, measured at a given level of gravity bias, the more the mode frequency estimates deviate from their Newtonian values.

As a guide to simulators, we provide an example of how a softening technique, in this case Plummer softening, can be used as a basis for developing new softening kernels whose gravity biases grow more slowly with increasing softening length. These then provide much more accurate estimates for mode frequencies than Plummer softening does.

Generally, the use of gravitational softening lowers the exponential growth rate of spiral modes. Therefore, strongly softened $N$-body simulations may risk “losing” some of these modes as their growth rates are overtaken by that of e.g. swing amplified noise (Romeo, 1995).

Acknowledgements

We wish to thank the organizers of the workshop “The secular evolution of self-gravitating systems over cosmic ages” at the Institut d’Astrophysique de Paris, May 24-27 2016, where this work was initiated. We are grateful to Prof. Jerry Sellwood for his very helpful insights into setting up and performing $N$-body simulations and to Prof. Scott Tremaine for his helpful suggestions. JBF acknowledges support from Program number HST-HF2-51374 which was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. SDR acknowledges financial support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 721463 to the SUNDIAL ITN network. WD acknowledges support by STFC grant ST/N000757/1.

Appendix A The gravity bias and the suppression factor

A.1 Gravity bias

Analogous to the three-dimensional case discussed in Dehnen (2001), in two dimensions, the expectation value of the gravitational potential is

[TABLE]

with $\Sigma$ the surface density that causes the gravitational potential $V$ through the inter-particle interaction potential $-\psi$. The integral covers the whole surface of the galaxy. The gravity bias can be obtained by rewriting the interaction potential as

[TABLE]

which leads to

[TABLE]

with $\varepsilon\boldsymbol{u}=\boldsymbol{r}-\boldsymbol{r}^{\prime}$. We replace the surface density by its Taylor expansion around the position $\boldsymbol{r}$,

[TABLE]

such that

[TABLE]

If we perform the integration in polar coordinates $(u,\theta)$ and replace $(\boldsymbol{u}\cdot\boldsymbol{\nabla})^{n}=u^{n}(\cos\theta\,\nabla_{\!x}+\sin\theta\,\nabla_{\!y})^{n}$ by its binomial expansion, we have

[TABLE]

The $\theta$ integral vanishes for odd $n$ or $l$, while for even $n$ and $l$

[TABLE]

such that we obtain (using $\Delta\equiv\boldsymbol{\nabla}^{2}$)

[TABLE]

with coefficients $a_{n}$ as given by equation (9).

A.2 Relation to the reduction factor

Because $f(\boldsymbol{r})\equiv r^{-1}-\psi(r)$ is an isotropic function in $\boldsymbol{r}$ space, its moments

[TABLE]

vanish for odd $l$ or $m$. For even $l$ and $m$,

[TABLE]

where $n=l+m$ and we have used relation (61). Let

[TABLE]

be the two-dimensional Fourier transform of $f(\boldsymbol{r})$. Because $f(\boldsymbol{r})$ is isotropic and real-valued, then so is $F(\boldsymbol{k})=F(k)$. From the equation (65)

[TABLE]

and hence the Taylor expansion of $F(\boldsymbol{k})$

[TABLE]

The reduction factor is related to $F(\boldsymbol{k})$ via

[TABLE]

Thus, the coefficients $c_{n}$ of the Taylor series

[TABLE]

are given by

[TABLE]

for even $n$ and if $a_{n}$ is finite. For odd $n$ and/or infinite $a_{n}$, the situation is more complicated. If the kernel $\psi$ has compact support or $\psi\sim r^{-1}$ exponentially fast as $r\to\infty$, then all the coefficients $a_{n}$ are finite, the Taylor series (67) of $F(\boldsymbol{k})$ converges, and relation (70) holds for all $n$, i.e. $R(k)$ is a function of $k^{2}$ only. This is the situation for the $F_{3}$ and $L_{2}$ kernels.

If the kernel does not satisfy the above conditions, but $f(r)\sim r^{-p}$ at $r\to\infty$, then $a_{n}=\infty$ for $n\geq p-2$ and the Taylor series (67) does not converge. However, the series of the non-divergent terms is still useful, only the remainder grows faster than $\varepsilon^{p-2}$. A typical example is Plummer softening, for which $a_{0}=2\pi$, while $a_{n>0}=\infty$ and (see equation 52)

[TABLE]

i.e. $c_{0}=a_{0}/2\pi$ as per relation (70), but $c_{1}=-\tfrac{1}{2}\neq 0$. The problem is that the two-dimensional Fourier transform (71) is not smooth at the origin, but has discontinuous gradient (such that its Taylor series fails), while the one-dimensional function $F(k)$ and hence the reduction factor $R(k)$ are well behaved for all $k\geq 0$. Thus, a divergent $a_{n}$ indicates a non-vanishing $c_{n-1}$ and, conversely, a $c_{n-1}\neq 0$ for even $n$ implies $a_{n}=\infty$.