DDH model EOS inference pipeline

This section outlines the implementation of the density-dependent Hartree (DDH) model within an inference framework. It discusses how to incorporate various likelihood calculations such as perturbative Quantum Chromodynamics (pQCD), Chiral Effective Field Theory (chiEFT), constraints from the maximum mass of neutron stars, and nuclear saturation properties.

Advantages of the DDH Approach

The DDH model provides a significant advantage over models with constant couplings by accounting for quantal fluctuations of the baryon fields even in the ground state, providing a more realistic depiction of nuclear matter under extreme conditions (refer to Furnstahl et al., 1994).

Literature Overview

Over the years, various authors have proposed different functional forms of density-dependent coupling schemes for the \(\sigma\), \(\omega\), and \(\rho\) mesons:

  1. Typel 1999 Model (Typel & Wolter, 1999):

    • The coupling constants \(\Gamma_i(\rho)\) for \(i = \sigma, \omega\) are parametrized as:

      \[\Gamma_i(\rho) = \Gamma_i(\rho_{\mathrm{sat}}) f_i(x)\]

    • The function \(f_i(x)\) is defined by:

      \[f_i(x) = a_i \frac{1 + b_i (x + d_i)^2}{1 + c_i (x + d_i)^2}\]

    • For the \(\rho\) meson:

      \[\Gamma_\rho(\rho) = \Gamma_\rho(\rho_{\mathrm{sat}}) \exp[-a_\rho (x - 1)]\]

  2. Malik 2023 Model (Malik, 2023):

    • This model extends the coupling parametrization to:

      \[\Gamma_M(\rho) = \Gamma_{M,0} h_M(x), \quad x = \rho / \rho_0\]

    • The functions are defined as:

      • For isoscalar (\(\sigma\) and \(\omega\)):

        \[h_M(x) = \exp[-(x^{a_M} - 1)]\]

      • For isovector (\(\rho\)):

        \[h_\rho(x) = \exp[-a_\rho (x - 1)]\]

  3. Gogelein 2008 Model (Gogelein, 2008; Utilized by Char 2023):

    • The density dependence for meson couplings is expressed as:

      \[\Gamma_i(n_B) = a_i + (b_i + d_i x^3) e^{-c_i x}, \quad x = n_B / n_0\]

    • This model allows flexible control across different density regimes, ideal for Bayesian inference applications.

These models are crucial in the study of dense nuclear matter, providing the groundwork for robust statistical analyses and theoretical predictions within the framework of nuclear astrophysics.

[ ]:

from time import perf_counter

import ultranest
import ultranest.stepsampler
import InferenceWorkflow.Likelihood as likelihood
from InferenceWorkflow.DDH_Likelihood import Likelihood as ddh_likelihood
import InferenceWorkflow.prior as prior
import TOVsolver.main as tov
import EOSgenerators.RMF_DDH as DDH

import pandas as pd
import numpy as np
import math
import matplotlib.pyplot as plt

from TOVsolver.unit import g_cm_3, dyn_cm_2, km, Msun, MeV, fm
from TOVsolver.maximum_central_density import maximum_central_density
from TOVsolver.constant import oneoverfm_MeV,G,c
from scipy.interpolate import CubicSpline, InterpolatedUnivariateSpline, interp1d

import warnings
warnings.filterwarnings('ignore')

Set up prior ranges (DDH model ([Malik, 2023] https://doi.org/10.3847/1538-4357/ac5d3c ))

[ ]:

parameters = ["as_", "av", "ar", "gs0", "gv0", "grho0", "rho0"]

def prior_transform(cube):
    params    = cube.copy()
    params[0] = prior.flat_prior(0,0.3 ,cube[0])
    params[1] = prior.flat_prior(0,0.3 ,cube[1])
    params[2] = prior.flat_prior(0,1.3 ,cube[2])
    params[3] = prior.flat_prior(6.5,13.5 ,cube[3])
    params[4] = prior.flat_prior(7.5,14.5 ,cube[4])
    params[5] = prior.flat_prior(5,16 ,cube[5])
    params[6] = prior.normal_Prior(0.153, 0.005,cube[6])

    return params

Setup Likelihood (DDH model ([Malik, 2023] https://doi.org/10.3847/1538-4357/ac5d3c ))

likelihood_transform Function

The likelihood_transform function computes the total log-likelihood for a given set of parameters that define neutron star and nuclear matter properties. It uses a combination of nuclear saturation properties, maximum mass constraints, chiEFT constraints, and pQCD constraints.

Parameters

  • para (array): An array of parameters that define the coupling constants and other relevant physical properties for the model.

Returns

  • float: The total log-likelihood of the parameter set, which combines contributions from nuclear saturation properties, neutron star mass constraints, chiEFT constraints, and pQCD constraints. If the resulting log-likelihood is -inf, it is converted to -1e20 for computational stability.

Function Details

1. Define Meson Masses

The function begins by defining the meson masses in femtometer (fm) units:

  • m_sig: Sigma meson mass (550 MeV)

  • m_w: Omega meson mass (783 MeV)

  • m_rho: Rho meson mass (763 MeV)

  • rho0: Saturation density, taken as the last element in the para array.

2. Density-Dependent Coupling for \(\sigma\), \(\omega\), and \(\rho\)

The density-dependent couplings are calculated using the DDH.Function method, which applies the Malik22 model type for density-dependent coupling in this example. To use a different model type, you would need to modify the 'prior_transform' function and specify the desired model type within 'likelihood_transform', as each model type may require a different number of parameters. The calculated coupling values from para are then assembled into a single array, theta, which is used for subsequent calculations.

3. Nuclear Matter Properties (NMP) Log-Likelihood

The nuclear matter properties (NMP) log-likelihood is computed by initializing a DDH_Log_Likelihood object, DDH_like, using theta.

  • The log-likelihood for NMP is computed using DDH_like.compute_NMP().

4. Neutron Star Maximum Mass Log-Likelihood

The log-likelihood contribution based on neutron star maximum mass is calculated using DDH_like.compute_MR(). This function returns:

  • 0 if the maximum mass (M_max) is above 2 solar masses.

  • -inf if M_max is below 2 solar masses.

The step function used for this calculation is smooth, similar to a Fermi-Dirac distribution, allowing for a gradual transition.

5. chiEFT Constraints

The function then applies chiEFT constraints on pure neutron matter (PNM) using:

  • theta1, which is theta with an additional element 0.0 for pure neutron matter.

  • DDH.compute_eos_alpha(theta1), which generates the equation of state (EOS) for PNM.

The chiEFT likelihood is computed for PNM EOS using log_EFT = likelihood.chiEFT_PNM(...). Units PNM EOS needs to be converted to MeV.fm^-3.

6. pQCD Constraints

Finally, the pQCD constraints are applied on the beta equilibrium EOS using DDH.compute_eos(...):

  • EOS units needs to be converted to GeV.fm^-3 as required by the pQCD likelihood.

  • log_pQCD is computed using likelihood.ln_pQCD(...), applied at a density of 8 * 0.16 fm^-3 to ensure the neutron star maximum mass stays below this density.

7. Total Log-Likelihood Calculation

The total log-likelihood log_tot is calculated by summing:

  • log_NMP, log_MR, log_EFT, and log_pQCD.

If log_tot is -inf, it is set to -1e20 to avoid numerical issues.

Example Usage

```python para = [0.086372, 0.054065, 0.509147, 9.180364, 10.981329, 7.652728, 0.150] log_likelihood = likelihood_transform(para) print(log_likelihood)

[ ]:

def likelihood_transform(para):
    """
    Compute the log-likelihood of a given parameter set for neutron star and nuclear matter properties.

    Parameters:
    para (array): Array of parameters defining the couplings and other physical properties.

    Returns:
    float: The total log-likelihood of the parameter set.
    """

    # Define the meson masses in fm scale
    m_sig = 550 / oneoverfm_MeV
    m_w = 783 / oneoverfm_MeV
    m_rho = 763 / oneoverfm_MeV
    rho0 = para[-1]

    # Calculate the density-dependent coupling factors
    gsf, gwf, grf, dgsf, dgwf, dgrf = DDH.Function(type='Malik22', couplings=para)
    theta = np.array([m_sig, m_w, m_rho, gsf, gwf, grf, dgsf, dgwf, dgrf, rho0])


    # Define log-likelihood for Nuclear Matter Properties (NMP) at nuclear saturation properties
    # Default Gaussian constraints: self.NMP_Base = {"e0":[-16, 0.2, 'g'], "k0":[230, 40, 'g'], "jsym0":[32.5, 1.8, 'g']}
    # You may update these values using the update_base and update_quantities methods in the DDH_Log_Likelihood class.
    DDH_like = ddh_likelihood(theta)

    #1. This block for compute log-likelihood for Nuclear Matter Properties (NMP)
    log_NMP = DDH_like.compute_NMP()

    #2. Compute log-likelihood for NS maximum mass, returning 0 if M_max ≥ 2 and -∞ if M_max < 2
    # The step function is smooth, like a Fermi-Dirac distribution, providing a smooth transition
    log_MR = DDH_like.compute_MR()

    #3. Apply chiEFT constraints to pure neutron matter EOS
    theta1 = np.concatenate((theta, np.array([0.0])))  # Alpha = 0 for pure neutron matter
    EOS_pnm = np.array(DDH.compute_eos_alpha(theta1))

    #4. chiEFT likelihood input requires MeV.fm^-3 units for both energy (e) and pressure (p)
    EOS_pnm[1] *= oneoverfm_MeV
    EOS_pnm[2] *= oneoverfm_MeV

    log_EFT = likelihood.chiEFT_PNM(EOS_pnm, type="Super Gaussian", contraint_quantity="p", enlargement=1.0)  # 100% enlargement for energy in PNM

    #5. Apply pQCD constraints on beta equilibrium EOS
    EOS_beta = np.array(DDH.compute_eos([0],[0],theta))
    EOS_beta[1] *= oneoverfm_MeV / 1000  # pQCD likelihood input requires GeV.fm^-3 units for both e and p
    EOS_beta[2] *= oneoverfm_MeV / 1000

    log_pQCD = likelihood.ln_pQCD(EOS_beta, [7 * 0.16])  # pQCD applied at 7 * 0.16 fm^-3, ensuring NS max mass is below this density


    # Total log-likelihood
    log_tot = log_NMP + log_MR + log_EFT + log_pQCD

    if log_tot == -np.inf:
        log_tot = -1e301

    return log_tot

Test likelihood_transform

The pQCD contribution uses deterministic midpoint quadrature, so repeated evaluations of the same parameters return the same value. This cell times one likelihood evaluation at a reference parameter point. The DDH EOS and Tolman–Oppenheimer–Volkoff (TOV) solves dominate the runtime. Evaluation times vary across the prior and with hardware and software versions. The optional NumbaMinpack dependency accelerates EOSgenerators.fastRMF_EoS and is unused by this DDH calculation.

[ ]:

para = np.array([0.086372, 0.054065, 0.509147, 9.180364, 10.981329, 3.826364 * 2, 0.150])
start = perf_counter()
log_likelihood = likelihood_transform(para)
elapsed = perf_counter() - start
assert np.isfinite(log_likelihood)
print(f"Reference log-likelihood: {log_likelihood:.6f} ({elapsed:.2f} s for one evaluation)")

Inference

Quick sampler check (default)

UltraNest evaluates the prior transform and the full DDH likelihood at one prior draw when it initializes. Together with the reference evaluation above, the default notebook makes two likelihood calls. Allow seconds to a few minutes on a single CPU, depending on the time shown above and the cost of the prior draw. For example, a serial test on an Apple Silicon Mac with Python 3.12, SciPy 1.14 and UltraNest 4.3 took about 19 seconds including imports. Five additional prior draws took 0.8–3.4 seconds each. This check validates the sampler interface. The optional run below generates posterior samples.

[ ]:

start = perf_counter()
ultranest.ReactiveNestedSampler(
    parameters, likelihood_transform, prior_transform, log_dir=None, num_test_samples=1
)
print(f"UltraNest callback check passed in {perf_counter() - start:.2f} s.")

Full nested sampling run (optional)

The reference settings use at least 2,000 live points, 14 slice-sampling steps for the seven parameters, and a budget of 500,000 likelihood calls. On one CPU, a rough estimate is the number of likelihood calls multiplied by the typical time per call. The timings above give two local measurements. Time several prior draws for a more representative estimate before submitting a full run.

Typical time per call

2,000 initial calls

500,000-call budget

1 second

33 minutes

5.8 days

72 seconds (reported by a reviewer)

40 hours

417 days

These figures are serial estimates that exclude sampler overhead. UltraNest can converge before using the full call budget. It may also stop at the call limit before achieving the requested precision. Inspect the convergence diagnostics and effective sample size before using the posterior or evidence.

For a long run, use a batch job on a compute node and keep output/ddh on persistent storage. Set run_full_inference = True below to start. UltraNest saves its progress in that directory and resumes from it when you rerun with the same likelihood, prior, and parameter order. Use a new directory when changing the model or prior. To estimate throughput, start with a small pilot run in a separate directory. UltraNest also supports MPI for running on multiple cores when the notebook workflow is run as a Python script.

[ ]:

run_full_inference = False

if run_full_inference:
    ddh_sampler = ultranest.ReactiveNestedSampler(
        parameters,
        likelihood_transform,
        prior_transform,
        log_dir="output/ddh",
        resume="resume",
    )
    ddh_sampler.stepsampler = ultranest.stepsampler.SliceSampler(
        nsteps=2 * len(parameters),
        generate_direction=ultranest.stepsampler.generate_mixture_random_direction,
    )
    result = ddh_sampler.run(min_num_live_points=2000, max_ncalls=500_000)
    ddh_sampler.print_results()

[ ]:

if run_full_inference:
    import corner

    figure = corner.corner(
        result["samples"],
        labels=[r"$a_{\sigma}$", r"$a_{\omega}$", r"$a_{\rho}$",
                r"$\Gamma_{\sigma}$", r"$\Gamma_{\omega}$", r"$\Gamma_{\rho}$", r"$\rho_0$"],
        smooth=0.9,
        label_kwargs=dict(fontsize=22),
        title_kwargs=dict(fontsize=22),
        quantiles=[0.16, 0.84],
        title_quantiles=[0.16, 0.50, 0.84],
        levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.)),
        plot_density=False,
        plot_datapoints=False,
        fill_contours=True,
        show_titles=True,
        max_n_ticks=3,
        title_fmt=".2f",
    )