MOND theories in galpy?

[edan0314]

Is it possible to use MOND theories within Galpy?

[jobovy]

There are currently no MOND models implemented in galpy, but as long as they could be written as a set of forces, they could easily be implemented. I tried implementing a simple example and learned that it is a bit clunky currently, because ideally one would implement this as a 'wrapper', but wrappers currently only really work for Potential instances, not more general Force instances. However, one can just do the wrapping oneself. The following is an example MOND force that uses the 'simple interpolation function' mu(x) = 1/[1+1/x] (the 0.5*(abar+numpy.sqrt(abar**2.+4.*self._a0*abar)) part is the solution for the total acceleration using the MOND equation):

import numpy
from galpy.potential import DissipativeForce, evaluateRforces, evaluatezforces, evaluatephitorques, _isNonAxi
from galpy.util import conversion
class MondForce(DissipativeForce.DissipativeForce):
    def __init__(self,amp=1.,baryon_force=None,a0=0.612,ro=8.,vo=220.):
        DissipativeForce.DissipativeForce.__init__(self,amp=amp,ro=ro,vo=vo)
        self._a0= conversion.parse_force(a0,ro=ro,vo=vo)
        self._baryon_force= baryon_force
        self.isNonAxi= _isNonAxi(self._baryon_force)
    def _calc_abar(self,R,z,t,phi):
        return numpy.sqrt(evaluateRforces(self._baryon_force,R,z,t=t,phi=phi,use_physical=False)**2.
                          +evaluatephitorques(self._baryon_force,R,z,t=t,phi=phi,use_physical=False)**2./R**2.
                          +evaluatezforces(self._baryon_force,R,z,t=t,phi=phi,use_physical=False)**2.)
    def _Rforce(self,R,z,t=0.,phi=0.,v=None):
        abar= self._calc_abar(R,z,t,phi)
        return 0.5*(abar+numpy.sqrt(abar**2.+4.*self._a0*abar))*evaluateRforces(self._baryon_force,R,z,t=t,phi=phi,use_physical=False)/abar
    def _phitorque(self,R,z,t=0.,phi=0.,v=None):
        abar= self._calc_abar(R,z,t,phi)
        return 0.5*(abar+numpy.sqrt(abar**2.+4.*self._a0*abar))*evaluatephitorques(self._baryon_force,R,z,t=t,phi=phi,use_physical=False)/abar
    def _zforce(self,R,z,t=0.,phi=0.,v=None):
        abar= self._calc_abar(R,z,t,phi)
        return 0.5*(abar+numpy.sqrt(abar**2.+4.*self._a0*abar))*evaluatezforces(self._baryon_force,R,z,t=t,phi=phi,use_physical=False)/abar

then you can for example plot the rotation curve of a MWPotential2014-like model that uses MOND rather than DM (MWPotential2014[:2] is just the baryonic part of MWPotential2014):

from galpy.potential import MWPotential2014, plotRotcurve
import matplotlib.pyplot as plt
MondMWPotential2014= MondForce(baryon_force=MWPotential2014[:2],a0=1.*1e-10*u.m/u.s**2)
plotRotcurve(MWPotential2014[:2],ro=8.,vo=220.,label=r'MWPotential2014: baryons')
plotRotcurve(MWPotential2014,ro=8.,vo=220.,label=r'MWPotential2014: baryons+DM',overplot=True)
rs= numpy.linspace(0.,40.,101)*u.kpc
plt.plot(rs,u.Quantity([numpy.sqrt(-r*evaluateRforces(MondMWPotential2014,r,0.,v=[0.,0.,0.],ro=8.,vo=220.,quantity=True)) for r in rs]).to(u.km/u.s),
         label=r'MWPotential2014: baryons+MOND')
plt.ylim(0.,370.)
plt.legend(frameon=False,fontsize=15.)

This gives

This is only a simple example and it makes an assumption about how the 3D MOND force emerges from the 1D MOND equation (i.e., that the 1D MOND equation relates the magnitude of the acceleration to the magnitude of the baryonic acceleration and that the acceleration points in the same direction as the baryonic acceleration). If there is serious interest in using MOND with galpy, we could try re-writing the code to at least make implementing this slightly easier (using wrapper forces, allow use of plotRotcurve rather than computing the circular velocity by hand as in the example above, ...).

[edan0314]

My goal is verify orbital integrations of ultra-faint dwarf galaxies and their uncertainty.
Does this MondMWPotential2014 that you created coul be used inside .integrate() as a potential like MWPotential2014 and other models for an arbitrary orbit?

[jobovy]

Yes, you can use this to do orbit integrations, although it will be quite slow, because it's not using the C back-end (since this new force isn't implemented in C). But this works:

from galpy.orbit import Orbit
o= Orbit()
ts= numpy.linspace(0.,1.,1001)*u.Gyr
o.integrate(ts,MWPotential2014)
o.plot(label=r'MWPotential2014: baryons+DM')
o.integrate(ts,MondMWPotential2014)
o.plot(overplot=True,label=r'MWPotential2014: baryons+MOND')
plt.xlim(7.7,9.5)
plt.ylim(-0.2,0.2)
plt.legend(frameon=False,fontsize=16.)

and gives the orbit of the Sun in both potentials:

[jobovy]

I ran the code above with astropy-units=True in the configuration file. Without that option, the evaluateRforces function in the line that computes the MondMWPotential2014line returns a number, not a Quantity. I've edited the code to work regardless of how galpy is configured and I repeat the updated code here:

from galpy.potential import MWPotential2014, plotRotcurve
import matplotlib.pyplot as plt
MondMWPotential2014= MondForce(baryon_force=MWPotential2014[:2],a0=1.*1e-10*u.m/u.s**2)
plotRotcurve(MWPotential2014[:2],ro=8.,vo=220.,label=r'MWPotential2014: baryons')
plotRotcurve(MWPotential2014,ro=8.,vo=220.,label=r'MWPotential2014: baryons+DM',overplot=True)
rs= numpy.linspace(0.,40.,101)*u.kpc
plt.plot(rs,u.Quantity([numpy.sqrt(-r*evaluateRforces(MondMWPotential2014,r,0.,v=[0.,0.,0.],ro=8.,vo=220.,quantity=True)) for r in rs]).to(u.km/u.s),
         label=r'MWPotential2014: baryons+MOND')
plt.ylim(0.,370.)
plt.legend(frameon=False,fontsize=15.)

The orbit integration code should work as it was.

[edan0314]

Now worked softly!
Thanks!
I'm working on this project for my thesis and I'm waiting my orientor to answer me If we'll use this mond potential on our analysis.
However we are definetly goin to use Galpy
May I mark this post as answerd and when he gives me the feedback I give you a return if we are going to use this mond or not?

[jobovy]

Sounds good! Because this was a topic in the 'General' discussion rather than Q&A, I don't think you can label a post as the answer, but thanks!

[edan0314]

Dear Jo Bovy,
I just speaked with my orientor and he said we can make use of this MoND potential for future parts of our work.
However I still hae some questions about this model.

1. What model of MoND this potential that you created is? MoND has various formulations, and one important fact is cite refences about the potential we are using to describe MoND.
2. It is possible for you to send an article that uses this parametrization for this MoND interpolation function?
3. What acctually the code you created do? New a new user of python, actually I just started using it because of Galpy, So I am not too familiar with the treatment you did to create this potential.
4. You said that we could make it easier to implement, what is actually on your mind about it?

[jobovy]

Hi!

I used the 'simple interpolation function' mu(x) = 1/[1+1/x] in the MOND equation

F_newton = m mu(a/a0) a. (1)

where a is the acceleration, a0 is the MOND acceleration (which I set to a0=1.*1e-10*u.m/u.s**2 above), and F_newton is the Newtonian force (from the baryons). To make this a 3D equation, I did

1. Assume Equation (1) holds between the magnitude of the Newtonian force and the magnitude of the acceleration
2. Solve for the magnitude a using the simple interpolation function: F_newton / m = g_newton (what galpy forces really are) = a/[1+a0/a], --> a^2/g_newton = a + a0 --> a = 0.5*(g_newton+numpy.sqrt(g_newton**2.+4.*a0*g_newton))
3. Then get components of the acceleration by assuming that vec{a} // vec{g_newton}

Note that these are assumptions about how to create a 3D MOND theory that is consistent with the simple 1D MOND theory that explains rotation curves (for which you only need the radial force in the plane of the galaxy).

I think the simple interpolation function was proposed in this paper: https://ui.adsabs.harvard.edu/abs/2005MNRAS.363..603F/abstract

The code I created implemented the MOND forces as a galpy force so it can be used within galpy to integrate orbits. Forces in galpy are either (a) conservative, in which case they are obtained from a potential (these derive from a base Potential class) or (b) not-necessarily-conservative, in which case the forces are directly implemented (these derive from a base DissipativeForce class, because one of their uses is implementing dissipative forces such as dynamical friction. Because this MOND force does not arise from a potential, I implemented it as a sub-class of a DissipativeForce. But basically it's just a way of implementing the 3D force which is then used to integrate orbits (and compute rotation curves etc.).

In terms of improving this code, one main issue is that every time a component of the force is requested, the force is computed from scratch. But during orbit integration, one normally needs all 3 components of the force, so one could just compute the magnitude g_newton once and then get all three components. This often comes up as an issue and I generally speed up the code by storing parts of the force calculation and re-using it if one asks for another component at exactly the same position. For example, see

galpy/galpy/potential/ChandrasekharDynamicalFrictionForce.py

Lines 246 to 252 in 97bbae3

	new_hash = hashlib.md5(
	numpy.array([R, phi, z, v[0], v[1], v[2], t])
	).hexdigest()
	if new_hash != self._force_hash:
	self._calc_force(R, phi, z, v, t)
	return self._cached_force * v[0]

This isn't crucial though, the code I shared above works, it might just be a bit slower than it could be (so if it's fast enough for you, don't worry about it). If it is very slow, then even more speed-ups would be gained by implementing the force in C, but that would be much more complicated.

I hope this answers all of your questions!

[edan0314]

It answerd all the questions! I'll certainly study more this code and the basic MoND theories. All this answers arrived from one question, which I never imagined that you would make a MoND potential. I am still learning many topics in astrophysics, and hope to learn more.
Since our work is to analyse UFD's and the uncertainty in its proper motions, the speed is still good enough for me! Even if we increase the number of UFD's anlysed, it still done in a satisfactory amount of time!
By the way I must thank and congrats you for the awesome tool you developed!
And for answer some questions, not only this, but a few more questions in the past!

[jobovy]

You're welcome, it was an interesting problem to play around with!