Can the gravitational potential in a spiral galaxy be positive?

Question

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

Answer 1

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."