Miyamoto–Nagai Potential

The Miyamoto–Nagai potential (or Miyamoto–Nagai disk model) is an analytical, axisymmetric three-dimensional gravitational potential–density pair introduced by Japanese astronomers Masanori Miyamoto and Ryuzaburo Nagai in 1975 to model the gravitational field and three-dimensional mass distribution of thickened stellar disks in spiral and lenticular galaxies. Before 1975, analytical models of galactic potentials were restricted either to spherically symmetric systems—such as the Plummer model—or to razor-thin two-dimensional disks with zero vertical thickness, such as the Kuzmin (1956) and Toomre (1963) disk models, whereas computing the gravitational potential of a three-dimensional thickened disk required numerical evaluation of Bessel-function integrals. By replacing the vertical coordinate |z| in Kuzmin's potential with the smoothed algebraic combination (z² + b²)^1/2, Miyamoto and Nagai produced a closed-form algebraic potential Φ(R, z) whose ratio of vertical scale height b to radial scale length a continuously interpolates from an infinitely flat disk (b/a = 0) to a spherical Plummer star cluster (a/b = 0).

In cylindrical galactocentric coordinates (R, ϕ, z), where R = (x² + y²)^1/2 is the projected distance in the galactic midplane and z is the vertical height perpendicular to the disk plane, a Miyamoto–Nagai model of total mass M, radial scale length a ≥ 0, and vertical scale height b ≥ 0 is defined by the gravitational potential:

Φ(R, z) = −GM / { R² + [a + (z² + b²)^1/2]² }^1/2

Every partial derivative of Φ(R, z) involves only elementary powers and square roots, making radial force F_R = −∂Φ/∂R and vertical restoring force K_z = −∂Φ/∂z rapid to evaluate during numerical orbit integrations:

∂Φ/∂R = GMR / { R² + [a + (z² + b²)^1/2]² }^3/2

∂Φ/∂z = *GMz*[a + (z² + b²)^1/2] / { (z² + b²)^1/2 [R² + (a + (z² + b²)^1/2)²]^3/2 }

Substituting Φ(R, z) into the axisymmetric Poisson equation in cylindrical coordinates, ∇²Φ = (1/R) ∂/∂R(R ∂Φ/∂R) + ∂²Φ/∂z² = 4πGρ(R, z), yields the closed-form three-dimensional Miyamoto–Nagai mass density distribution:

ρ(R, z) = (b²M/4π) { aR² + [a + 3(z² + b²)^1/2][a + (z² + b²)^1/2]² } / { [R² + (a + (z² + b²)^1/2)²]^5/2 (z² + b²)^3/2 }

Because the numerator of ρ(R, z) is positive for all R ≥ 0, z ∈ ℝ whenever a ≥ 0 and b > 0, the Miyamoto–Nagai density is non-negative everywhere in space (ρ(R, z) > 0) and integrates over all space to the finite total mass M. At large distances (r = (R² + z²)^1/2 ≫ a + b), the density falls off as ρ ∝ R²/r⁵ near the midplane (z ≈ 0) and as ρ ∝ |z|⁻⁵ along the rotation axis (R = 0).

The algebraic structure of Φ(R, z) unifies two classic models of gravitational potential theory as endpoints of the flatness ratio b/a:

1. The Spherical Plummer Limit (a = 0, b > 0): When the radial disk scale vanishes (a = 0), R² + [a + (z² + b²)^1/2]² simplifies to r² + b², where r is the three-dimensional spherical radius. The potential and density reduce to the spherical Plummer model with scale radius b: Φ(r) = −GM/(r² + b²)^1/2 and ρ(r) = (3M/4πb³)(1 + r²/b²)^−5/2.
2. The Razor-Thin Kuzmin–Toomre Limit (b → 0, a > 0): When the vertical thickness vanishes (b → 0), (z² + b²)^1/2 becomes |z|, recovering Grigori Kuzmin's (1956) and Alar Toomre's (1963) Model 1 potential, Φ_Kuzmin(R, z) = −GM/[R² + (a + |z|)²]^1/2. Above the disk (z > 0), Φ_Kuzmin(R, z) matches the gravitational potential of a point mass M located on the symmetry axis at a fictitious mirror point (R = 0, z = −a) below the plane; across the midplane z = 0, the vertical derivative ∂Φ/∂z jumps discontinuously by 4πGΣ_Kuzmin(R), where Σ_Kuzmin(R) = (Ma/2π)(R² + a²)^−3/2 is the surface density of the razor-thin disk (which is also a separable axisymmetric Stäckel potential).

In the equatorial plane (z = 0) of a Miyamoto–Nagai disk, the gravitational potential simplifies to Φ(R, 0) = −GM/[R² + (a + b)²]^1/2, which depends on a and b only through the combined scale length s ≡ a + b:

• Circular Velocity Curve v_c(R): Balancing centrifugal force against the radial potential gradient gives v_c²(R) = GMR²/[R² + (a + b)²]^3/2, which peaks at R_max = 2^1/2(a + b) with v_c,max ≈ 0.6204[GM/(a + b)]^1/2.
• Radial and Vertical Epicyclic Frequencies κ(R) and ν_z(R): The circular angular speed is Ω(R) = (GM)^1/2/[R² + (a + b)²]^3/4, the radial epicyclic frequency governing Lindblad resonances, Oort constants, and the Toomre stability criterion is κ²(R) = *GM*[R² + 4(a + b)²]/[R² + (a + b)²]^5/2, and the vertical oscillation frequency is ν_z²(R) = [(a + b)/b]Ω²(R).

In composite models of the Milky Way—such as the Paczyński (1990), Allen–Santillán (1991), and Bovy (2015) potentials used to integrate globular cluster orbits, tidal streams, and asymmetric drift kinematics via the Jeans equations and collisionless Boltzmann equation—a Miyamoto–Nagai stellar disk (M_disk ≈ 6.8 × 10¹⁰ M_☉, a ≈ 3.0 kpc, b ≈ 0.28 kpc) or a three-component superposition approximating an exp(−R/R_d) radial disk profile (Smith et al., 2015) is paired with a central Hernquist profile bulge and an outer Navarro–Frenk–White or Einasto profile dark-matter halo.

Why does the midplane rotation curve of a Miyamoto–Nagai disk depend only on a + b?

When z = 0 is substituted into the Miyamoto–Nagai potential, the inner term [a + (z² + b²)^1/2] collapses to (a + b). Because the radial force −∂Φ/∂R in the midplane involves derivatives with respect to R at fixed z = 0, a spherical Plummer model with a = 0, b = 4 kpc and a thin disk with a = 3.7 kpc, b = 0.3 kpc share the same midplane rotation curve v_c(R), even though their vertical restoring forces K_z differ substantially.

Can the density of a Miyamoto–Nagai disk ever become negative?

For a single Miyamoto–Nagai potential with M > 0, a ≥ 0, and b > 0, the three-dimensional mass density ρ(R, z) is positive everywhere in space. When b/a < 0.2, the isodensity contours near the galactic center develop a pinched or "peanut" shape along the R-axis because the vertical gradient steepens near z = 0.

How is the phase-space distribution function constructed for a Miyamoto–Nagai disk?

In an axisymmetric potential Φ(R, z), energy E and vertical angular momentum L_z = Rv_ϕ are conserved integrals of motion. Miyamoto and Nagai (1975) showed that ρ(R, z) can be written as a function of R and Φ(R, z), allowing the two-integral distribution function f(E, L_z) to be recovered via axisymmetric generalizations of the Eddington inversion formula.

1. Miyamoto, M., & Nagai, R. (1975). "Three-dimensional models for the distribution of mass in galaxies." Publications of the Astronomical Society of Japan, 27(4), 533–543. https://doi.org/10.1093/pasj/27.4.533
2. Toomre, A. (1963). "On the Distribution of Matter Within Flattened Galaxies." The Astrophysical Journal, 138, 385–392. https://doi.org/10.1086/147653
3. Paczyński, B. (1990). "A test of the galactic origin of gamma-ray bursts." The Astrophysical Journal, 348, 485–494. https://doi.org/10.1086/168257
4. Smith, R., Flynn, C., Candlish, G. N., Fellhauer, M., & Gibson, B. K. (2015). "Modelling the gravitational potential of thick and thin discs." Monthly Notices of the Royal Astronomical Society, 448(3), 2934–2940. https://doi.org/10.1093/mnras/stv228
5. Bovy, J. (2015). "galpy: A Python Library for Galactic Dynamics." The Astrophysical Journal Supplement Series, 216(2), 29. https://doi.org/10.1088/0067-0049/216/2/29