How to account the effect of the LMC and SMC

[edan0314]

I'm analysing some UFDs and some of then have great impact of the LMC and/or SMC.
How could I take into account the effect of the Milky Way along with the LMC and SMC all at once?
This is the phase space of one of the dwrfs I'm studying, with respect to the Milky Way, asi stands for the Irrgang13I from galpy.potential, eroor

o=Orbit([354.347*units.deg,-63.266*units.deg,55*units.kpc,-0.161*units.mas/units.yr,-1.157*units.mas/units.yr,-36.2*units.km/units.s],radec=True,**get_physical(asi))
o.integrate(ts,asi,method="rk4_c")

This way I can integrate this orbit with respect to the Milky Way,'s potential but I'm not sure how to take into account the effects of LMC and SMC both at the same time, along with the effect of the Milky Way? Specially because the phase space of this dwarf with respect to the LMC/SMC is different in positions and velocities.

[jobovy]

You can do this by using a MovingObjectPotential for the LMC and the SMC. Because galpy isn't an N-body code, you can't account for the full set of interactions, but luckily the hierarchy of the masses is such that the following is a very good approximation:

1. Integrate the orbit of the LMC in the potential of the MW (asi in your case
2. Set up a MovingObjectPotential for the LMC using this LMC orbit. Because the LMC is so massive, you might also want to account for its effect on the MW barycenter
3. Integrate the orbit of the SMC in the combined potential of the MW + the moving LMC (+LMC-effect-on-barycenter)
4. Set up a MovingObjectPotential for the SMC using this SMC orbit
5. Integrate further orbits in the MW + LMC + SMC (+LMC-effect-on-barycenter)

You can look at this section of the documentation for an example of how to do steps 1 and 2. You can get the SMC orbit using Orbit.from_name('SMC'). I think if you use the example from the documentation, you'll see that the SMC is bound to the LMC, that is, orbiting it as you integrate backwards in time. Setting up a MovingObjectPotential for the SMC is similar to the LMC, you'll need a mass and mass model for the SMC (I don't have any recommendation). Because the SMC has low mass, there's no need to take into account its effect on the MW barycenter (although you could).

Happy to provide further feedback as you try to implement this. If you want feedback on code, full examples that I can ran are most useful.

[edan0314]

Dear Jo Bovy,
After some time, I'm ready to implement this idea, actually I just did. However, I'm having some issues:

1. The console returns the following message:

“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)”.

However, as exhibited in the code below, all time intervals I add units in Gyr.

2. When I make the orbital integration for the Sun, the results from the inertial and non-inertial are considerably different. Unlike what is expected from the tutorial in the documentation. I imagine I'm making mistakes, when trying to use the Irrgang13I potential. Since this potential has internal units, I changed some parts from the tutorial.

Here's my code:

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

plt.rcParams['text.usetex'] = True
plt.rc('text', usetex=True)
plt.rc('font', family='serif')

RFsun=get_physical(asi)
RFsun_num=[]
for k, v in RFsun.items():
    RFsun_num.append(float(v))

msol_to_kg=1.9885*(10**30)

"""
######################################################################################################
######################################################################################################
Defining the MW+CDF+moving LMC potential
######################################################################################################
######################################################################################################
"""
#LMC orbit
oLMC= Orbit.from_name('LMC', ro=RFsun_num[0], vo=RFsun_num[1])
cdf= ChandrasekharDynamicalFrictionForce(GMs=1.38**11*units.Msun, rhm=5.*units.kpc, dens=asi, ro=RFsun_num[0], vo=RFsun_num[1])
ts=np.linspace(0,2,10001)*units.Gyr
oLMC.integrate(ts, asi+cdf, method='rk4_c')

#moving potential
lmcpot= HernquistPotential(amp=1.38*10.**11*units.Msun, a=5.*units.kpc/(1.+np.sqrt(2.)))
moving_lmcpot= MovingObjectPotential(oLMC,pot=lmcpot)
#LMC acceleration over the MW
loc_origin= 1e-4
ax= lambda t: evaluateRforces(moving_lmcpot,loc_origin,0.,phi=0.,t=t, ro=RFsun_num[0], vo=RFsun_num[1], use_physical=True)
ay= lambda t: evaluatephitorques(moving_lmcpot,loc_origin,0.,phi=0.,t=t, ro=RFsun_num[0], vo=RFsun_num[1], use_physical=True)/loc_origin
az= lambda t: evaluatezforces(moving_lmcpot,loc_origin,0.,phi=0.,t=t, ro=RFsun_num[0], vo=RFsun_num[1], use_physical=True)
if oLMC.time(use_physical=True)[0] > oLMC.time(use_physical=True)[1]:
    t_intunits= oLMC.time(use_physical=True)[::-1] # need to reverse the order for interp
else:
    t_intunits= oLMC.time(use_physical=True)
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= NonInertialFrameForce(a0=[ax,ay,az])
#final potential
mw_lmc=asi+nip+moving_lmcpot

sunts= np.linspace(0.,1.,10001)*units.Gyr
osun_inertial= 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], ro=RFsun_num[0], vo=RFsun_num[1])
osun_inertial.integrate(sunts, asi,  method='rk4_c')
osun_inertial.plotx(label=r'$\mathrm{Inertial}$')
osun_noninertial= 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], ro=RFsun_num[0], vo=RFsun_num[1])
osun_noninertial.integrate(sunts, mw_lmc, method='rk4_c')
osun_noninertial.plotx(overplot=True, label=r'$\mathrm{Non-inertial}$')
plt.legend(fontsize=18.,loc='upper left',framealpha=0.8)

The value for the LMC mass (1.38×10^(11) Msun) is taken from: https://academic.oup.com/mnras/article/487/2/2685/5491315
And here is the result:

[jobovy]

I think the main issue with this code is that you should actually use use_physical=False everywhere you use use_physical=True. In particular, the accelerations given to the NonInertialFrameForce must be in internal units and be a function of time in internal units, so that's why doing everything in internal units using use_physical=False is necessary (and that should be enough to accomplish this). I think your code runs and gives reasonable results with that change.

[edan0314]

By the way, I just executed the code from the tutorial. It ran smoothly, fast, and the results were as expected.

[edan0314]

I was making some tests with the MW potential, and to implement the SMC I need its half-mass radius.
However, I just can't find any reference to it. I wonder if you have any estimate value (along with a reference, if possible) for this quantity. As you mentioned that is small, compared to the MW or LMC, I'm considering modeling it as a KeplerPotential. What do you think about it?

[jobovy]

I don't have a standard reference for the SMC's mass model, but I found this paper which seems okay: https://ui.adsabs.harvard.edu/abs/2009MNRAS.395..342B/abstract. Their Fig. 1 implies a half-mass radius ~2 kpc. Modeling as a KeplerPotential may be good enough, but if there's a close encounter that always has the possibility of giving a large deflection, so a slightly softened potential might be better, e.g., a HernquistPotential with a=2.*u.kpc/(1+np.sqrt(2)) (so it has a 2 kpc half-mass radius).

[edan0314]

Good morning, I just tested and it worked!
However, now I'm leading with another problem. With the addition of the SMC moving potential.
When I try to integrate SMC orbit in the potential composed of Irrgang13I + Chandrasekhar DF + Non-inertial forces of LMC + Moving LMC, it gives back:

'Physical unit conversion parameters (ro,vo) are not compatible between potentials to be combined'

I tried to use:
o_smc_asi=Orbit.from_name('SMC')
o_smc_asi=Orbit.from_name('SMC', **get_physical(asi))
o_smc_asi=Orbit.from_name('SMC', **get_physical(asi + cdf2_asi + nip_asi + moving_lmcpot_asi))
None of these are working.

I'll send you the code again. It's a little slow, but I imagine it is because I'm using in the same code both the MWPotential2014 and Irrgang13I to make this new potential that adds LMC & SMC influence.

`import numpy as np
import math
import random
import matplotlib.pyplot as plt
from galpy.util.conversion import get_physical
from astropy import units
from galpy.potential.mwpotentials import Irrgang13I as asi
from galpy.potential import MWPotential2014 as mwp, ChandrasekharDynamicalFrictionForce
from galpy.orbit import Orbit
from galpy.potential import HernquistPotential, MovingObjectPotential
from galpy.potential import (evaluateRforces, evaluatephitorques, evaluatezforces)
from galpy.potential import NonInertialFrameForce

plt.rcParams['text.usetex'] = True
plt.rc('text', usetex=True)
plt.rc('font', family='serif')

RFsun=get_physical(asi)
RFsun_num=[]
for k, v in RFsun.items():
RFsun_num.append(float(v))

ts= np.linspace(0,2.1,10001)*units.Gyr
ts_ufd= np.linspace(0,2,10001)*units.Gyr

"""
######################################################################################################
######################################################################################################
MWPotential2014 -> General potential
######################################################################################################
######################################################################################################
"""
#LMC moving potential
o_LMC_mwp= Orbit.from_name('LMC')
mwp[2]*= 1.5 # Don't run this if you've already run it before in the session
cdf_mwp= ChandrasekharDynamicalFrictionForce(GMs=(1.38/2)10.**11.units.Msun,rhm=5units.kpc, dens=mwp)
o_LMC_mwp.integrate(ts,mwp + cdf_mwp, method="rk4_c")
lmcpot= HernquistPotential(amp=1.3810.**11.units.Msun, a=5units.kpc/(1.+np.sqrt(2.)))
moving_lmcpot_mwp= MovingObjectPotential(o_LMC_mwp,pot=lmcpot)

#LMC on MW influence
loc_origin= 1e-4
ax= lambda t: evaluateRforces(moving_lmcpot_mwp,loc_origin,0.,phi=0.,t=t, use_physical=False)
ay= lambda t: evaluatephitorques(moving_lmcpot_mwp,loc_origin,0.,phi=0.,t=t, use_physical=False)/loc_origin
az= lambda t: evaluatezforces(moving_lmcpot_mwp,loc_origin,0.,phi=0.,t=t, use_physical=False)
if o_LMC_mwp.time(use_physical=False)[0] > o_LMC_mwp.time(use_physical=False)[1]:
t_intunits= o_LMC_mwp.time(use_physical=False)[::-1] # need to reverse the order for interp
else:
t_intunits= o_LMC_mwp.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_mwp.time(use_physical=False)[0] > o_LMC_mwp.time(use_physical=False)[1]:
t_intunits= o_LMC_mwp.time(use_physical=False)[::-1] # need to reverse the order for interp
else:
t_intunits= o_LMC_mwp.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_mwp= NonInertialFrameForce(a0=[ax_int,ay_int,az_int])

#SMC moving potential
o_SMC_mwp=Orbit.from_name('SMC')
cdf2_mwp= ChandrasekharDynamicalFrictionForce(GMs=(7/2)10.**9.units.Msun, rhm=2units.kpc, dens=mwp)
o_SMC_mwp.integrate(ts, mwp + cdf2_mwp + nip_LMC_mwp + moving_lmcpot_mwp, method="rk4_c")
smcpot_SMC_mwp= HernquistPotential(amp=710**9units.Msun, a=2units.kpc/(1.+np.sqrt(2.)))
moving_smcpot_mwp= MovingObjectPotential(o_SMC_mwp, pot=smcpot_SMC_mwp)

#final potential
pot_mw=mwp + nip_LMC_mwp + moving_lmcpot_mwp + moving_smcpot_mwp

"""
######################################################################################################
######################################################################################################
ASI -> General potential
######################################################################################################
######################################################################################################
"""
#LMC orbit
o_lmc_asi= Orbit.from_name('LMC', ro=RFsun_num[0], vo=RFsun_num[1])
cdf_asi= ChandrasekharDynamicalFrictionForce(GMs=(1.38/2)10**11units.Msun, rhm=5.*units.kpc, dens=asi, ro=RFsun_num[0], vo=RFsun_num[1])
o_lmc_asi.integrate(ts, asi+cdf_asi, method='rk4_c')

#moving potential
lmcpot= HernquistPotential(amp=1.3810**11units.Msun, a=5*units.kpc/(1.+np.sqrt(2.)))
moving_lmcpot_asi= MovingObjectPotential(o_lmc_asi,pot=lmcpot)
#LMC acceleration over the MW
loc_origin= 1e-4
ax_asi= lambda t: evaluateRforces(moving_lmcpot_asi,loc_origin,0.,phi=0.,t=t, ro=RFsun_num[0], vo=RFsun_num[1], use_physical=False)
ay_asi= lambda t: evaluatephitorques(moving_lmcpot_asi,loc_origin,0.,phi=0.,t=t, ro=RFsun_num[0], vo=RFsun_num[1], use_physical=False)/loc_origin
az_asi= lambda t: evaluatezforces(moving_lmcpot_asi,loc_origin,0.,phi=0.,t=t, ro=RFsun_num[0], vo=RFsun_num[1], use_physical=False)
if o_lmc_asi.time(use_physical=False)[0] > o_lmc_asi.time(use_physical=False)[1]:
t_intunits= o_lmc_asi.time(use_physical=False)[::-1] # need to reverse the order for interp
else:
t_intunits= o_lmc_asi.time(use_physical=False)
ax4int_asi= np.array([ax_asi(t) for t in t_intunits])
ax_int_asi= lambda t: np.interp(t,t_intunits,ax4int_asi)
ay4int_asi= np.array([ay_asi(t) for t in t_intunits])
ay_int_asi= lambda t: np.interp(t,t_intunits,ay4int_asi)
az4int_asi= np.array([az_asi(t) for t in t_intunits])
az_int_asi= lambda t: np.interp(t,t_intunits,az4int_asi)
nip_asi= NonInertialFrameForce(a0=[ax_int_asi,ay_int_asi,az_int_asi])

#SMC moving potential
cdf2_asi= ChandrasekharDynamicalFrictionForce(GMs=(7/2)10.**9.units.Msun, rhm=2units.kpc, dens=asi)
o_smc_asi=Orbit.from_name('SMC', ro=RFsun_num[0], vo=RFsun_num[1])
"""
R_SMC = o_smc_asi.R()
vR_SMC = o_smc_asi.vR()
vT_SMC = o_smc_asi.vT()
z_SMC = o_smc_asi.z()
vz_SMC = o_smc_asi.vz()
phi_SMC = o_smc_asi.phi()(180/math.pi)
o_smc_asi=Orbit([R_SMCunits.kpc, vR_SMCunits.km/units.s, vT_SMCunits.km/units.s, z_SMCunits.kpc, vz_SMCunits.km/units.s, phi_SMCunits.deg], ro=RFsun_num[0], vo=RFsun_num[1])
"""
o_smc_asi.integrate(ts,asi + cdf2_asi + nip_asi + moving_lmcpot_asi, method="rk4_c")
smcpot_asi= HernquistPotential(amp=710**9units.Msun, a=2*units.kpc/(1.+np.sqrt(2.)))
moving_smcpot_asi= MovingObjectPotential(o_smc_asi, pot=smcpot_asi)

#final potential
pot_asi=asi + nip_asi + moving_lmcpot_asi + moving_smcpot_asi

"""
Testing with the Sun
"""
sunts= np.linspace(0.,0.5,10001)units.Gyr
osun_inertial_mwp= Orbit([8.4units.kpc, 22.units.km/units.s, 242units.km/units.s, 0.*units.kpc, 22.*units.km/units.s, 0.units.deg], **get_physical(mwp))
osun_inertial_mwp.integrate(sunts, mwp, method='rk4_c')
osun_noninertial_mwp= Orbit([8.4units.kpc, 22.units.km/units.s, 242units.km/units.s, 0.*units.kpc, 22.*units.km/units.s, 0.*units.deg], **get_physical(pot_mw))
osun_noninertial_mwp.integrate(sunts, pot_mw, method='rk4_c')

osun_inertial_asi= Orbit([8.4*units.kpc, 22.units.km/units.s, 242units.km/units.s, 0.*units.kpc, 22.*units.km/units.s, 0.units.deg], **get_physical(asi))
osun_inertial_asi.integrate(sunts, asi, method='rk4_c')
osun_noninertial_asi= Orbit([8.4units.kpc, 22.units.km/units.s, 242units.km/units.s, 0.*units.kpc, 22.*units.km/units.s, 0.*units.deg], **get_physical(pot_asi))
osun_noninertial_asi.integrate(sunts, pot_asi, method='rk4_c')

osun_inertial_mwp.plotx(color='r', label=r'$\mathrm{Inertial - MWP}$')
osun_noninertial_mwp.plotx(overplot=True, color='b', label=r'$\mathrm{Non-inertial - MWP}$', linestyle='--')
osun_inertial_asi.plotx(overplot=True, color='m', label=r'$\mathrm{Inertial - ASI}$')
osun_noninertial_asi.plotx(overplot=True, color='c', label=r'$\mathrm{Non-inertial - ASI}$')
plt.legend(bbox_to_anchor=(1.005,0.9), loc='upper left', borderaxespad=0)`

[edan0314]

I found the problem!
I forgot to define ro and vo for the Chandrasekhar DF acting on SMC, this was leading to errors with the physical units!
This is the result:

[jobovy]

That's great!

By the way, I forgot to address your question about

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)”

This is not a problem, it just results from the way MovingObjectPotential is implemented triggering this warning, but you can ignore it in this case (it's been on my to-do list for a while to suppress this warning).