2

Can someone point me to the numerical methods to solve the Poisson equation for the galaxy rotation curves? I've heard of some tools like PETSc or OpenFOAM. But I'm not sure if these are the right tools for calculating galaxy rotation curves.

Calculation of velocity as $v^2(R) = \frac{ G}{R}M(R)$, where $M(R)$ is the galactic mass till radius $R$, seems too simplified with the assumption of mass distribution having spatial symmetry. So, I'm looking for numerical solution of the Poisson equation, for more accurate results.

Also, I believe, one needs to solve the galaxy rotation curve n 3D to get reasonable results. 2D solutions may be misleading, even though the galaxy is mostly a disc. I'm thinking this because the gravitational field would be in 3D, and restricting the gravitational field to 2D may be erroneous.

Kyle Kanos
  • 29,127
Angela
  • 1,073

1 Answers1

1

When solving the Poisson equation numerically, you can use any programming language that allows for multi-dimensional arrays (C/C++, Fortran, Rust, Python, ....). While I imagine OpenFOAM and PETSc can be used, they are not required.1

You have a relatively simple PDE that you can apply discretization to use in a numerical model. In this way, your PDE would appear as, assuming a cylindrical model,2 \begin{align}\nabla^2\Phi\triangleq & \frac{1}{r_i}\left(r_{i+1/2}\frac{\Phi_{i+1,j,k}-\Phi_{i,j,k}}{dr^2}-r_{i-1/2}\frac{\Phi_{i,j,k}-\Phi_{i-1,j,k}}{dr^2}\right) \\ &\quad+\frac{1}{r_i^2}\frac{\Phi_{i,j+1,k}+\Phi_{i,j-1,k}-2\Phi_{i,j,k}}{d\phi^2}\\ &\quad+\frac{\Phi_{i,j,k+1}+\Phi_{i,j,k-1}-2\Phi_{i,j,k}}{dz^2}=\rho(r_i,\,\phi_j,\,z_k)\tag{1} \end{align} This would require a 3D array for $\Phi$ and some tricks for handling $r\simeq0$. Note that under axial symmetry, the middle term vanishes as, $$\Phi_{i,j+1,k}=\Phi_{i,j-1,k}=\Phi_{i,j,k},$$ which would require only a 2D array to solve since $j$ index wouldn't matter anymore. In either case, you need to solve Equation (1) for $\Phi_{i,j,k}$, which can then be used in a relaxation method (or any of the "See also" links in this Wikipedia entry) to find the potential to the given density profile.

Once you have that solution, you can use the same discretization to determine the velocity profile, $$v(r_i)^2\simeq r_i\frac{\Phi_{i+1,j,k}-\Phi_{i-1,j,k}}{2\cdot dr}.$$ where the $j$ and $k$ indices on $\Phi$ would either be used to average over all heights or for a single slice.


1. My proverbial \$0.02 would be to not use those modules as knowing how to numerically solve PDEs and ODEs is a good tool to have in your possession.

2. For spherical systems, the radial term is of the form $r^{-2}\partial_r\left(r^2\partial_r\Phi\right)$; most sources I've seen prefer writing it as the split operator, $\partial^2_r\Phi+\frac{2}{r}\partial_r\Phi$, and deriving the finite difference stencil from that. It might be easier to just use Cartesian coordinates and compute the spherical/cylindrical/etc coordinates from the positions when computing the mass density.

Kyle Kanos
  • 29,127