MeanField/experiments/DENSITY_ORDER_STUDY_2026-09-19.md

Independent density order: controlled study

Branch: research/density-order-study-20260919. Parent: 7ed4d23, the completed negative-density investigation. The separate log-coordinate assessment is on research/log-density-assessment-20260919 at 90580b2.

Finding: Q3 removes the resolved n=1 density undershoot, and Q4 substantially improves higher-index density accuracy. Neither provides positivity at higher n or resolves the n=3 coupled-solver stall in the common eight-step budget. All n=1 and n=1.5 controls converge and remain well virialized. Density remains an independent unknown. Log density was assessed mathematically and in small 1D controls, not implemented in the coupled stellar solver.

Question and controls

Does increasing only the independent density approximation order reduce the observed unphysical undershoots, improve internal EOS consistency, and preserve virial accuracy at acceptable coupled-solver cost?

The density field remains an independent unknown in every equilibrium solve. The frozen EOS projections are observations of approximation spaces; none is installed as a replacement field or used to alter a solver state.

The baseline is discontinuous tensor Q2 density. The alternatives are Q3 and Q4 in the same Gauss-Legendre nodal family. Enthalpy remains continuous Q3; gravity retains its RT2/Q2 pair; displacement remains Q3; the mesh and geometry are unchanged. The surface uses the existing degree-four spherical harmonics with all azimuthal orders and simultaneous interior/surface solution.

The active density block has 832 times (p+1)^3 coefficients: 22,464 at Q2, 53,248 at Q3 and 104,000 at Q4. Other field blocks retain their dimensions. The central-density constraint is present only in the n=3 cases.

Implementation and mathematical scope

MEAN_FIELD_DENSITY_POLYNOMIAL_ORDER selects the base density order at build time; its default remains two. The existing uniform order increment still applies, and is zero throughout this study. The field registry supplies the selected order to typed field spaces, quadrature form queries, adapters, normalization and compiled block layouts. The physics-facing API is unchanged.

Under the study's default policy the closure selects a rule covering both its density mass form and EOS source form. Their nominal degrees differ once density can have higher order than enthalpy. Selecting the greater resolved rule preserves policy overrides and the Q2 baseline. Adjacent requested orders may use the same MFEM Gauss rule, so a missing nominal bound does not by itself prove actual underintegration.

Q3 is a useful n=1 control because the linear EOS image of Q3 enthalpy lies in the density space. For n=3, a general Q3 enthalpy cubed is Q9. For n=1.5 it is nonpolynomial. Raising density order therefore does not promise positivity for the higher-index models.

There is also a discrete consistency benefit at Q3 and above. On each cell, Q3 enthalpy belongs to the DG density test space, even through the common curved VALUE mapping. If a represents that enthalpy in the density basis and R is the physical closure residual, the implemented polytropic EOS gives sum_q w_q h_h rho_h - (n+1) sum_q w_q P(h_h) = a^T R. A zero closure residual therefore enforces this enthalpy-weighted EOS moment cell by cell under the closure quadrature. Q2 does not generally enforce it. This algebraic property does not establish positivity or virial balance; other quadrature rules and finite residuals introduce their own defects, and the central-density border does not change this identity.

Frozen and coupled observations

Frozen observations retain the saved enthalpy and geometry from the n=1 and n=1.5 nonrotating equilibria and the n=3 initial seed. Each target space receives an independently assembled physical-volume-weighted L2 projection, with one common q26 quadrature, and a nodal interpolant of rho_EOS(h_h). The latter is explicitly not the original independently tabulated Lane-Emden seed.

The squared EOS error has a higher degree than the projection equations. A separate shared q30 evaluation checks it against q26 without changing or re-solving the projected fields. Both rules are also used on coupled final states. At fractional n this is a numerical sensitivity check rather than a polynomial exactness argument.

The coupled runs retain the same production quadrature policy, whose resolved rule can change when density order changes. Actual closure rules are saved for each case. The frozen common-rule controls isolate approximation order more strictly; the coupled comparison measures the complete consistent discretization, including its policy-selected quadrature adjustments.

Density order does not enter the pressure-force rule or its enthalpy admissibility guard: their field operands are enthalpy and the displacement gradient. Pressure-preparation rejection counts therefore measure changed reached states and coupled directions, rather than a density-order change in pressure sampling. The current preparation order evaluates displacement and pressure before the barotropic closure. Closure therefore cannot preempt that pressure stage in these runs; all rejection sources are nevertheless retained separately rather than interpreting a stage count as total admissibility.

Coupled runs use the original Lane-Emden initialization, the same normalization and default preconditioner, an FGMRES restart of 160, at most 200 linear iterations per step, relative linear tolerance 0.03, and at most eight Newton steps. Nonlinear relative tolerance is 1e-8 and absolute tolerance is zero. Cases cover n=1,1.5,3 without rotation and the established n=1,J=0.01 rotating control. A failure within this budget is a recorded trajectory outcome, not a proof that the discretization has no root.

The initializer tabulates Lane–Emden density and enthalpy before constructing either FEM field, then interpolates them independently into their spaces (libmeanfield/impl/seed/lane_emden.cpp, lines 231–236, and libmeanfield/impl/seed/stellar_equilibrium_projection.cpp, lines 164–180). Gravity uses the represented density; it does not replace enthalpy. Thus the production density seed is I_p(rho_LE_tabulated), whereas the frozen nodal and L2 controls are I_p(rho_EOS(h_h)) and Pi_p(rho_EOS(h_h)). Nonlinear EOS evaluation does not commute with either radial or FEM interpolation. Increasing density order does not remove the discrepancy between the tabulated density target and the EOS image of represented Q3 enthalpy. Source structure supports unchanged seeded enthalpy across these orders; no cross-run bitwise identity of the initial states is asserted.

For n=3, the existing central-density phase condition constrains central enthalpy to its EOS value at the requested central density; it does not directly pin the independent DG density's point value. Bordered convergence is assessed alongside the residual with artificial central forcing removed and the normalized size of that forcing. That unbordered diagnostic retains the scalar constraint rows. A small bordered residual or virial alone cannot establish that the artificial forcing has vanished.

All expensive builds, numerical cases and artifact audits run serially. Reported solver costs must exclude the post-solve observer and distinguish Newton steps, Krylov work, preparation and preconditioner refresh.

The normalized residual and its relative stopping threshold belong to each order's discrete system. Changing density order changes its test space, diagonal Riesz weights and initial interpolant as well as possibly its resolved quadrature. A smaller final residual across orders is not a physical-error ranking. Krylov totals also reflect different Newton trajectories, warm starts and per-action costs; they do not by themselves measure preconditioner quality. Physical diagnostics and measured Newton time are therefore reported alongside the iteration counts and density block sizes. Timings are single observations on this workstation, with serial workloads; they are not repeated benchmark means or estimates with error bars.

Density sign regions are integrated by isolating polynomial roots of degree two through four along reference-cell lines. The observer also checks surface traces, EOS mismatch, signed and negative mass, moments and unchanged source fields. It separates resolved negative witnesses from values at reconstruction uncertainty. Tangential sampling and mapped integration remain numerical observations, not global sign certificates.

Here M = integral rho_h dV is signed mass, M_- = -integral min(rho_h,0) dV, and V_- is the volume where density is negative. Reported fractions are M_-/abs(M) and V_-/V; the denominator of the mass fraction is not the integral of absolute density. EOS L2 error means sqrt(integral (rho_h-rho_EOS(h_h))^2 dV) in the model's units, without an additional relative normalization.

If the stored-state observer encounters negative enthalpy, complete EOS error integrals are marked undefined while density integrals remain reported. No density is clipped. Both virial definitions from pressure evaluated at enthalpy and, where admissible, at density are retained; the required virial tolerance remains 1e-5.

The summary enthalpy minimum and EOS-availability flag refer to integration nodes. Surface enthalpy samples and the enthalpy at exported negative-density witnesses are reported separately. Sliced mass and volume are cross-checked against the physical observer; first moments and inertia are only compared, because their additional coordinate factors need not be integrated exactly by the fixed eight-point normal rule. Tangential refinement alone cannot establish accuracy in that normal direction.

Results

Frozen approximation controls

The frozen campaign completed all three cases and its final artifact audit passed. The following uses the common-q26 L2 projection; EOS L2 errors use the separate unsliced q30 measurement. Negative measures use the finest root-sliced tangential rule. These are approximation controls, not new equilibrium solves.

Frozen model	Density order	M_- / abs(M)	Negative volume fraction	EOS L2 error
n=1 equilibrium	Q2	6.09585e-8	3.27346e-4	1.25266e-4
n=1 equilibrium	Q3	1.05e-31 (unresolved roundoff)	2.49e-16	1.56e-15
n=1 equilibrium	Q4	1.25e-31 (unresolved roundoff)	2.73e-16	1.60e-15
n=1.5 equilibrium	Q2	2.24299e-5	1.49181e-2	3.67965e-4
n=1.5 equilibrium	Q3	9.33602e-6	1.04476e-2	1.12325e-4
n=1.5 equilibrium	Q4	3.18362e-6	6.83440e-3	4.80663e-5
n=3 seed	Q2	2.92063e-5	1.21064e-1	1.06409e-2
n=3 seed	Q3	2.63486e-6	3.73200e-2	6.86223e-4
n=3 seed	Q4	1.02066e-8	1.69628e-2	3.84980e-5

The n=1 control matches the representability argument: Q3 and Q4 contain the linear EOS image of Q3 enthalpy, and no resolved negative witnesses remain. Higher-order density does not guarantee positivity for either higher-index case. The n=3 Q4 negative mass is about 2,862 times smaller than Q2, yet negative density still occupies about 1.70% of the sampled stellar volume.

The difference between n=1.5 and n=3 is consistent with surface regularity. Where enthalpy vanishes linearly, the fractional-index density behaves as distance-to-surface raised to 1.5 and has a singular second derivative. For n=3, cubing the represented enthalpy gives a polynomial Q9 target. Increasing density order can therefore improve the integer-index approximation very rapidly without resolving the conditioning of the coupled nonlinear equations.

At n=1.5, changing the error rule from q26 to q30 changes the Q4 EOS L2 norm by about 0.18%; the corresponding fixed eight-point root-sliced norm differs by about 2.2%. This is why the unsliced q30 values are used above. Tangential refinement gives stable negative masses but cannot establish convergence of the fractional EOS error in the normal direction. The integral-index norms agree across q26/q30 to arithmetic accuracy, except for relative comparisons of n=1 Q3/Q4 errors that are already at roundoff.

Coupled results

All nine n=1 and n=1.5 controls converged with full Newton steps and no rejected trials. Their energy and force virial errors from P(h) are all below 6.6e-7. The following times sum Newton iterations and exclude initialization and the post-solve density observer.

Model	Density	Newton updates	Krylov iterations	Newton seconds	M_- / abs(M)	EOS L2 error, q30
n=1, J=0	Q2	5	212	147.5	6.09585e-8	1.25266e-4
n=1, J=0	Q3	6	286	204.3	unresolved roundoff	3.13581e-16
n=1, J=0	Q4	6	286	215.1	unresolved roundoff	3.51021e-16
n=1, J=0.01	Q2	6	296	196.1	6.11256e-8	1.25327e-4
n=1, J=0.01	Q3	6	305	217.8	unresolved roundoff	7.26203e-16
n=1, J=0.01	Q4	6	304	229.1	unresolved roundoff	3.90198e-16
n=1.5, J=0	Q2	6	335	222.0	2.24535e-5	3.67965e-4
n=1.5, J=0	Q3	6	329	241.7	9.35454e-6	1.12195e-4
n=1.5, J=0	Q4	7	397	303.7	3.20115e-6	4.80640e-5

For n=1, Q3/Q4 remove the resolved undershoot and reduce the EOS mismatch to roundoff. Q3 and Q4 also agree in their virial diagnostics to about 1e-14. The nonrotating Q3/Q4 total times are 39%/46% above Q2, but both took an extra Newton step: their fifth residuals were 1.91e-12 and 1.80e-12 against thresholds of about 1.55e-12 and 1.54e-12. Q2's fifth residual narrowly passed its own threshold. Those total-time differences should not be interpreted as pure per-iteration order overhead. The rotating control takes six updates at every order; Q3/Q4 cost about 11%/17% more there.

For n=1.5, Q3 reduces negative mass by 2.40 times and EOS L2 error by 3.28 times for 8.9% more Newton time. Q4 reduces them by 7.01 and 7.66 times, respectively, for 36.8% more time. Its sixth residual, 2.26e-12, narrowly missed its 1.93e-12 threshold; the seventh update cost 58.8 seconds. Negative volume fractions remain 1.493%, 1.046% and 0.685% at Q2, Q3 and Q4. These are converged discrete roots with a resolved density defect, not merely incomplete Newton iterates.

Density/EOS consistency is not the only internal accuracy measure. At n=1.5, hydrostatic_gradient_rms is 0.00727475, 0.00731970 and 0.00729901 at Q2/Q3/Q4. It is the unscaled volume RMS of grad(h) + g_h - grad(Psi) in model acceleration units, using the independent gravity-gradient field g_h. It is neither a relative error nor the normalized algebraic residual. Thus density enrichment has not appreciably improved this local balance diagnostic.

None of the three n=3 cases converged within eight updates. Their endpoints must not be treated as equally accurate equilibrium models or their times as time to solution.

Density	F_final / F_initial	Smallest accepted alpha	Krylov iterations	Newton seconds	Energy virial, P(h)	Force virial, P(h)
Q2	0.430062	1/64	221	301.6	1.77163e-4	1.34597e-4
Q3	0.985301	1/256	452	478.8	8.07504e-6	2.28169e-5
Q4	0.987949	1/512	458	508.2	1.37586e-5	2.84372e-5

Q3 reaches the energy-virial milestone alone; its force virial and nonlinear residual still fail their respective criteria. Q4's better approximation space does not produce a better virial endpoint than Q3. Q2 recovers larger steps during this window, while Q3/Q4 make only about 1–1.5% residual progress. This is evidence of a stalled short trajectory, not a proof against eventual convergence with a longer run or different initialization.

All recorded n=3 linear solves meet their true-residual tolerance. Q2 has four pressure-force rejections and four insufficient-decrease rejections; Q3 and Q4 each have five and nine. No other rejection stages occur in these runs. Both enriched cases first accept only 1/256 of their Newton direction. The unchanged pressure quadrature makes this a change in reached states and directions, not a change in pressure sampling. Increasing density order alone has not removed pressure-admissibility and insufficient-decrease failures. This campaign records those failures but does not independently localize the nonlinear mechanism responsible for insufficient decrease.

The n=3 endpoint negative-mass fractions are 2.22117e-6, 9.13593e-7 and 5.21562e-7, with negative volume fractions 3.440%, 2.931% and 2.619%. EOS L2 errors are 9.78088e-3, 8.99701e-4 and 5.08742e-4. These still-unconverged fields show much less improvement than the frozen projections; their different trajectories prevent interpreting those ratios as equilibrium convergence rates.

The unbordered n=3 residuals remain 1.33327e-3, 7.18640e-4 and 7.26009e-4. Normalized central-border forcing is 3.17039e-5, 1.27669e-5 and 1.20966e-5. Thus neither a physical root nor vanishing artificial forcing has been demonstrated. Signed multipliers and both virial definitions are retained in the full tables.

The initializer distinction also matters. At n=3, the production Q4 seed has an observed EOS L2 mismatch of about 5.26e-4 under q18, while the frozen Q4 projection of represented enthalpy has 3.85e-5 under q30. These are different fields and quadratures, not a controlled improvement ratio. Together with the initializer's source structure, they show why refining the density space need not make the independently interpolated seed close to the weak closure manifold. A closure-consistent initial guess is a separate hypothesis; it could populate an independent density unknown without deriving it from enthalpy during the subsequent coupled solve.

The actual closure rules also differ only at Q4 in this study:

Model index	Q2 rule: achieved order / points per cell	Q3 rule	Q4 rule
n=1	17 / 729	17 / 729	19 / 1,000
n=1.5	19 / 1,000	19 / 1,000	21 / 1,331
n=3	23 / 1,728	23 / 1,728	25 / 2,197

All 832 stellar elements use the listed rule in each case. Q2/Q3 therefore share closure sampling, while Q4's coupled accuracy/cost includes the policy-selected quadrature enrichment as well as the larger space.

Implications and recommendation

Q3 is the economical consistency control for n=1: it contains the EOS image and Q4 adds no resolved accuracy benefit there. Q4 is a useful higher-index accuracy reference, but neither higher order is a positivity mechanism or a demonstrated n=3 convergence remedy. The production default remains Q2; these are isolated, reproducible research variants.

The next density-representation experiment should test a pointwise exponential with an independently chosen vacuum reference, retaining density as an independent unknown and solving surface and interior together. Judge it by EOS consistency, accepted Newton steps, both virial diagnostics and cost—not positivity alone. A separate conservative, closure-consistent seed control would distinguish initialization effects from representation effects. Neither of these further stellar experiments was performed here.

The complete tables and checked PNG/PDF figure are in ../../diagnostics/density_order_2026-09-19/final_analysis/ relative to the worktree. trial_stages.csv retains every source/disposition count, including accepted trials; coupled_outcomes.csv retains all endpoint diagnostics.

Comparison figure and coupled-result table.

Relation to the independent log-density hypothesis

The separate LOG_DENSITY_ASSESSMENT_2026-09-19.md in the log-density worktree contains the derivation, primary-source MESA comparison and reproducible one-dimensional controls. It does not contain a coupled stellar log solve.

Exponentiating polynomial nodal coefficients cannot remove the observed between-node undershoots: the current negative fields can already have positive nodal values. A pointwise representation rho = rho_scale * exp(u_h) changes the trial family and gives positive interior density in exact arithmetic while keeping density independent. A finite polynomial log field cannot also have an exact zero-density trace.

For the present vacuum polytrope, the more suitable hypothesis is rho = rho_ref * exp(u_h), where a reference chosen independently of live enthalpy is positive inside and vanishes at the surface with the correct leading exponent. That exponent is a polytropic assumption to test, not a universal prescription for the future general EOS. Geometry dependencies of the reference and all mapped operators must be differentiated consistently.

The log-coordinate closure derivative is a density-weighted mass matrix. Consequently normalization, local solves and globalization must account for its small response near vacuum. Exponentiation alone does not fix those conditioning issues, nor does it make the separately represented enthalpy admissible in every Newton trial. Floating-point underflow also needs an explicit rejection/range policy; silently clipping density would change the equations. A typed physical-density evaluator and its directional action would allow polynomial and log representations to be exchanged internally without exposing representation details in the physics-facing API.

Validation

All three builds passed the same 48 focused contract tests with the same RNG seed. They cover field registration and spaces, quadrature policy, prepared closure values and derivatives, gravity coupling, physical mass diagonals, and the normalized complete operator.

Density order	Test cases passed	Assertions passed
Q2	48	9,766
Q3	48	9,822
Q4	48	9,894

The final frozen artifact audit and all three coupled artifact audits pass. All twelve coupled cases completed with their declared outcome. The audit helper has nine passing small tests; the comparison helper has seven, including a check that separately audited orders retain the same mesh and non-density blocks. The independent log-coordinate control has four passing small tests.

No ESTER comparison or moderately rotating n=3 validation is part of this study. The n=1,J=0.01 case is a rotation control; its numerical J value is not a fraction of critical rotation.

Reproduction and retained evidence

Raw campaigns, immutable diagnostic executables, input/output hashes, build logs and test logs are retained under ../../diagnostics/density_order_2026-09-19/ relative to this worktree. Every campaign records its source commit, mesh hash, complete CMake cache, arguments and runtime environment. The mesh SHA256 is ce11ae99e21e6c3bbfbee402acd2e191c1da0d8261d2227b4203f73f2f337a74.

Producer commits are 241621c for frozen controls, c0e8551 for coupled Q2, cd71f37 for coupled Q3 and e054c74 for coupled Q4. Later changes to reporting scripts and this document do not replace those immutable executables or their recorded provenance. The final comparison generator is from c0a496e. The log control's source and CSV hashes were separately verified in log_reference_validation.json.

The committed presets are density_order_frozen_cases_2026-09-19.json and density_order_coupled_cases_2026-09-19.json in experiments/. Frozen replay requires the earlier harmonic-surface campaign retained alongside these data; the runner verifies the saved state, layout, metadata and original provenance.

For example, the Q3 configuration uses:

/Applications/CLion.app/Contents/bin/cmake/mac/aarch64/bin/cmake -S . -B cmake-build-release-density-q3 -G Ninja \
  -DCMAKE_BUILD_TYPE=Release \
  -DCMAKE_CXX_COMPILER=/opt/homebrew/opt/llvm/bin/clang++ \
  -DCMAKE_C_COMPILER=/opt/homebrew/opt/llvm/bin/clang \
  -DCMAKE_MAKE_PROGRAM=/opt/homebrew/bin/ninja \
  -DBoost_DIR=/opt/homebrew/anaconda3/lib/cmake/Boost-1.82.0 \
  -DMEAN_FIELD_DENSITY_POLYNOMIAL_ORDER=3 \
  -DMEAN_FIELD_UNIFORM_POLYNOMIAL_ORDER_INCREMENT=0 -DBUILD_TESTING=ON

Build plateau_diagnostic and density_order_contract_tests with two compile workers. The focused test filter includes all non-normalization cases in that target plus the reference-mass-diagonal comparison and complete normalized operator test:

env OMPI_MCA_btl=self,sm cmake-build-release-density-q3/density_order_contract_tests \
  --catch2 --reporter compact --rng-seed 1218739852 \
  '~[normalization],Direct Reference Mass Diagonals Match Full True DOF Assembly,Normalized Stellar Equilibrium Preserves Physical Residuals Jacobians And Frozen Solve Coordinates'

The same fixed seed is used for the other density orders. Campaign runtime settings are OMPI_MCA_btl=self,sm, OPENBLAS_NUM_THREADS=1, OMP_NUM_THREADS=1 and one MPI rank. The campaign runner uses a fresh output directory for each invocation; its --timeout 1800 applies per case. Use the order-specific build directory and the coupled preset with experiments/run_plateau_campaign.py; all other arguments are preserved in each provenance.json and *_run.json.

experiments/audit_density_order.py checks the retained evidence independently of the producer's summary tables. It validates saved polynomial algebra, bindings, sums, trial histories and declared gates, but cannot reconstruct the full FEM integrals because local matrices, RHS vectors and all mapped weights are not exported. experiments/summarize_density_order.py consumes only passing audit JSON files and emits comparison tables and scientific figures. It keeps undefined values and unresolved sign witnesses explicit.