Logarithmic potential energy

[adrn]

I was playing with the logarithmic potential in Agama and comparing to my implementation in Gala but found an inconsistency. The definition in Agama states that the potential corresponds to (1/2) \sigma^2 \ln[ r_c^2 + x^2 + (y/p)^2 + (z/q)^2 ], but it looks like the parameter it accepts isv0. Is that the same as \sigma? If so, I would expect this:

import agama
agama_pot = agama.Potential(
    type="logarithmic",
    v0=0.1,
    scaleradius=2.5,
    axisRatioY=1.,
    axisRatioZ=1.,
)
print(agama_pot.potential([5., 0., 0.]))
# 0.007210030031772178

to return the same as

import numpy as np
0.5 * 0.1**2 * np.log(2.5**2 + 5**2)
# 0.017210096880912056

Am I misunderstanding the input parameters?

[adrn]

BTW: I'm assuming, based on this code in the potential factory that what is input to agama.Potential(..., scaleradius) actually becomes the coreradius in the defined logarithmic potential..

[eugvas]

hmm, when I run the first snippet of your code, it prints out 0.0172100968809 – same as the analytic expression..
Did I say that it uses "sigma" anywhere? if so, this must be a mistake; sigma (the velocity dispersion in a self-consistent isothermal model) is \sqrt{2} v0, but the parameters passed to the potential constructor are indeed v0 and scaleRadius (= core radius).

[adrn]

Oh, sorry, I forgot something in the snippet -- I had also used:

import astropy.units as u
agama.setUnits(length=u.kpc, time=u.Myr, mass=u.Msun)

I agree -- if I don't do the setUnits, I get what you see. But not after setUnits.

[adrn]

Did I say that it uses "sigma" anywhere?

The variable names and comment reference sigma here:

Agama/src/potential_analytic.h

Line 92 in 3052049

\f$ \Phi(r) = (1/2) \sigma^2 \ln[ r_c^2 + x^2 + (y/p)^2 + (z/q)^2 ] \f$. */

[eugvas]

ah, I see.
the problem is actually in the definition of the logarithmic potential:
$\Phi = 0.5 v_0^2 \ln(r_c^2+r^2)$
clearly, a dimensional quantity under the logarithm paves a pathway to nasty surprises! it should be understood (implicitly) as taking $\ln([r_c^2+r^2]/L^2)$, where $L$ is the length unit. Since only the potential difference between two points is important in practice, the choice of $L$ rarely matters – unless, of course, you just want to compute the value of the potential for some reason. By calling setUnits(...), you change the value of $L$ expressed in some fixed internal code units, so the potential gets a constant offset. You may verify that the potential difference between any two points remains unchanged.
I don't know what is the best way to proceed here: one could redefine $\Phi = 0.5 v_0^2 \ln(1 + r^2/r_c^2)$, but this would not work if $r_c=0$ (which is allowed). And since the unit conversion occurs at a higher level than the potential evaluation itself, there isn't really any way to make the potential instance aware of the units choice (quite the opposite – the whole code is designed to be agnostic to units, and this works fine everywhere except this peculiar case).

Thanks for catching the wrong usage of sigma – it is indeed intended to be v0, I will update the code shortly. The docs (reference.pdf) are correct though.

[adrn]

Ah, of course - that makes sense.

[eugvas]

ok, it turned out to be not such a big deal to add an extra dimensional parameter (length unit) to the Logarithmic potential, and restore the invariance of its zero-point w.r.t. the choice of units (available in the latest commit).