Once we add a third dimension, the initial conditions for an isothermal disk get a bit more complex. In addition to balancing the radial components of the gravitational potential with circular velocity, I now need to balance the vertical acceleration via a pressure gradient, i.e. I need to make sure the disk is in hydrostatic balance. Given that the pressure and density can be related for an isothermal gas via the expression , this translates to calculating the vertical density profile at each radial location in the disk.
The gravitational potential
The static potential for the disk will once again be based on the Milky Way. Because it is spherically symmetric, the NFW halo can thus be expressed in exactly the same way as the 2D case, but for reasons that will become clear momentarily, I will express it in cylindrical coordinates, :
Again, , , and is the halo scale length (see previous post for variable values).
The potential of the disk now has a vertical component, and becomes a Miyamoto-Nagai profile:
where kpc is the disk scale length and is the disk scale height. From these potentials, I can calculate the radial acceleration and associated circular velocity at every value of , exactly as I did in the 2D case, as well as the vertical acceleration, , at any point in .
The potential described above gives the circular motion of the gas in the disk. For the disk surface density profile, I will again use an exponential in the radial direction:
where is the gas disk scale length. I thus have the boundary conditions required in order to calculate the vertical density profile at each location in radius, as follows. I start with the equation for hydrostatic equilibrium in the vertical direction,
where is the vertical acceleration due to the combined halo+disk potential (). Given , I can solve for , the density profile as a function of , at a given . First, I substitue in for ,
Then, I rearrange and integrate both sides:
where I have substituted for the integral of . Exponentiating and rewriting the constant of integration, I get
To solve for the constants, we need boundary conditions. We can set by saying that when , where is the central density at a given . Then we have
where is just the value of evaluated at . To solve for , we apply the constraint that the integral of the density profile over must equal the surface density,
Generally, this integral will need to be evaluated numerically.