¿Existe una forma matemática de calcular la densidad, el volumen, etc. de un planeta simplemente usando la masa como datos?

esquema de algoritmo

Entonces, como dijo L. Dutch, si asumes el volumen o la densidad promedio del planeta, puedes hacer esto. Pero el volumen y la densidad promedio no son valores básicos; ellos mismos dependen del material y la composición de una manera que describiré en breve. Como una palabra de advertencia, voy a suponer que entiendes de cálculo para esto, si no lo haces, la respuesta de L. Dutch es lo mejor que puedes hacer.

Entonces, la dificultad fundamental de este problema es que para conocer la masa del planeta, necesitas conocer la distribución de la densidad. La otra respuesta asume que es constante, pero esta no es una suposición bien fundada: el interior del planeta tiene mucha más presión, lo que conduce a densidades más altas. Entonces, ¿cómo lidiamos con esta dificultad? ¡Por qué, con ecuaciones diferenciales, por supuesto! Salvo cualquier asunto divertido como un planeta que gira rápidamente o uno que todavía está acumulando, la ecuación principal que usaría es la ecuación de equilibrio hidrostático que básicamente dice que la presión que empuja contra una capa delgada de material debe cancelar exactamente la fuerza de gravedad tirando hacia abajo :

dónde$r$es la distancia radial desde el centro del planeta,$P(r)$es la presión a esa profundidad,$\rho(r)$la densidad y$G$la constante gravitatoria.$m(r)$es una función que te da la masa dentro de la esfera de radio$r$, y como tal$m(r) = \int_0^r 4 \pi r^2 \rho(r') dr'$. Pero queremos resolver un sistema de ecuaciones diferenciales, así que derivamos esto para obtener

Ahora, solo queremos integrar este sistema de ecuaciones incorporando su condición de contorno. Sin embargo, espere, no estamos del todo listos: hay tres funciones que queremos determinar:$\rho$,$P$, y$m$, pero sólo dos ecuaciones. ¡Necesitamos otra ecuación para resolver esto!

Esa ecuación viene en forma de la llamada ecuación de estado, que es donde entra la composición planetaria. Resulta que para un material dado, hay una función que relaciona la densidad y la presión llamada ecuación de estado (que denotaremos por$g$):

Debo mencionar en primer lugar que el EOS también relaciona otra variable con estos dos, generalmente la temperatura, la entropía o la energía interna específica. Sin embargo, para determinar los perfiles de densidad, los efectos son bastante pequeños (~ un poco %) por lo que generalmente puede establecer la temperatura en un valor razonable como 5,000 K y extraer g de esa isoterma.

Ahora, las cosas pueden ponerse un poco complicadas aquí: no se garantiza que la ecuación de estado dé lugar a una relación uno a uno entre la densidad y la presión. En otras palabras, no necesariamente podemos escribir$P = f(\rho)$y$\rho = f^{-1}(P)$. Si este es el caso, tenemos que usar un solucionador para lo que se conoce como una ecuación algebraica diferencial . Puede imitar este comportamiento utilizando un paquete ODE normal disponible en la mayoría de los idiomas, pero puede ser bastante lento. Otra opción es el lenguaje Julia, que tiene un gran paquete de ecuaciones diferenciales que puede resolver ecuaciones diferenciales algebraicas . Sin embargo, por lo general para este régimen, al menos podemos escribir$P = f(\rho)$aunque no podemos invertirlo, lo que para este sistema de ecuaciones específico es lo suficientemente bueno como para reescribirlas como:

Ahora, ¿cómo encontramos un EOS para un material determinado? Pues hay varias formas:

1. Lo más común en el campo de la ciencia planetaria es usar una tabla prefabricada generada por alguien que sabe más sobre EOS que tú. SESAME y FPEOS son dos que conozco: este último es de fácil acceso pero, lamentablemente, no tiene muchos elementos importantes (sobre todo hierro, que es importante para los planetas rocosos). El primero es más completo pero hay que solicitar el acceso a través del enlace que di. De cualquier manera, con esta opción, tendrá que pasar algún tiempo aprendiendo otros scripts o escribiendo sus propios scripts para analizar las tablas.
2. Puede hojear la literatura para encontrar un artículo con una fórmula o diagrama que pueda digitalizar , como este .
3. Simplemente puede dar algún tipo de fórmula empírica que crea que parece lo suficientemente razonable. Relaciones de ley de potencia como$P = C\rho^{\gamma}$son a menudo populares por esto.

¡Uf, finalmente tenemos nuestras ecuaciones configuradas! Ahora solo necesitamos usar nuestras condiciones de contorno: si$R$es el radio de tu planeta, esta condición es que$m(R) = M$, dónde$M$es la masa que deseas que tenga el planeta. Pero desafortunadamente, dado que no sabemos qué$R$debería ser, tenemos que adoptar un enfoque iterativo descrito por el siguiente flujo de control:

1. Colocar$m(r = \epsilon) = 0$y$\rho(r = \epsilon) = \rho_{guess}$. Aquí$\epsilon$es un pequeño número positivo pero no del todo$0$porque nuestras ecuaciones son singulares allí. Probablemente pueda hacer un cambio de variable para arreglar esto, pero es menos claro físicamente lo que está pasando y también soy perezoso.
2. Resuelva las dos ecuaciones acopladas anteriores hasta obtener el radio$R$tal que$P(R) = f(\rho(R)) = 0$(ya que no hay presión en la superficie del planeta). Es posible que tenga que usar una pequeña presión positiva distinta de cero por razones numéricas.
3. Controlar$M(R)$.

• Si está lo suficientemente cerca de la masa que desea, ¡ya está! Extraiga los perfiles de densidad y presión de la solución ODE.
• Si es demasiado grande, haz$\rho_{guess}$más pequeño y repita el paso 2
• Si es demasiado pequeño, haga$\rho_{guess}$más grande y repita el paso 2

¡Y eso es todo! Si tengo tiempo más tarde, intentaré agregar un script de ejemplo simple para usar, pero espero que sea suficiente para comenzar.

EDITAR: Imágenes y un código simple.

Así que hice un código Python simple para calcular las masas planetarias, el resultado está a continuación. Debería ser bastante fácil de usar, todo lo que necesita hacer es tener los paquetes apropiados instalados y cambiar Mdesired y la función de ecuación de estado a su gusto:

import scipy.interpolate as interp
import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt
import scipy.integrate as sciint
from math import pi

G = 6.67408e-11   # gravitational constant in [m^3][kg^-1][s^-2]
rho_guess = 1e4   # initial guess for center of planet density in [kg][m^-3]
Mdesired = 6e24          # Desired planet mass in [kg]

# Now we normalize units so that G, M, rho_guess=1.
# This will help numerical stability.
rhostar = rho_guess
Mstar = Mdesired
Tstar = (G * rhostar)**(-1/2)
Lstar = (Mstar / rhostar)**(1/3)
Pstar = Mstar / (Lstar * Tstar**2)
Mdesired = Mdesired / Mstar
rho_guess = rho_guess / rhostar
M0 = 0.0
eps = 1e-4  # start ODE at this value of r to avoid singularity

def dudr(r, u):
    # u[0] = rho, u[1] = m
    drhodr = -(u[1]/(r**2 * dPdrho(u[0])) * u[0])
    dmdr = 4*pi*r**2*u[0]
    return [drhodr, dmdr]

def P(rho):
    """EOS of epsilon iron phase from OSTI 6345571 
    """
    rho0 = 8.430*1e3 / rhostar
    beta0 = 182*1e9 / Pstar
    betaprime0 = 5.0
    eta = rho/rho0

    P = 1.5*beta0*((eta**(7/3) - eta**(5/3))
                   *(1 + (3/4)*(eta**(2/3) - 1)*(betaprime0 - 4)))
    return P

rho_EOS_arr = np.linspace(0, 50, 1000)
dPdrho_arr = np.gradient(P(rho_EOS_arr), rho_EOS_arr)
dPdrho = interp.interp1d(rho_EOS_arr, dPdrho_arr, bounds_error=False, 
                         fill_value=(0.0, 0.0))

def found_surface(r, u):
    """Event that terminates ODE when surface is reached
    """
    Psurf = 1e-7 
    return P(u[0]) - Psurf
found_surface.terminal = True

def M(rho_guess, plot=False, obj = False):
    """Find mass of planet given guess for density at center.
    Can also plot density profiles and give ODE solution object
    """
    sol = sciint.solve_ivp(dudr, (eps, 100), (rho_guess, M0),
                           events = found_surface, 
                           t_eval = np.linspace(eps, 50, 10000))

    if plot:
        plt.figure()
        plt.plot(sol.t*Lstar/1e3, sol.y[0,:]*rhostar/1e3)
        plt.title(r"$\rho$ vs $r$")
        plt.xlabel(r"$r$ $(km)$")
        plt.ylabel(r"$\rho$ $(g/cm^3)$")
        
        plt.figure()
        plt.plot(sol.t*Lstar/1e3, sol.y[1,:]*Mstar/1e3)
        plt.title(r"$M$ vs $r$")
        plt.xlabel(r"$r$ $(km)$")
        plt.ylabel(r"$M$ $(kg)$")
    
    if obj:    
        return sol
    else:
        return sol.y[1,-1]

rhoc = opt.newton(lambda rho: M(rho) - Mdesired, rho_guess, tol = 1e-3, maxiter=200)
                      
sol = M(rhoc, plot=True, obj=True)
print(f"radius is {sol.t[-1]*Lstar/1e3} km")
print(f"mass is {sol.y[1,-1]*Mstar:e} kg")

plt.show()

En este momento, lo tengo configurado para que Mdesiredsea la masa de la Tierra para compararla. Esto es lo que predice para el perfil de densidad de la Tierra . No está tan mal: predice que el radio de la Tierra es de ~5000 km en lugar de ~6400 km. También tiene un perfil de densidad algo similar al del núcleo.

La razón principal de la discrepancia es que este cálculo asume que toda la Tierra es Hierro en el$\epsilon$fase: esto realmente solo es cierto para el núcleo, por lo que describe la corteza y el manto como demasiado densos y el resultado final es que predice el radio de la Tierra. Para arreglar esto, tendrías que modificar la ecuación de estado para que su valor cambie con$m(r)$. Este es un buen ejercicio y si lo entendiste del todo, te recomiendo que lo pruebes.

EDIT 2: Un código más preciso

Así que realmente fui y me disparé como un nerd en este caso. A continuación se muestra un código que desarrollé que, en principio, creo que puede brindarle los perfiles de densidad exactos que ve en la literatura (sin correcciones de temperatura), siempre que tenga las ecuaciones de estado adecuadas. Me di cuenta de que cometí un error en mi última edición al describir el camino a seguir para múltiples EOS, si solo hace que el EOS varíe abruptamente por$m(r)$en los límites materiales, dará resultados no físicos. La razón es que el código asume implícitamente un perfil de densidad continuo, cuando en realidad es el perfil de presión el que es continuo.

Hay un par de formas de lidiar con esto, pero lo que hace mi código es resolver el perfil de densidad de un material a la vez, uniendo las interfaces con las condiciones de contorno correctas. Como entradas, necesitas:

1. Una masa deseada para el planeta en$kg$
2. Una estimación de la densidad en el centro del planeta en$kg/m^3$. Esto no tiene que ser súper preciso: ~ 10000 probablemente funcionará, pero si está trabajando con masas muy altas o bajas, es posible que deba ajustar para que los solucionadores converjan
3. Una colección de materiales que definen la composición de tu planeta. Hay dos listas correspondientes a este elemento: una da la fracción de la masa planetaria total para cada material y otra es una lista de ecuaciones de estado para cada material. He incluido los de perovskita de hierro y silicato, que es una buena aproximación para la Tierra, pero si desea otros materiales, simplemente examine la literatura y encuentre más fórmulas.

¡Y eso es! Todo lo que necesita cambiar son las variables en la sección de entrada, todo lo demás debería hacer el trabajo pesado. Sin más preámbulos, aquí está el código:

import scipy.interpolate as interp
import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt
import scipy.integrate as sciint
from math import pi

#####################################
###### INPUTS #######################
#####################################

def P_silicates(rho):
    """EOS of cold silicate perovskite from Bina 1995
    """
    Ktau0 = 262*1e9 / Pstar
    Kprimetau0 = 4
    rho0 = 4.1 * 1e3 / rhostar
    f = 0.5*((rho/rho0)**(2/3) - 1)
    xi = -(3/4)*(Kprimetau0 - 4)
    Pc = 3*Ktau0*f*(1+2*f)**(5/2)*(1-xi*f)
   
    return Pc

def P_iron(rho):
    """EOS of epsilon iron phase from OSTI 6345571 
    """
    rho0 = 8.430*1e3 / rhostar
    beta0 = 182*1e9 / Pstar
    betaprime0 = 5.0
    eta = rho/rho0
    
    Pc = 1.5*beta0*((eta**(7/3) - eta**(5/3))
                    *(1 + (3/4)*(eta**(2/3) - 1)*(betaprime0 - 4)))
    return Pc

G = 6.67408e-11   # gravitational constant in [m^3][kg^-1][s^-2]
rho0_guess = 1e4   # initial guess for center of planet density in [kg][m^-3]
Mdesired = 6e24          # Desired planet mass in [kg]
mat_Mfracs = [0.3, 0.7]  # Mass fraction for each material type, starting from
                         # core outward
mat_EOSs = [P_iron, P_silicates]  # EOS for each material type


##############################################
######## SOLVER ##############################
##############################################

# Now we normalize units so that G, M, rho_guess=1.
# This will help numerical stability.
msum = sum(mat_Mfracs)
mat_Mfracs = [val/msum for val in mat_Mfracs]
rhostar = rho0_guess
Mstar = Mdesired
Tstar = (G * rhostar)**(-1/2)
Lstar = (Mstar / rhostar)**(1/3)
Pstar = Mstar / (Lstar * Tstar**2)
Mdesired = Mdesired / Mstar
rho0_guess = rho0_guess / rhostar
M0 = 0.0

def dudr(r, u, Mdesired, mat_props):
    # u[0] = rho, u[1] = m
    drhodr = -(u[1]/(r**2 * dPdrho(*u, Mdesired, mat_props)) * u[0])
    dmdr = 4*pi*r**2*u[0]
    return [drhodr, dmdr]

def dPdrho(rho, m, Mdesired, mat_props):
    Pfunc = mat_props["Pfunc"]
    drho = 1e-9
    P1 = Pfunc(rho+drho)
    Pn1 = Pfunc(rho-drho)
    return (P1 - Pn1) / (2*drho)

def end_of_material(r, u, Mdesired, mat_props):
    """Event that terminates ODE when the end of a shell for a given material 
    is reached.
    """
    cum_mass_frac_end = mat_props["cum_mass_frac_end"]
    return u[1]/Mdesired - cum_mass_frac_end
end_of_material.terminal = True    

def found_surface(r, u, Mdesired, mat_props):
    """Event that terminates ODE when surface is reached
    """
    Psurf = 1e-7 
    Pfunc = mat_props["Pfunc"]
    return Pfunc(u[0]) - Psurf
found_surface.terminal = True

def solve_one_mat(rho0, m0, Mdesired, mat_props):
    r0 = mat_props["r0"]
    last_mat = mat_props["last_mat"]
    if last_mat:
        sol = sciint.solve_ivp(
            dudr, (r0, 100), (rho0, m0),
            args = (Mdesired, mat_props),
            events = found_surface, 
            t_eval = np.linspace(r0, 100, 100000),
            tol = 1e-8
            )
    else:
        sol = sciint.solve_ivp(
            dudr, (r0, 100), (rho0, m0),
            args = (Mdesired, mat_props),
            events = (found_surface, end_of_material), 
            t_eval = np.linspace(r0, 100, 100000),
            tol = 1e-8
            )
    return sol

def M(rho0, Mdesired, mat_Mfracs, mat_EOSs, plot=False):
    """Find mass of planet given guess for density at center.
    Can also plot density profiles and give ODE solution object
    """
    cumsum_Mfracs = np.cumsum(mat_Mfracs)
    rho0, m0, r0 = [rho0, 0.0, 1e-7]  # solver starts at 1e-7 radius to avoid
                                      # singularity
    P0 = mat_EOSs[0](rho0)
    r_arr = np.array([])
    rho_arr = np.array([])
    m_arr = np.array([])
    
    for i in range(len(mat_Mfracs)):
        # Build the mat_props dict for a given material, this tells the
        # solver what material we're using
        mat_props = {}
        mat_props["Pfunc"] = mat_EOSs[i]
        mat_props["cum_mass_frac_end"] = cumsum_Mfracs[i]
        mat_props["r0"] = r0
        mat_props["last_mat"] = (i == len(mat_EOSs)-1)
        rho0 = opt.newton(lambda rho: mat_EOSs[i](rho) - P0, rho0, tol = 1e-8)
        
        sol = solve_one_mat(rho0, m0, Mdesired, mat_props)        
        rho0, m0, r0 = [sol.y[0,-1], sol.y[1,-1], sol.t[-1]]
        P0 = mat_EOSs[i](rho0)
        r_arr = np.concatenate((r_arr, sol.t))
        rho_arr = np.concatenate((rho_arr, sol.y[0,:]))
        m_arr = np.concatenate((m_arr, sol.y[1,:]))
        
        # If surface is reached before we can get through all materials,
        # we must end the loop
        if len(sol.t_events[0]) != 0 and not(mat_props["last_mat"]):
            missed_layers = True
            break
        else: 
            missed_layers = False
    
    if plot:
        plt.figure()
        plt.plot(r_arr*Lstar/1e3, rho_arr*rhostar/1e3)
        plt.title(r"$\rho$ vs $r$")
        plt.xlabel(r"$r$ $(km)$")
        plt.ylabel(r"$\rho$ $(g/cm^3)$")
        
        plt.figure()
        plt.plot(r_arr*Lstar/1e3, m_arr*Mstar)
        plt.title(r"$M$ vs $r$")
        plt.xlabel(r"$r$ $(km)$")
        plt.ylabel(r"$M$ $(kg)$")
    
    return (m_arr[-1], r_arr, rho_arr, m_arr, missed_layers)

# Solve for density at planet center that gives the desired planetary mass
rhoc = opt.newton(
    lambda rho0: M(rho0, Mdesired, mat_Mfracs, mat_EOSs)[0] - Mdesired, 
    rho0_guess, tol = 1e-4, maxiter = 200
    )
# Solve for and plot final values                          
R, r_arr, rho_arr, m_arr, missed_layers = M(
    rhoc, Mdesired, mat_Mfracs, mat_EOSs, plot=True
    )
print(f"radius is {r_arr[-1]*Lstar/1e3} km")
print(f"mass is {m_arr[-1]*Mstar:e} kg")

En este momento tengo algunos parámetros de entrada que son más o menos correctos para la Tierra, a saber$M_{desired} = 6\times 10^{24} kg$, con una composición de 30% Hierro y 70% Silicatos. Entonces, ¿cómo lo hace? ¡Bastante bien, en realidad! Aquí está el perfil de densidad que genera:

¡Lo cual es bastante parecido a los tipos de perfiles que verás si haces una búsqueda rápida en Google! También estima el radio de la Tierra en 6250 km, muy cerca del valor real de 6370 km. ¡Espero que esto haya sido útil para ti!

Un par de simplificaciones que he hecho si quieres profundizar más

1. Ignoro la dependencia de la temperatura. En realidad, el EOS relacionará 3 variables y necesitará una tercera ecuación que le diga cómo varía la temperatura con la profundidad. Esta es una ecuación más complicada porque depende del tiempo de acumulación y decaimiento radiactivo y un montón de otras cosas que debe especificar.
2. En realidad, la forma de la EOS no será la EOS de un solo material: a medida que se profundiza, la composición planetaria cambia. Tbh, no estoy seguro de cómo calcula esto a partir de los primeros principios: probablemente haya algún tipo de argumento termodinámico / flotabilidad, pero no estoy seguro. EDITAR: mi segundo código aborda esto cuando hay límites claros entre los tipos de materiales. En realidad, habrá una zona de transición de cierto espesor, pero parece que la mayoría de la literatura ignora esto de todos modos.
3. Si hay otras fuerzas en su planeta, como las fuerzas centrípetas del giro, eso puede cambiar la ecuación de equilibrio de presión.
4. Para objetos súper pesados ​​como estrellas de neutrones, debe tener en cuenta GR a través de la ecuación Tolman-Oppenheimer-Volkoff

Para calcular la densidad$\rho$o el volumen$V$de un cuerpo, puedes usar la relación$\rho=M/V$.

Al ser dos valores desconocidos y una sola ecuación, no puede resolverlo a menos que tenga alguna otra restricción para agregar. Si puede dar, por ejemplo, el radio, puede calcular el volumen como$V=4/3 \pi r^3$y de ahí la densidad.

Otra posible extrapolación es que, si asume que la densidad promedio es la misma que la de la Tierra o cualquier otro planeta que desee considerar, a partir de la densidad puede calcular el volumen.