<p>I think your recreation is essentially correct. For example, if you look at Figure 1 in that paper, you can see that the potential goes from positive to negative as a function of the azimuthal angle ("phase").</p>

<p>What you're missing is that the potential and density functions they define are <em>perturbations</em>, which are intended to be <em>added to</em> an axisymmetric galaxy model. The idea is that the axisymmetric disk is modulated by their perturbation, so that the total density is less than average (but not less than zero!) where their perturbation is <span class="math-container">$&lt; 0$</span> and greater than average where their perturbation is <span class="math-container">$&gt; 0$</span>. (Similarly, the total potential is <span class="math-container">$&lt; 0$</span> everywhere; in the regions where their perturbation is positive, the total potential becomes less negative, but never <span class="math-container">$&gt; 0$</span>.)</p>

<p>As they point out (pp.4-5), "The impressions invoked by the density distributions of Figures 3 and 5 can be somewhat misleading. These densities must be regarded as perturbations to an azimuthally uniform stellar disk with the same vertical scale height."</p>

<p>And: "In Figure 10, a disk component with the same radial dropoff and scale height as the perturbation, and just sufficient amplitude to make the net density everywhere positive, has been added to the perturbation density." and "In Figures 12 and 13, various amounts of perturbation density are shown added to a representative full stellar disk."</p>

CC BY-SA 4.0
2020-02-04T22:22:56.360Z
Okay, this is quite a relief, though it's a temporary one. I means that I need to generate a non-modulated density profile and then calculate it's potential field and add that to this, which means I need to actually think a bit. ;-) Thanks for the speedy answer!
uhoh
CC BY-SA 4.0
2020-07-01T13:46:13.223Z
@uhoh Hi - is there a reason you've not marked this answer as accepted? Just noticed it has been a while with no activity on this question. Perhaps you require more information?
Lio Elbammalf
CC BY-SA 4.0
2020-07-01T13:54:54.880Z
Yes, some combination of inattention and procrastination; I'd planned on accepting and then asking a new follow-up question in a few days then lost track. Thanks for the reminder!
uhoh
CC BY-SA 4.0
2020-02-04T16:57:44.270Z
2020-02-04T16:57:44.270Z
Peter Erwin
<p>Thinking about <a href="https://space.stackexchange.com/q/41189/12102">this question</a> I wanted to start with a rough model of the average gravitational potential of the Milky way. I ran across D. P. Cox and G. C. Gomez 2002 <a href="https://arxiv.org/abs/astro-ph/0207635" rel="nofollow noreferrer">Analytical Expressions for Spiral Arm Gravitational Potential and Density</a> which I think I can understand at least enough to script it. They started with an analytical density distribution and approximated its potential with an analytical form, later they solved for the exact density which generates it, but I am not using that.</p>

<p>I calculated the first example discussed in Section 4 and at least at first glance it seems to agree with their figures. At large distances the potential tends to zero which is good, but the modulation is both positive and negative! This is true both from my script and in their figures.</p>

<p>Without a repulsive force, I don't think the potential can be positive. What am I missing?</p>

<p><strong>Note 1:</strong> I've plotted for <span class="math-container">$z=0$</span></p>

<p><strong>Note 2:</strong> I'm using kg meters and seconds for units, so the plotted potential is in m^2/s^2.</p>

<p><a href="https://i.sstatic.net/yEoyg.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/yEoyg.png" alt="gravitational potential of a spiral galaxy"></a></p>

<pre><code>import numpy as np
import matplotlib.pyplot as plt

def PHI(r, phi, z):
    term_1 = -4 * pi * G * H * rho_0
    term_2 = np.exp(-(r-r_0)/Rs)
    gamma  = N * (phi - phi_0 - np.log(r/r_0)/np.tan(alpha))
    K      = n * N / (r * np.sin(alpha))
    KH     = K * H
    beta   = KH * (1 + 0.4*KH)
    D      = (1 + KH + 0.3*KH**2) / (1 + 0.3*KH)
    term_3 = ((C/(K*D)) * np.cos(n*gamma)) * (np.cosh(K*z/beta))**-beta # sech is just 1/cosh

    return term_1 * term_2 * (term_3.sum(axis=0)) # sum over n

G      = 6.67430E-11 # m^3 / kg s^2
parsec = 3.0857E+16 # meters
mH     = 1.007825 * 1.660539E-27 # kg
pi     = np.pi

N     = 2                        # number of arms
alpha = 15 * pi/180.          # pitch angle
Rs    = 7000 * parsec            # radial dropoff
rho_0 = 1E+06 * (14./11) * mH    # midplane arm density 
r_0   = 8000 * parsec            # at fiducial radius
H     = 180 * parsec             # scale height of perturbation

C = np.array([8/(3*pi), 0.5, 8/(15*pi)])[:, None, None]
n = np.array([1, 2, 3])[:, None, None]

# plot it
hw = 30000 * parsec
x = np.linspace(-hw, hw, 200)
X, Y = np.meshgrid(x, x)
r = np.sqrt(X**2 + Y**2)
phi = np.arctan2(Y, X)
z = 0.

phi_0 = 0.

potential = PHI(r, phi, z)

if True:
    plt.figure()
    plt.imshow(potential)
    plt.colorbar()
    plt.gca().axes.xaxis.set_ticklabels([])
    plt.gca().axes.yaxis.set_ticklabels([])
    plt.title('+/- 30 kpc')
    plt.show()
</code></pre>

CC BY-SA 4.0
2020-02-04T13:49:08.800Z
related: [Have parameters been determined for the Fourier components of a Milky Way spiral density-potential model?](https://astronomy.stackexchange.com/q/19795/7982)
uhoh
<galaxy><gravity><galactic-dynamics>
Can the gravitational potential in a spiral galaxy be positive?
uhoh
# Can the gravitational potential in a spiral galaxy be positive?

(Asked by: username_0 on 2020-02-04)

Thinking about this question I wanted to start with a rough model of the average gravitational potential of the Milky way. I ran across D. P. Cox and G. C. Gomez 2002 Analytical Expressions for Spiral Arm Gravitational Potential and Density which I think I can understand at least enough to script it. They started with an analytical density distribution and approximated its potential with an analytical form, later they solved for the exact density which generates it, but I am not using that.
I calculated the first example discussed in Section 4 and at least at first glance it seems to agree with their figures. At large distances the potential tends to zero which is good, but the modulation is both positive and negative! This is true both from my script and in their figures.
Without a repulsive force, I don't think the potential can be positive. What am I missing?
Note 1: I've plotted for $z=0$
Note 2: I'm using kg meters and seconds for units, so the plotted potential is in m^2/s^2.

import numpy as np
import matplotlib.pyplot as plt

def PHI(r, phi, z):
    term_1 = -4 * pi * G * H * rho_0
    term_2 = np.exp(-(r-r_0)/Rs)
    gamma  = N * (phi - phi_0 - np.log(r/r_0)/np.tan(alpha))
    K      = n * N / (r * np.sin(alpha))
    KH     = K * H
    beta   = KH * (1 + 0.4*KH)
    D      = (1 + KH + 0.3*KH**2) / (1 + 0.3*KH)
    term_3 = ((C/(K*D)) * np.cos(n*gamma)) * (np.cosh(K*z/beta))**-beta # sech is just 1/cosh

    return term_1 * term_2 * (term_3.sum(axis=0)) # sum over n

G      = 6.67430E-11 # m^3 / kg s^2
parsec = 3.0857E+16 # meters
mH     = 1.007825 * 1.660539E-27 # kg
pi     = np.pi

N     = 2                        # number of arms
alpha = 15 * pi/180.          # pitch angle
Rs    = 7000 * parsec            # radial dropoff
rho_0 = 1E+06 * (14./11) * mH    # midplane arm density 
r_0   = 8000 * parsec            # at fiducial radius
H     = 180 * parsec             # scale height of perturbation

C = np.array([8/(3*pi), 0.5, 8/(15*pi)])[:, None, None]
n = np.array([1, 2, 3])[:, None, None]

# plot it
hw = 30000 * parsec
x = np.linspace(-hw, hw, 200)
X, Y = np.meshgrid(x, x)
r = np.sqrt(X**2 + Y**2)
phi = np.arctan2(Y, X)
z = 0.

phi_0 = 0.

potential = PHI(r, phi, z)

if True:
    plt.figure()
    plt.imshow(potential)
    plt.colorbar()
    plt.gca().axes.xaxis.set_ticklabels([])
    plt.gca().axes.yaxis.set_ticklabels([])
    plt.title('+/- 30 kpc')
    plt.show()

**Tags:** galaxy, gravity, galactic-dynamics

> **Comments:**
> **username_0**: related: [Have parameters been determined for the Fourier components of a Milky Way spiral density-potential model?](https://astronomy.stackexchange.com/q/19795/7982)  

---

## 1 answers

### Answer by username_1 on 2020-02-04. Score: 4 | Accepted answer

I think your recreation is essentially correct. For example, if you look at Figure 1 in that paper, you can see that the potential goes from positive to negative as a function of the azimuthal angle ("phase").

What you're missing is that the potential and density functions they define are *perturbations*, which are intended to be *added to* an axisymmetric galaxy model. The idea is that the axisymmetric disk is modulated by their perturbation, so that the total density is less than average (but not less than zero!) where their perturbation is $< 0$ and greater than average where their perturbation is $> 0$. (Similarly, the total potential is $< 0$ everywhere; in the regions where their perturbation is positive, the total potential becomes less negative, but never $> 0$.)

As they point out (pp.4-5), "The impressions invoked by the density distributions of Figures 3 and 5 can be somewhat misleading. These densities must be regarded as perturbations to an azimuthally uniform stellar disk with the same vertical scale height."

And: "In Figure 10, a disk component with the same radial dropoff and scale height as the perturbation, and just sufficient amplitude to make the net density everywhere positive, has been added to the perturbation density." and "In Figures 12 and 13, various amounts of perturbation density are shown added to a representative full stellar disk."

> **Comments:**
> **username_0**: Okay, this is quite a relief, though it's a temporary one. I means that I need to generate a non-modulated density profile and then calculate it's potential field and add that to this, which means I need to actually think a bit. ;-) Thanks for the speedy answer!  
> **username_2**: @username_0 Hi - is there a reason you've not marked this answer as accepted? Just noticed it has been a while with no activity on this question. Perhaps you require more information?  
> **username_0**: Yes, some combination of inattention and procrastination; I'd planned on accepting and then asking a new follow-up question in a few days then lost track. Thanks for the reminder!  

---

