handling symbolical and numerical nested integrals

Question

I'm working on reproducing a calculation from this paper https://arxiv.org/abs/1712.03972 for my thesis but I'm struggling with how to properly implement the nested integrals in Python, specifically how to handle them symbolically and numerically.

I want to reproduce the left panel in Figure 1, which requires computing the energy densities $\rho_\phi(T)$ (eq. 3.4) and $\rho_N(T)$ (eq. 3.7).

$$ \rho_\phi(T) = \frac{g_\phi}{2\pi} \int_0^\infty dp \ p^2 \sqrt{p^2 + m_\phi^2} f_\phi(p,T) $$

$$ \rho_N(T) = \frac{1}{\tau_\phi} \left( g_{*s}(T)^{1/3}T \right)^4 \int_T^{T_{cd}} \frac{m_\phi n_\phi(\lambda)}{\left( g_{*s}(\lambda)^{1/3}\lambda\right)^4} \frac{1 +\Delta_{*s}(\lambda)}{H(\lambda)\lambda} d\lambda $$

First, I compute the distribution function $f_\phi(T,p)$ as in eq. 3.2.

$$ f_\phi(p, T)|_{T<T_{cd}} = \left[ \exp\left( \frac{\sqrt{m_\phi^2+p_*^2(T,T_{cd})}}{T_D(T_{cd})} \right) -1 \right]^{-1}\exp\left( -\frac{1}{\tau_\phi} \int_T^{T_{cd}} \frac{m_\phi}{\sqrt{m_\phi^2 + p_*^2(T,\lambda)}} \frac{1+\Delta_{*s}(\lambda)}{H(\lambda)\lambda} d\lambda \right) $$

where $\Delta_{*s} \equiv \frac{T}{3g_{*s}(T)} \frac{d g_{*s}(T)}{dT}$ and $p_*^2(T,\lambda) \equiv p^2\left( \frac{g_{*s}(\lambda)^{1/3} \lambda}{g_{*s}(T)^{1/3}T} \right)^2$.

From my understanding the the integral in $\lambda$ should be solved symbolically because the lower limit is $T$, and I want to have a function of $T$ and $p$ to later integrate numerically or symbolically over $p$ to get $\rho_\phi(T)$ and $n_\phi(T)$ (eq. 3.5)

$$ n_\phi(T) = \frac{g_\phi}{2\pi} \int_0^\infty dp \ p^2 f_\phi(p,T) $$

And then I have to integrate $n_\phi$ in eq. 3.7 to obtain $\rho_N(T)$.

The first problem is given by the degrees of freedom $g_{*s}(T)$, which I take from this file http://www.laine.itp.unibe.ch/eos15/. If I keep $g_{*s}(T)$ as a symbolic function in sympy, it does not solve the integral and returns an Integral object instead. If I instead plug in a numerical value for $T$, the integral gets solved (because $g_{*s}(T)$ is tabulated). I then numerically integrate my $f_\phi$ over $p$ to get $\rho_\phi$.

To generalize for all the temperatures I simply loop over all my values of $T$. However, since I do not obtain the same results of the paper I do not know if this approach is incorrect or of there are other mistakes. Furthermore, with this approach I can not compute $\rho_N(T)$ since I do not have an expression for $n_\phi(T)$ to integrate.

To summarize, is there a more efficient way to handle the symbolic integration involving $g_{∗s}(T)$? How can I compute $\rho_N​(T)$? Is there a more effective approach to integrate these functions both symbolically and numerically?

Any help will be very much appreciated! Thank you!

Here is the code

import numpy as np
from tqdm import tqdm
from scipy.interpolate import interp1d
from sympy import symbols, pi, sqrt, integrate, Function, exp, lambdify, oo
from scipy.integrate import quad

# Reading data
data = np.loadtxt('temperatures_data_SM.dat', skiprows=1)
T_values = data[:, 0]
g_s = data[:, 7]

# necessary to compute dg_*s/dT in Delta_*s
e = data[:, 2]  # energy density / T^4
c = data[:, 4]  # heat capacity / T^3

# Interpolating functions for g_s, c, and e
g_s_func = interp1d(T_values, g_s, kind='cubic', fill_value='extrapolate')
c_func = interp1d(T_values, c, kind='cubic', fill_value='extrapolate')
e_func = interp1d(T_values, e, kind='cubic', fill_value='extrapolate')

# Define symbols for SymPy
T, lam, p = symbols('T lam p')

m_phi = 10.0  # MeV
tau_phi = 20.0  # s 
T_cd = 5.0e3  # MeV
T_D_cd = 5.0e3  # MeV

Grav = 6.70881e-45 # hbar c^4 MeV^-2
hbar = 6.58211899e-22 # MeV s
G8P3 = sqrt(8*pi*Grav/3)/hbar

T_values = T_values[T_values<T_cd]

# Define symbolic functions for the interpolations
g_s_sym = Function('g_s_sym')(T)
c_sym = Function('c_sym')(T)
e_sym = Function('e_sym')(T)

def p_star_sym(lam, T, p):
    return p * (g_s_sym**(1/3) * lam / (g_s_sym**(1/3) * T))

def Delta_s_sym(lam, T):
    return lam / (3 * g_s_sym) * (30 / lam * c_sym - 120 / lam * e_sym)

def Hubble_sym(lam, T):
    rho_SM = (pi**2 / 30) * g_s_sym * lam**4  
    return G8P3 * sqrt(rho_SM)

# Define the function to integrate symbolically
def integral_term_sym(lam, p, T):
    return (m_phi / sqrt(m_phi**2 + p_star_sym(lam, T, p)**2)) * (1 + Delta_s_sym(lam, T)) / (Hubble_sym(lam, T) * lam)

# Define a function that substitutes numerical values after symbolic integration
def evaluate_expression(T_value, inner_int_res):
    g_s_sub = g_s_func(T_value)
    c_sub = c_func(T_value)
    e_sub = e_func(T_value)
    
    return inner_int_res.subs({
        g_s_sym: g_s_sub,
        c_sym: c_sub,
        e_sym: e_sub,
        T: T_value
    })

def tot_func(T_example):

    # Perform the symbolic integration in f_phi(p,T)
    int_res_1 = integrate(integral_term_sym(lam, p, T), (lam, T_example, T_cd))

    # Evaluate the integration result at T = T_example
    mid_res_1 = evaluate_expression(T_example, int_res_1)

    # Compute final f_phi(p, T_example)
    factor = exp(sqrt(m_phi**2 + p_star_sym(T_example, T_cd, p)**2) / T_D_cd) - 1 
    f_phi = exp(- mid_res_1 / tau_phi) / factor

    # Compute  and return rho_phi
    rho_phi = p**2 * sqrt(p**2 + m_phi**2) * f_phi

    integrand_function = lambdify(p, rho_phi, 'numpy')

    integral_p_value, error = quad(integrand_function, 0, np.inf, epsabs= np.inf)

    return float(integral_p_value / (2 * pi**2))


lst = []
for i in tqdm(T_values[T_values<T_cd]):
    rhos = tot_func(i)
    lst.append(rhos)