MeanField_SERiF/docs/static/laneemden.md

Lane–Emden profiles and radial projection

serif::utils::laneemden solves the dimensionless, regular-center Lane–Emden equation,

[ \theta''+\frac{2}{\xi}\theta'+\theta^n=0, \qquad \theta(0)=1,\qquad\theta'(0)=0. ]

The dimensionless calculation is independent of MFEM, the equation of state, and any choice of physical units. Seed construction and normalization can therefore share it.

Dimensionless solution

#include "serif/utils/laneemden/laneemden.hpp"

const auto solution = serif::utils::laneemden::solve(3.0, {
    .step = 0.01,
    .maximum_radius = 100.0,
    .maximum_steps = 100000,
    .surface_tolerance = 1e-10
});

The result is std::expected<Solution, std::string>. Check success before dereferencing. A solution contains n, an ordered vector of Sample{xi, theta, derivative}, and an optional surface. The final sample is the first zero when found. Otherwise it is the requested maximum radius: this is a valid finite-interval solution, including the infinite-radius profiles with n >= 5. Exhausting the step budget is an error. evaluate(solution, xi) interpolates within the stored interval and rejects extrapolation.

The solver retains the central power-series recurrence from Boudreaux's original PolytropicStellarModel, evaluating the coefficients iteratively. Both theta and its derivative are initialized from that series and advance together through classical RK4. The final approach uses shorter steps and a short root bracket, avoiding nonreal powers of negative theta for fractional indices. surface_tolerance controls the final root localization, not the accumulated integration error. Reduce step to establish convergence. Fractional-index surface derivatives need not converge at fourth order, because the solution loses smoothness there. The regular analytic solutions at n = 0, 1, 5 are covered by the tests.

Choose a physical star

Three functions produce the same small PhysicalScales type:

Function	Physical inputs, besides the solution and explicit G
scale_from_mass	Total mass and central density
scale_from_enthalpy	Central specific enthalpy and central density
scale_from_polytropic_constant	Polytropic EOS constant K and central density

Each requires a found finite surface. The index and surface give mu = -xi_surface^2 * theta_prime_surface. Mass and central density then fix

[ \alpha=\left(\frac{M}{4\pi\rho_c\mu}\right)^{1/3},\qquad R=\alpha\xi_{\mathrm{surface}},\qquad h_c=4\pi G\rho_c\alpha^2. ]

For an EOS prescription, h_c = (n + 1) K rho_c^(1/n) instead fixes alpha. scale_from_polytropic_constant requires n > 0; the incompressible n = 0 limit can use either other prescription. K is a plain number because its dimensions depend on n. Other physical inputs use SERiF's distinct quantity types.

Choosing mass and central density creates a reference star and reports its implied K; it does not modify an existing EOS. To preserve a prescribed polytropic EOS, use its K. In particular, at n = 3 a fixed K and G already fix the mass: changing central density changes the radius, not that mass.

Use one consistent unit system throughout. SERiF does not convert units implicitly. For example, CGS inputs use grams, centimetres, seconds, density in grams per cubic centimetre, and G in cubic centimetres per gram per square second. Neither physical scaling nor integration contains a hardcoded gravitational constant.

evaluate_physical(solution, scales, radius) returns density, pressure, specific enthalpy, gravitational potential, enclosed mass, and gravity. Here gravity is positive outward dPhi/dr; inward gravitational acceleration is its negative. Potential is zero at infinity. Beyond the surface the evaluator uses vacuum values: zero density, pressure, and enthalpy, constant enclosed mass, Phi = -GM/r, and gravity = GM/r^2. At the surface it uses the vacuum side of the density jump for n = 0.

Project a physical profile

Projection accepts a radial callable returning either double or std::expected<double, std::string>. The latter propagates sampling failures. This complete function prepares a reference polytrope and projects its density:

#include "serif/utils/laneemden/laneemden.hpp"
#include "serif/utils/projection/projection.hpp"

/** @brief Construct a polytropic reference density on a supplied scalar H1 or L2 space. */
serif::utils::projection::backend::mfem::ProjectionResult polytropic_density(
    mfem::FiniteElementSpace& space,
    serif::dimensions::QuantityValue<serif::dimensions::Mass> mass,
    serif::dimensions::QuantityValue<serif::dimensions::MassDensity> central_density,
    serif::dimensions::QuantityValue<serif::dimensions::GravitationalConstant> gravitational_constant
) {
    namespace laneemden = serif::utils::laneemden;
    namespace projection = serif::utils::projection;
    namespace dimensions = serif::dimensions;

    const auto solution = laneemden::solve(3.0);
    if (!solution) {
        return std::unexpected(solution.error());
    }
    const auto scales = laneemden::scale_from_mass(
        *solution, mass, central_density, gravitational_constant
    );
    if (!scales) {
        return std::unexpected(scales.error());
    }
    const auto density = [&](double radius) -> projection::RadialValue {
        return laneemden::evaluate_physical(
            *solution, *scales, dimensions::QuantityValue<dimensions::Length>{radius}
        ).transform([](const laneemden::PhysicalSample& sample) {
            return sample.density.value();
        });
    };
    return projection::backend::mfem::project_scalar(space, density);
}

The returned grid function owns its coefficients and borrows the supplied finite-element space. That space and its mesh must remain alive. Projection completes synchronously, so the profile and scaling objects in the example do not need to outlive the call.

Use project_vector with the radial gravity amplitude for outward potential gradients in component-vector spaces. The projection handles MFEM component ordering. For a displaced star, options can supply .center = {x, y, z}. An optional .mapping = &evaluator samples the physical geometry through SERiF's mapping and applies the required scalar-integral pullback. The mapping must belong to the same mesh and be refreshed after its fields change.

These are the finite element's native coefficient projections, usually interpolation; hanging-node constraints are enforced by conforming restriction before publishing the field. They do not enforce a global mass constraint or perform a global L2 fit. There is no implicit rescaling of the mesh or adjustment of the sampled profile. Numerical sample failures discard the entire result. Incompatible space or mapping bindings throw a caller-error exception. Compactified infinity needs an explicit limiting-value prescription, as shown below.

Compactified infinity

Finite physical samples are still taken from the radial profile. At points identified by the mapping as infinity, scalar_at_infinity or vector_at_infinity can supply an explicit limit. Each callback receives the original mesh position and returns the limit already pulled back into mesh coordinates, before MFEM's own element mapping. Without a callback, projection rejects those points. A failed, non-finite, or wrong-dimensional limit also rejects the result.

For example, a star's vacuum potential tends to zero at infinity:

namespace projection = serif::utils::projection;
const projection::backend::mfem::RadialProjectionOptions exterior{
    .mapping = &evaluator,
    .scalar_at_infinity = [](const mfem::Vector&) -> projection::RadialValue {
        return 0.0;
    }
};

Pass exterior as the last argument to the projection call. A displaced exterior or a different asymptotic field must provide its own limit.

The callback is used only at points the mapper classifies as infinity. Finite points extremely close to that boundary still use the existing Jacobian determinant and inverse, whose conditioning can reduce accuracy.