artificial density oscillations in the outskirt of Spheroid potential

[syrte]

Hi Eugene,

I noticed that the density profile of agama spheroid potential sometimes shows oscillations in the outskirt, see e.g., the blue curve below. It is generated by the following code.

pot1 = agama.Potential(
    type='Spheroid',
    densityNorm=1, scaleRadius=1,
    alpha=1, beta=3, gamma=1,
    outerCutoffRadius=4, cutoffStrength=2,
)

In the beginning, I was surprised, because I thought the density profile is analytical. Later, I realized it is because of the limited grid number of multipole decomposition (right?).

For potential directly generated by specifying a density function with type=Multipole, we may use keywords rmin, rmax, gridSizeR etc to specify the desired precision (see the orange curve with gridSizeR=2000). But these keywords seem to not work for the Spheroid potential object, as I have tried. It is a bit wired since Spheroid is implemented by Multipole as well (right?).

While I understand people usually do not care about the very low-density region, it might be good to allow users to control the precision of Spheroid (and maybe other default available potentials) when strange behavior appears.

[syrte]

I think it is my misuse of the potential object. If I really care about the precision of density, I should use a Density object rather than a Potential object.

After some more exploration, I found the default gridSizeR=25 of the Multipole potential is indeed precise enough for calculating potential in most cases. The problem is the density.

When one constructs a Multipole potential pot with a density function (either user defined or agama analytical density), one might naively expect pot.density to give exactly the same result as the original density function. But it is actually not due to the Multipole expansion approximation. It might be good to clarify this in the document (I am sorry if I have overlooked it!)

[eugvas]

hi Zhaozhou, You inferred correctly that the same type of profile used to construct a Density or a Potential object do not always produce identical results. More specifically, models which have analytically expressed potentials, such as Logarithmic, MiyamotoNagai, Ferrers and Dehnen (*), should behave identically in both cases. OTOH some density models (Spheroid, Sersic, Nuker and Disk) do not have analytic potential counterparts, and are automatically converted to a Multipole potential, or in the case of Disk models, to a sum of Multipole and DiskAnsatz (the latter is a rather technical auxiliary object that cannot be instantiated directly). The Multipole potential is capable of reproducing the original density profile accurately enough in the region where it is high, but in the exponentially declining tail it usually cannot do a good job (unless the grid is dense enough, and even then it will probably fail, just a bit later). So if you need an accurate density as well as potential, perhaps it may be worth creating two objects (Density and Potential) with the same parameters (I remind that one can create a single dictionary object with named parameters and pass it to the constructor of both classes). Now there is another twist of story, which produces the counter-intuitive behaviour you noticed. When you create a potential as

p = agama.Potential(type='Spheroid', ...)

it is automatically converted to a Multipole with the default grid parameters that are beyond user control (namely, GridSizeR=50, lmax=32, mmax=6, and automatically selected rmin, rmax). Thus providing these parameters explicitly has no effect. The reason is that one can provide several density components in a single call to Potential constructor (either as several dict objects, or listing them in separate sections of an ini file), and they will be merged into one Multipole potential, plus whatever many DiskAnsatz auxiliary objects); thus the grid parameters provided for any one of these components are ignored. The primary motivation for this behaviour is computational efficiency. OTOH if you first create a Density object and then a Multipole potential from it, as

d = agama.Density(type='Spheroid', ...)
p = agama.Potential(type='Multipole', density=d, gridSizeR=..., rmin=..., rmax=..., lmax=...)

or even merging the two calls into a single line

p = agama.Potential(type='Multipole', density='Spheroid', gridSizeR=..., rmin=..., rmax=..., lmax=...)

then you have explicit control over the grid parameters. The first variant also allows you to use the Density object separately if you need exact density evaluation; in the second variant it is created internally and temporarily but is not user-accessible. In principle, you can even use this approach with multiple density components, creating a single Multipole potential from all of them, while still retaining control over grid parameters:

d1 = agama.Density(type='Spheroid', ...)
d2 = agama.Density(type='Sersic', ...)
p = agama.Potential(type='Multipole', density=agama.Density(d1, d2), gridSizeR=..., rmin=..., rmax=..., lmax=...)

But this cannot be used with Disk density models (or rather, they would not be separated into DiskAnsatz and Multipole-represented residual, and as a result, the entire disk potential would be poorly approximated by just a Multipole). To summarize, different way of creating objects lead to efficient or convenient usage, and may sometimes be combined to achieve both goals, but this flexibility sometimes goes at odds with intuitive expectations. I'll think more about how to explain it better in the docs (as of now, this is all covered but perhaps not highlighted clearly enough).

() NB1: the Dehnen model has "analytic" potential which is actually calculated using 1d numerical integration, and as a result, might be less accurate than a Multipole approximation in practice (and certainly more costly). Thus it is advisable to replace type='Dehnen' with 'Spheroid' in most cases, the latter by default reduces to Dehnen gamma=1 if not explicitly specified otherwise. NB2: Plummer and NFW models have simple analytic potentials in the spherical case, but when created with axisRatioY/Z different from unity, they are automatically* converted to corresponding Spheroid models and then treated as usual (i.e. forwarded to a Multipole potential if necessary). This is a sort of syntactic sugar, but perhaps produces counter-intuitive behaviour that does more damage than benefit - I might decide to eliminate it and demand explicit specification of the Spheroid density type in these cases.

[syrte]

Hi Eugene, Thank you very much for the very detailed explanation and instruction! Now I understand better the rationale and how to make it work. agama is really powerful and flexible!