The density of stars drops off exponentially from the center of the galaxy and can be approximated as follows using a cylindrical coordinate system:

$$ \rho(R, z) = \rho_0 \exp \frac{-R}{h_R} \ \exp \frac{-|z|}{h_z} $$


The radial scale length and height of the thin disk are ~ 3000 pc and ~ 200 pc respectively. At the moment rho_0 will just be taken as 1. We observe at a given l and b of 90 and 15 respectively using the galactic coordinate system (just look at l and b for now, the rest are annotations I added that I will describe shortly). The galactic coordinate system uses spherical coordinates:

enter image description here

We want to know where the density of stars peak along at given LOS (from Earth) per 1 sq degree at a given l and b. There are 41253 sq degrees in the entire sky, so 1/41253 will be the fraction we want to multiply the density equation by. However, we must also consider how the volume element changes by distance along the given LOS, so we need to also multiply rho and 1/41253 by this change in volume. This volume element can be represented by $4 \pi r^2 dr$ (the surface area of the sphere and the small change in the distance along the LOS). We are almost done but need a way to transform between coordinate systems by relating galactic coordinates to distance r and R (as seen in the diagram above). Mathematically, they are related as follows (Eq. 3 in Where Does the Disk Stop and the Halo Begin? Kinematics in a Rotation Field Morrison, Flynn and Freeman 1990 Ast. J. 100 (4) 1191-1222):

$$ r^2 = R_☉^2 + R^2 - 2 R R_☉ \cos b \cos l $$


So, for a given r (distance from star or LOS to sun), l and b, we can find R (the distance from the G.C. to the star or LOS). This is in cylindrical coordinates, so we need to use cos(b)R and sin(b)R to find R and z in the density equation above (note these are different r's). The volume element needs to be the LOS from Earth, which corresponds is the R in the diagram and equation for the paper (for the given l,b, and r). Putting this all together, I have the following code. Please let me know if this attempt seems reasonable, or if not, where I went wrong.

import numpy as np
import matplotlib.pyplot as plt

C=8500 #Distance from center of galaxy to sun (R_sun_

#Galactic Coordinates
b=np.deg2rad(15) #above plane

#Scale heights for the thin disk

nsqdeg_in_sky=41253 #square deg in whole sky

#R: Distance from star to sun

#distance from star to G.C.
#r: distance from star to G.C. kpc

#Spherical Galactic Coordinates to cylindrical coordinates
#z component of R
#r component of R

#density equation describing thin disk

#desnity equation combining volume element and portion of sky
N_thin = (1./nsqdeg_in_sky)*(4*np.pi)*R*R*rho_thin

plt.title("Number of stars per 1 square degree")
plt.xlabel(r'$distance$ [pc]')
plt.ylabel(r'$N_\ast$ [pc$^{-1}$]')

edit to add plot for given l and b above: enter image description here

  • 1
    $\begingroup$ I converted the screenshots of equations to MathJax and mad a few other format changes Please feel free to edit further. $\endgroup$
    – uhoh
    Aug 2, 2020 at 22:44
  • 1
    $\begingroup$ This plot cannot be what you seek if the star has a high elevation over the disk, since the starting density near the sun is near maximal for this$R$, declining with $z$ over the disk. Given the monotonic decline of density with $R$ and $z$ the maximum along a ray will essentially be the closest point to the centre. $\endgroup$ Aug 3, 2020 at 7:10
  • $\begingroup$ @Anders Sandberg You're seeing more stars in your FOV as you look out because of the volume element. That's an important part to consider. It's not just a function of density. $\endgroup$
    – Astroturf
    Aug 3, 2020 at 8:39
  • $\begingroup$ We are only interested in where contribution from the thin disk stars peak. $\endgroup$
    – Astroturf
    Aug 3, 2020 at 9:10

1 Answer 1


If you have a ray from an origin point $\mathbf{x}_0$ with unit direction $|\mathbf{v}|=1$, $\mathbf{x}(s)=[x(s),y(s),z(s)]=\mathbf{x}_0+\mathbf{v}s$ and look at the density times distance $\rho(s)=\rho_0 e^{-R(\mathbf{x})/h_R}e^{-z(\mathbf{x})/h_z}s^2$ the derivative is $\rho'(s) =\rho_0\left\{ [-e^{-R(\mathbf{x})/h_R}R'(\mathbf{x})/h_R]e^{-z(\mathbf{x})/h_z}s^2 + [-e^{-z(\mathbf{x})/h_z}z'(\mathbf{x})/h_z]e^{-R(\mathbf{x})/h_R}s^2 +[2s]e^{-R(\mathbf{x})/h_R}e^{-z(\mathbf{x})/h_z} \right\}$ $=\rho(s)\left\{[-R'(\mathbf{x})/h_R]+[-z'(\mathbf{x})/h_z]+[2s]/s^2\right\}$.

Now, $R'(s) = (x(s)v_x+y(s)v_y)/R(s)$ and $z'(s)=v_z$ so the criterion for a maximum or minium is $0 = -(x(s)v_x+y(s)v_y)/h_RR(s)-v_z/h_z + 2/s$ or $$\frac{2}{s} = \frac{[x_0v_x+y_0v_y] +[v_x^2 + v_y^2]s}{h_R\sqrt{(x_0+v_xs)^2+(y_0+v_ys)^2}}+\frac{v_z}{h_z}$$

I can't see how to simplify it further. I think the right approach here is just to do a numeric solution using golden section search rather than do a lot of algebra that anyway ends with a root finding algorithm. Doing that also allows you to use any stellar density model.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .