Mixing MoND + LMC & SMC moving potential

[edan0314]

Dear Jo Bovy,
I'm the author of both of these questions:

(1) - MoND theories in Galpy: #570
(2) - How to account the effect of LMC and SMC: #594

I'd like to mixing those two questions. I implemented the Effects of LMC & SMC for different potentials, but now I'd like to include also a simulation of MoND in my analysis, along with the effects of LMC & SMC. As we discussed in (2) this week, I implemented the effects of those satellites for other potentials different from the MWPotential2014. I tryed to do the same with the potential you developed in (1).
The ChandrasekharDinamicalFriction (CDF) used to calculate the orbit of LMC does not accept the MoND potential in (1). To bypass this issue, I used only the baryon part of the MWPotential2014 to evaluate the dens required to calculate CDF. Which is reasonable to me, since CDF is an emergent force which arises from the movement of a massive body M across a population of stars. This generates an anisotropy of the total gravitational field exerted on M. Since MoND does not have dark matter, all the dynamical friction comes from the bulge and disk. Can you provide feedback about it?
Now the actual problem: Everything runs normally and smooth for the evaluation of the LMC moving potential. But in the evaluation of the SMC, the simulation breaks down. To be more precise in the integration of the SMC orbit. It gives the following warning:

"ValueError: Found time value not in the integration time domain".

Which I don't know what it means. However the same procedure is adopted for other potentials, and nothing bad happens.
The code follows:

import numpy as np
from galpy.potential import MovingObjectPotential, NonInertialFrameForce, HernquistPotential, ChandrasekharDynamicalFrictionForce, DissipativeForce, evaluateRforces, evaluatezforces, evaluatephitorques, _isNonAxi
from galpy.util import conversion
from galpy.potential import MWPotential2014
import matplotlib.pyplot as plt
import matplotlib
from galpy.util.conversion import get_physical
from galpy.orbit import Orbit
from galpy.potential.mwpotentials import Irrgang13I as asi
from astropy import units

MWPotential2014[2]*= 1.5
ts= np.linspace(0, 2, 10001)*units.Gyr

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 np.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+np.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+np.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+np.sqrt(abar**2.+4.*self._a0*abar))*evaluatezforces(self._baryon_force, R, z, t=t, phi=phi, use_physical=False)/abar
mond= MondForce(baryon_force=MWPotential2014[:2], a0=1.2*1e-10*units.m/units.s**2)

o_LMC_mond= Orbit.from_name('LMC')
cdf_mond= ChandrasekharDynamicalFrictionForce(GMs=(1.88/2)*10.**11.*units.Msun,rhm=5*units.kpc, dens=MWPotential2014[:2])
o_LMC_mond.integrate(ts,mond + cdf_mond, method="odeint")
lmcpot= HernquistPotential(amp=1.88*10.**11.*units.Msun, a=5*units.kpc/(1.+np.sqrt(2.)))
moving_lmcpot_mond= MovingObjectPotential(o_LMC_mond,pot=lmcpot)
#LMC on MW influence
loc_origin= 1e-4
ax= lambda t: evaluateRforces(moving_lmcpot_mond,loc_origin,0.,phi=0.,t=t, use_physical=False)
ay= lambda t: evaluatephitorques(moving_lmcpot_mond,loc_origin,0.,phi=0.,t=t, use_physical=False)/loc_origin
az= lambda t: evaluatezforces(moving_lmcpot_mond,loc_origin,0.,phi=0.,t=t, use_physical=False)
if o_LMC_mond.time(use_physical=False)[0] > o_LMC_mond.time(use_physical=False)[1]:
    t_intunits= o_LMC_mond.time(use_physical=False)[::-1] # need to reverse the order for interp
else:
    t_intunits= o_LMC_mond.time(use_physical=False)
ax4int= np.array([ax(t) for t in t_intunits])
ax_int= lambda t: np.interp(t,t_intunits,ax4int)
ay4int= np.array([ay(t) for t in t_intunits])
ay_int= lambda t: np.interp(t,t_intunits,ay4int)
az4int= np.array([az(t) for t in t_intunits])
az_int= lambda t: np.interp(t,t_intunits,az4int)
if o_LMC_mond.time(use_physical=False)[0] > o_LMC_mond.time(use_physical=False)[1]:
    t_intunits= o_LMC_mond.time(use_physical=False)[::-1] # need to reverse the order for interp
else:
    t_intunits= o_LMC_mond.time(use_physical=False)
ax4int= np.array([ax(t) for t in t_intunits])
ax_int= lambda t: np.interp(t, t_intunits, ax4int)
ay4int= np.array([ay(t) for t in t_intunits])
ay_int= lambda t: np.interp(t, t_intunits, ay4int)
az4int= np.array([az(t) for t in t_intunits])
az_int= lambda t: np.interp(t, t_intunits, az4int)
nip_LMC_mond= NonInertialFrameForce(a0=[ax_int, ay_int, az_int])
#SMC moving potential
o_SMC_mond=Orbit.from_name('SMC')
cdf2_mond= ChandrasekharDynamicalFrictionForce(GMs=(7/2)*10.**9.*units.Msun, rhm=2*units.kpc, dens=MWPotential2014[:2])
o_SMC_mond.integrate(ts, mond + cdf2_mond + nip_LMC_mond + moving_lmcpot_mond, method="odeint")
smcpot_SMC_mond= HernquistPotential(amp=7*10**9*units.Msun, a=2*units.kpc/(1.+np.sqrt(2.)))
moving_smcpot_mond= MovingObjectPotential(o_SMC_mond, pot=smcpot_SMC_mond)
#Final potential
pot_mond=mond + nip_LMC_mond + moving_lmcpot_mond + moving_smcpot_mond

rs= np.linspace(0.,50.,1001)*units.kpc
plt.plot(rs,units.Quantity([np.sqrt(-r*evaluateRforces(MWPotential2014,r,0.,v=[0.,0.,0.],ro=8.,vo=220.,quantity=True)) for r in rs]).to(units.km/units.s), label=r'MW2014',color='r')
plt.plot(rs,units.Quantity([np.sqrt(-r*evaluateRforces(asi,r,0.,v=[0.,0.,0.],ro=8.,vo=220.,quantity=True)) for r in rs]).to(units.km/units.s), label=r'ASI',color='lime')
plt.plot(rs,units.Quantity([np.sqrt(-r*evaluateRforces(mond,r,0.,v=[0.,0.,0.],ro=8.,vo=220.,quantity=True)) for r in rs]).to(units.km/units.s), label=r'MoND',color='slategray')
plt.ylim(0.,300.)
plt.legend(bbox_to_anchor=(1.05,1), loc='upper left', borderaxespad=0)

"""
sunts= np.linspace(0.,0.5,10001)*units.Gyr
osun_inertial_mond= Orbit([8.4*units.kpc, 22.*units.km/units.s, 242*units.km/units.s, 0.*units.kpc, 22.*units.km/units.s, 0.*units.deg], **get_physical(mond))
osun_inertial_mond.integrate(sunts, mond,  method='rk4_c')
osun_noninertial_mond= Orbit([8.4*units.kpc, 22.*units.km/units.s, 242*units.km/units.s, 0.*units.kpc, 22.*units.km/units.s, 0.*units.deg], **get_physical(pot_mond))
osun_noninertial_mond.integrate(sunts, pot_mond, method='rk4_c')

osun_inertial_mond.plotx(color='r', label=r'$\text{Inertial}$')
osun_noninertial_mond.plotx(overplot=True, color='b', label=r'$\text{Non-inertial}$', linestyle='--')
plt.legend(bbox_to_anchor=(1.005,0.9), loc='upper left', borderaxespad=0)
"""

In advance, thanks for the help!

[jobovy]

To address your two questions:

1. Yes, in your MOND model, the dynamical friction can only come from the disk and the bulge, because there is no dark matter. However, because the dynamical friction depends on the local density and the LMC/SMC are very far from the disk and the bulge, in practice the dynamical friction force will be very small and I think it would be better to ignore it in this case (because it will just slow down your code).
2. I think this error occurs because the odeint integrator tries to evaluate the force from the LMC just outside the integration-time interval as part of its dynamic time-stepping, but the LMC force is only defined within the integration-time interval. This probably wouldn't occur if you used leapfrog instead, because it doesn't dynamically change the time step, but it will be much slower. You could also try integrating the LMC orbit for a bit longer than the SMC time integration so the evaluation just outside the SMC time-integration interval is within the time-integration interval of the LMC. I hope that makes sense! I'm trying to test it, but the code is running very slow, so I haven't been able to confirm this, but I'm pretty sure this is the issue.

[edan0314]

I followed your recommendation about ignoring the CDF. This surely made the code run faster.
About the time integration, I followed your recommendation about increasing the time integration of LMC. The code is not breaking down no more, And I could evaluate the plots for the inertial and non-inertial potential for MoND. It is something reasonable. This is the result:

However, now galpy shows the following warnings:

lmc_orb ok
galpyWarning: You specified integration times as a Quantity, but are evaluating at times not specified as a Quantity; assuming that time given is in natural (internal) units (multiply time by unit to get output at physical time)
mw_nip ok
smc_orb ok
C:\Users\ericd\anaconda3\Lib\site-packages\galpy\potential\PowerSphericalPotentialwCutoff.py:101: RuntimeWarning: invalid value encountered in scalar divide
return -self._mass(r) * R / r**3.0

C:\Users\ericd\anaconda3\Lib\site-packages\galpy\potential\TwoPowerSphericalPotential.py:891: RuntimeWarning: divide by zero encountered in scalar divide
1.0 / Rz / (self.a + sqrtRz)

C:\Users\ericd\anaconda3\Lib\site-packages\galpy\potential\TwoPowerSphericalPotential.py:892: RuntimeWarning: invalid value encountered in scalar divide

• numpy.log(1.0 + sqrtRz / self.a) / sqrtRz / Rz

C:\Users\ericd\anaconda3\Lib\site-packages\galpy\potential\SCFPotential.py:603: RuntimeWarning: invalid value encountered in divide
dr_dR = numpy.divide(R, r)

C:\Users\ericd\anaconda3\Lib\site-packages\galpy\potential\SCFPotential.py:604: RuntimeWarning: invalid value encountered in divide
dtheta_dR = numpy.divide(z, r**2)

C:\Users\ericd\anaconda3\Lib\site-packages\galpy\potential\SCFPotential.py:499: RuntimeWarning: divide by zero encountered in power

• (l * (a + r) * numpy.power(r, -1) - (2 * l + 1))

C:\Users\ericd\anaconda3\Lib\site-packages\galpy\potential\SCFPotential.py:499: RuntimeWarning: invalid value encountered in multiply

• (l * (a + r) * numpy.power(r, -1) - (2 * l + 1))

c:\users\ericd\desktop\lsst\python\mond\incertezas_mond.py:26: RuntimeWarning: invalid value encountered in scalar divide
+evaluatephitorques(self._baryon_force,R,z,t=t,phi=phi,use_physical=False)2./R2.

C:\Users\ericd\anaconda3\Lib\site-packages\galpy\potential\PowerSphericalPotentialwCutoff.py:105: RuntimeWarning: invalid value encountered in scalar divide
return -self._mass(r) * z / r**3.0

galpyWarning: Cannot use C integration because some potentials are not implemented in C (using odeint instead)
Error in callback <function _draw_all_if_interactive at

I don't know if these warnings could generate errors. So I'll wait for your feedback. As soon as possible as you answer, I'll add this potential in the main file and start the orbital integrations I want to analyze.

By the way, I did this code late this night, and I was sleepy, haha. Looking more tenaciously, and following your instructions, I verified I did the evaluation of the non-inertial acceleration of LMC on MW twice! This was also making the code slower. Already corrected this in my test file!

[jobovy]

I think these are all fine, just warnings about the numerics that one quite often sees in this. The galpyWarning: Cannot use C integration because some potentials are not implemented in C (using odeint instead) warning occurs because your MOND potential is not implemented in C, so that's fine. Plot looks okay as well.