Según entiendo tu pregunta, y después de ver tu simulación, parece que realmente sólo tener materia oscura (DM) en su simulación. Un código de N-cuerpos simula la DM, es decir, las partículas que interactúan sólo por gravitación. Si quieres ir un paso más allá, debes incluir gas es decir, partículas que interactúan tanto por gravitación como por hidrodinámica.
Esto no es trivial, pero la idea básica es que cuando las partículas de gas se acercan lo suficiente, acumulan presión, son frenadas por las fuerzas de viscosidad, experimentan turbulencias, se calientan y pierden energía por radiación, etc. Si una región se vuelve lo suficientemente densa, puede formar estrellas, que a su vez vuelven a calentar el gas.
Cuando un código está basado en partículas (en lugar de en cuadrículas), esto se llama hidrodinámica de partículas suavizada (SPH).
Sin embargo, creo que puedes empezar haciendo algo más sencillo. Puede que me equivoque, pero viendo tu simulación, me parece que estás calculando las fuerzas gravitacionales entre la partícula $i$ y la partícula $j$ a través de $$ F_{ij} = \frac{m_i m_j}{r_{ij}^2}, \qquad\qquad(\mathrm{correct \, in \, real \, life}) $$ lo que, en principio, es correcto. Sin embargo, en las simulaciones numéricas, debido a los tamaños de paso finitos y a la precisión de los puntos flotantes, las partículas se acercan fácilmente, lo que da lugar a fuerzas divergentes, que tienden a arrojar partículas con velocidades anormalmente altas. Por esta razón, se utiliza el llamado longitud de ablandamiento gravitacional $\varepsilon$ , de manera que la fuerza entre las partículas se convierte en $$ F_{ij} = \frac{m_i m_j}{r_{ij}^2 + \varepsilon^2} \qquad\qquad(\mathrm{correct \, in \, simulations}). $$ El valor óptimo de $\varepsilon$ es una cuestión de elección. Puede probar diferentes valores, pero debería ser del orden de (o un poco menor que) la distancia media entre las partículas en las regiones que quiere resolver. Sus pasos de tiempo también deben adaptarse a esto, de modo que las partículas no viajen en promedio más de $\varepsilon/2$ por paso de tiempo (por ejemplo Rodionov y Sotnikova 2005 ). Por ejemplo, pruebe y empiece por poner $\varepsilon = L/100$ , donde $L$ es la longitud de su caja.
La inclusión de este parámetro debería hacer que las partículas de DM se agrupen más. Sin embargo, si quieres patrones espirales hermosos, tienes que incluir gas.
