MeanField/experiments/NEGATIVE_DENSITY_INVESTIGATION_2026-09-19.md

Negative density: location, mechanism, and implications

Branch: research/negative-density-diagnosis-20260919. Parent: 1148c09 (native spherical harmonic surface validation). Scope: frozen-state observations and independent reconstructions only. No equilibrium equation, field order, surface parameterization, preconditioner, normalization, or Newton strategy was changed. No new equilibrium was solved in this investigation.

Main finding

At the measured converged n=1.5 endpoint, the negative density is an undershoot of the quadratic discontinuous density representation required by the weak EOS closure. An independently assembled elementwise projection of the nonnegative EOS target sampled by the closure reproduces the stored coefficients within 1.40e-14, while the negative density is of order 1e-4. This identifies a property of the measured discrete root, rather than an unfinished density-equation solve.

The negative region is in the outermost stellar cells. All 96 surface-adjacent cells contain it; all stored density DOF values remain positive. The initial n=1.5 seed is already negative there before Newton, through a related but distinct interpolation mechanism. A nonrotating n=3 seed has a negative pocket inside the same outer cell layer. A thinner negative layer also exists in the converged n=1 control, despite being missed by q18 and q26 volume sampling.

The practical implication is that a better preconditioner or tighter solve cannot remove the negative polynomial at this same root. We need to assess the density representation/closure before treating convergence as physical validation. This does not prove that the previous Newton plateaus, or the separate hydrostatic-gradient error, have the same cause.

What the code actually solves

Density and enthalpy are separate unknowns. The current defaults are discontinuous L2/Q2 density and continuous H1/Q3 enthalpy; the uniform order increment is zero in these runs. The surface is the existing L=4 all-m candidate, solved simultaneously with the interior in the saved equilibrium runs.

On each stellar element, the density equation is

[ R_i=\sum_q w_q,\varphi_i(q)\left[\rho_h(q)-\rho_{\rm EOS}(h_h(q))\right],\qquad \rho_{\rm EOS}(h)=\left[\frac{h}{(n+1)K}\right]^n\quad(h\ge0). ]

The weights include the physical deformation. Consequently, a zero density residual gives

[ M\boldsymbol\rho=\boldsymbol b(h),\qquad M_{ij}=\sum_q w_q\varphi_i\varphi_j,\quad b_i=\sum_q w_q\varphi_i\rho_{\rm EOS}(h). ]

This is a quadrature-defined weighted L2 projection. It is not pointwise EOS equality, and its unique polynomial need not be nonnegative. The local mass matrices were independently assembled and solved with Cholesky; no regularization or clipping was used.

In particular, the implemented dependency is not “negative density produces negative enthalpy.” Preparation requires finite density and rejects sampled negative enthalpy. Pressure is evaluated from enthalpy. Hydrostatic balance constrains h + phi - centrifugal_potential - C; density affects it indirectly through the coupled problem. At every negative-density q18/q26 volume sample in this study, enthalpy was positive. The surface traces have exactly zero enthalpy even where the density trace is negative.

Source locations:

• Independent density and enthalpy declarations, Density and Enthalpy.
• Prepared closure, state checks and weightedResidual assembly around lines 603–652; BuildResidual around line 681.
• Pressure preparation, pressure from enthalpy around line 787.
• Hydrostatic preparation, Bernoulli imbalance around line 863.
• Seed initialization, independent ProjectCoefficient calls around lines 164–178.

These relative links are repository references. The implementation paths and line numbers refer to the recorded diagnostic branch; the numerical findings are tied to immutable producer commits below.

Cases and controls

Four frozen states were examined on the existing modest mesh:

State	Origin	Central-density constraint	Purpose
n=1, J=0 equilibrium	Saved converged harmonic campaign	Absent	Lower-index control
n=1.5, J=0 seed	Original Lane–Emden initialization	Absent	Determine whether Newton creates the defect
n=1.5, J=0 equilibrium	Saved converged harmonic campaign	Absent	Identify the endpoint closure mechanism
n=3, J=0 seed	Original Lane–Emden initialization	Present	Assess relevance to the target index

The n=3 state is not a converged n=3 model, and no rotating n=3 model was tested. Its normalized seed residual is about 3.10e-3. The two saved equilibria retain their previous normalized residuals, approximately 1.53e-12 and 1.30e-12.

The independent projection used the production closure query and separate q18/q26 rules. The production operator's actual rules and assembled residual were exported separately. Independent versus production closure residual differences have L2 norms of roughly 1.2e-17 to 2.1e-17 across the four states. Constant reconstruction errors were also checked independently. These checks distinguish a correct projection reconstruction from an accidental reuse of the production inverse.

For the n=1.5 equilibrium, the production-rule projected density differs from the saved field by only 2.25e-15 in the physical L2 norm. Raising projection quadrature to q26 changes its coefficients by at most 2.36e-7; the negative density persists. Quadrature error in the EOS projection is therefore far too small to account for the measured undershoot.

Location and sampling limits

The q26 volume samples locate n=1.5 equilibrium negativity at physical radii approximately 0.99839–0.99993 in unit-reference-radius coordinates. All negative samples belong to the 96 cells touching the stellar/nonstellar interface; none occur in the other 736 stellar cells at these quadrature points. The surface is slightly deformed, so a global radius is not the same as distance from the local surface.

For the n=3 seed, q26 samples locate a negative pocket at radii approximately 0.93615–0.97070, with positive enthalpy about 0.0302–0.0682 there. Its surface density trace is positive even though the density becomes negative farther into the same cell. Thus checking only the surface or density DOFs would miss this defect.

State	Minimum density at q26 samples	Minimum density DOF value	Minimum density on sampled surface traces
n=1 equilibrium	+9.14281e-4	+1.70519e-2	-1.14780e-4
n=1.5 seed	-7.46708e-4	+2.90301e-3	-9.26358e-4
n=1.5 equilibrium	-5.59609e-4	+2.83809e-3	-7.21193e-4
n=3 seed	-8.48632e-5	+3.26982e-6	+2.49230e-4

All 2,400 sampled surface trace points are negative in the n=1 and n=1.5 cases; none are negative on the n=3 surface itself. These counts include coincident cell-edge points and are not surface-area fractions or global minimum certificates.

Ordinary quadrature of the sign indicator is not reliable for the volume of a thin negative layer. For the n=1.5 equilibrium, the negative mass estimate changes from 2.89e-5 at q18 to 2.18e-5 at q26; its estimated negative volume nearly halves. The field is unchanged. We therefore added a separate diagnostic that solves the density polynomial's zero crossings along reference-cell lines and integrates its negative intervals directly.

Integrated negative regions

Define negative-mass magnitude as integral_{rho<0} |rho| dV, divided below by the absolute total stellar mass. Volume fractions use the entire stellar volume. The finest interval results are:

Frozen state	Negative mass / stellar mass	Negative volume / stellar volume	Physical radius range of measured negative intervals / reference radius	Inward reference-cell depth range
n=1 equilibrium	6.09585e-8	0.0327346%	0.999643–0.999886	0–0.0007655
n=1.5 seed	3.66771e-5	1.83795%	0.993813–1.000019	0–0.0329556
n=1.5 equilibrium	2.24535e-5	1.49252%	0.994917–1.001222	0–0.0286204
n=3 seed	2.88372e-5	12.0835%	0.934409–0.978469	0.114925–0.349767

The inward coordinate is zero on each surface face and one on the opposite face of that cell; it is not a fraction of the whole stellar radius. Physical radii are measured from the center and do not identify the local surface location. Small deviations above unit radius arise in the represented surface, including the curved reference mesh.

Each state has measured negative intervals in all 96 surface-adjacent cells and none on the sampled lines through the other 736 stellar cells. The n=1/n=1.5 intervals reach the surface; the n=3 intervals form an interior pocket. All 24,576 exported finest negative intervals per state lie above the independent auditor's coefficient-uncertainty threshold.

The minimum-density witnesses in this integration are approximately -1.14756e-4, -9.26332e-4, -7.21171e-4, and -8.54870e-5 respectively. These are different tangential samples from the earlier surface-trace table, not certified global minima. Enthalpy at the first three witnesses is zero on the surface; at the n=3 pocket minimum it is +0.0455581.

Tangential rules with 8, 12, and 16 points per direction agree on negative integrals to roundoff in these cases; the two finest negative-mass integrals differ by at most 1.0e-14 relative. Full mass and volume agree with the independently sampled q18/q26 observations within 1.77e-13 relative. Each root-partitioned interval uses eight integration points. Eight analytic C++ controls and a separate decimal-arithmetic audit check roots, signs, minima, weights, sums, source identity, and provenance. The auditor also has four synthetic test methods, including deliberate corruptions.

This is strong convergence evidence for the measured integrals, not a rigorous three-dimensional sign certificate or a bound on every integration error. Tangential locations remain sampled; a sufficiently narrow angular pocket could be missed. Near-double roots below reconstruction uncertainty are not certified. The auditor verifies saved mapped integrals and their agreement with separate whole-volume observations; it is not a second FEM mapping implementation.

The lines connect saved samples along representative surface-cell traces. The n=1.5 endpoint and independent closure projection overlap; the n=3 panel is an initial seed. The displayed samples illustrate location and are not global extrema.

Why positive data can produce a negative polynomial

The density DOFs are values at interior Gauss–Legendre points. Positive values at those points do not imply that the interpolating polynomial is positive everywhere in the element. This explains why simply constraining the existing density coefficients would not address the observed defect.

The initial seed and the converged closure use different approximation operations:

1. The seed interpolates an independently tabulated, nonnegative Lane–Emden density into its nodal density space.
2. At equilibrium, the weak EOS closure determines the weighted L2 projection of the separately represented enthalpy's EOS density.

Both operations can undershoot. The nonrotating seeds establish that neither rotational distortion nor an accepted Newton update is necessary for the defect. Its occurrence throughout the outer cell layer points to a radial approximation mechanism, rather than an isolated bad angular mode. This does not establish independence from mesh quality or deformation in more difficult runs. In a simple one-dimensional surface cell, take positive enthalpy proportional to inward distance and normalize the density to t^n, with 0 <= t <= 1. The exact quadratic L2 projections are

[ p_{1.5}(t)=-\frac{2}{105}+\frac{16}{35}t+\frac47t^2, \qquad p_3(t)=\frac1{20}-\frac35t+\frac32t^2. ]

The first is negative for 0 <= t < 0.0396969. The cubic-density example has an interior negative pocket, 0.118350 < t < 0.281650, and minimum -0.01 at t=0.2. Both retain the exact cell mean. The n=3 example is a smooth cubic, so the problem is not restricted to fractional powers or a nonsmooth EOS derivative. Three-point Gauss–Legendre interpolation also undershoots: it agrees with the cubic example exactly, and produces a slightly deeper surface undershoot for n=1.5. All interpolated DOF values are positive.

This independently explains the two observed kinds of location, without Newton, gravity, surface deformation, rotation, or exterior-domain approximation. It is a mechanism control, not a quantitative prediction for the curved three-dimensional mesh. The n=1 leading-order example is linear and exact. In the actual discretization, Q3 enthalpy need not lie in the Q2 density space; its projection uses the curved physical-volume weights. The toy model omits that higher-order reference-cell variation. Geometry weights alone would not spoil exact reproduction of an EOS density already in Q2. The n=1 endpoint is nevertheless reproduced by the same independent projection to roundoff.

For this fixed-slope toy model, refining the surface-cell width by a factor of two reduces undershoot amplitude by 2^-n and negative mass per unit area by 2^-(n+1). It does not remove the negative fraction of that cell. Seven exact/algebraic self-tests verify the projections, roots, positive DOFs, mean preservation, and refinement laws.

Why virial balance does not exclude the defect

The constant test function in the density space preserves the EOS density's cell integral under the projection quadrature. Positive and negative pointwise errors can therefore cancel in mass even when they remain locally substantial. Pressure is computed from h, while gravitational and other density-dependent terms use the independent density representation. A small enthalpy-based virial error is compatible with this discrepancy.

The n=1.5 endpoint still has enthalpy-based energy and force virial errors approximately 6.55e-7 and 6.01e-7, respectively. Its density-based pressure and virial are undefined wherever density is negative. Neither successful algebraic convergence nor the enthalpy virial establishes pointwise EOS consistency. Negative-mass fractions describe the defect; they are not virial-error estimates and must not be compared directly with the 1e-5 virial tolerance.

The newly detected n=1 face layer also means that the earlier statement “no negative density at the sampled volume points” must remain qualified. Positive volume-quadrature observations are not a positivity certificate, and density-based virials computed from those samples do not prove that the complete field admits a physical density EOS everywhere.

Candidate resolutions, not implemented here

A small controlled order/refinement study is the next inexpensive discriminator. Change density order alone while retaining the surface, enthalpy, gravity, and displacement choices. Q3 density can represent Q3 enthalpy's linear n=1 EOS exactly on compatible reference cells, making the n=1 layer a clean control. For n=3, Q3 density exactly represents the leading linear-enthalpy surface model, but a generic Q3 enthalpy cubed is Q9; Q3 density is not a positivity guarantee. For n=1.5 the composition remains nonpolynomial. Local refinement should reduce the defect but does not generally enforce its sign.

There is no declared stability contract requiring density order to equal gravity-potential order: the existing RT/L2 pairing constraint belongs to gravity flux and potential. The available uniform-order build option would raise multiple fields together and confound a density-only study. Any eventual independent choice should remain in the field registry and preserve its compile-time checks.

Pointwise EOS-derived density is the preferred architectural candidate for the barotropic model. Evaluate physical rho(h) consistently wherever density is required, as pressure is already evaluated from enthalpy. This removes the particular independent-density projection undershoot, conditional on admissible enthalpy, and gives zero density wherever the surface enthalpy is zero. It also removes an independent field/closure pair. It is a formulation change: gravity sources and forces, centrifugal forces, mass, angular momentum/inertia, the central-density condition, diagnostics, and output must all use the same constitutive density and correct derivatives.

The compiler's state blocks and dependency graph, normalization and preconditioner blocks must change together through a typed formulation choice. Replacing density values in a cache while leaving density-column dependencies intact would be incorrect. This can remain an internal numerical choice without adding complexity to the user's EOS/mass/rotation specification. Enthalpy admissibility still needs to hold; composing the EOS does not by itself prevent an enthalpy polynomial from undershooting. Its convergence, cost, and virial quality require measurement; none is established by this investigation.

Some weak quantities could remain unchanged at fixed enthalpy and geometry: with matched quadrature, the existing projection already preserves the constant mass moment and any gravity-source moments whose test functions lie in the density space. The benefit is a consistent constitutive density throughout the model, not a guarantee that every source term changes or that virial balance automatically improves.

A positivity-constrained conservative projection is an alternative if an independent polynomial density is needed. A local weighted projection could enforce nonnegativity and preserve the constant EOS moment. It cannot retain every original Q2 closure moment when the unique unconstrained solution is negative. Nor does preserving cell mass automatically preserve inertia, forces or virial balance. The altered residual and derivatives through its active constraints would need consistent coupling; sampling positivity constraints only guarantees the chosen samples unless a stronger polynomial certificate is supplied.

Mean-preserving scaling constructions are well established for time-dependent high-order DG methods, but their time-step and admissibility hypotheses do not supply a convergence theorem for this steady stellar system. Zhang and Shu's review describes that construction. Xu and Shu's stationary analysis further illustrates why preserving cell averages alone can fail to preserve a stationary discrete balance. Those are methodological references, not ready-made stellar solvers.

Changes that would not resolve the identified mechanism

• A stronger preconditioner or tighter Newton tolerance would reproduce the same negative density at this root.
• Exact algebraic elimination rho = M(q)^-1 b(h,q) also reproduces that polynomial. Its correct derivative is delta rho = M^-1(delta b - delta M rho); eliminating a block does not make this a pointwise EOS formulation.
• A basis change alone leaves the represented polynomial unchanged. Nonnegative coefficients in the current nodal basis are already present.
• Clipping density in diagnostics would conceal the error. Clipping or rescaling it after solving would alter the solved equations and density-dependent physical quantities.

Implication for the n=3 rotating target

This is useful progress: the investigation isolates a concrete density-representation defect that exists before Newton and remains at a converged root, including a directly relevant n=3 control. It also exposes a validation blind spot in ordinary quadrature sampling. It does not yet establish a valid n=3 equilibrium, explain every nonlinear plateau, or solve the hydrostatic-gradient discrepancy.

Any chosen resolution still needs an accurate n=3 nonrotating anchor, the intended central-density normalization for its homology freedom, then a rotating sequence to the agreed moderate-rotation target. Global virial error <=1e-5 remains required, together with EOS-consistent internal fields and resolution checks suitable for comparison with ESTER.

Reproducibility and validation

Initial spatial/projection producer: 9c22934d52b369f96364e781d6a4843939ec9523. Root-interval producer: d1a594808f960bbaeeb6df8a8ce6a95458765b19. Evidence directory: ../../diagnostics/negative_density_2026-09-19 relative to this checkout.

• campaign/ contains four serial observations, immutable binary, source/build/mesh provenance, exact commands, input/output hashes, and retained CSVs.
• campaign_audit.json verifies the saved algebra, independent production/projection residual agreement, physical-observer agreement, sums, coefficient errors, and hash bindings.
• regions_campaign/ contains the same four states, observed serially with the root-interval diagnostic. All four completed successfully; no Newton steps were taken.
• regions_campaign_audit_v1.json and its log record the independent interval audit. Ninety-six consumed files are hash-bound. Observed state vectors, full density coefficients, and block layouts are byte-identical between the two campaigns.
• audit_self_checks.json verifies acceptance of a valid case and rejection of missing hash coverage and a nonmaximal coefficient corruption that does not change the reported maximum error.
• cli_scope_checks.json records rejection of accidental Newton execution and incomplete snapshot options.
• reference_model/ and reference_model.log record the independent one-dimensional controls and seven passing self-tests.
• figures/surface_density_traces.png displays representative stored-field and independent-projection traces; the n=3 panel is an initial seed, not a solution.

The exact case arguments are committed in negative_density_cases_2026-09-19.json and negative_density_region_cases_2026-09-19.json. Each campaign contains its private executable, exact command records, build settings and source commit. Both use sandbox.smesh, SHA-256 ce11ae99e21e6c3bbfbee402acd2e191c1da0d8261d2227b4203f73f2f337a74, with uniform field-order increment zero. Run audit_negative_density.py on campaign/ and audit_negative_density_regions.py on regions_campaign/; the latter binds automatically to its sibling campaign/. The analytic control and region auditor each expose --self-test.

The numerical observations are single-rank. No production test suite was rerun because production code was unchanged; the compiled diagnostic, independent mathematical controls, actual-operator residual comparisons, frozen-state invariants and artifact audits are the validation specific to this investigation.