Cosmology Grand Potential G Star S From Potential

/-- **CAPSTONE (g*s from the potential).** The present-day entropy density
equals `(2π²/45)·(43/11)·T_γ³` with the same reduced hypothesis list: the
effective entropy dof in the η_B dynamical prefactor now rests on FRW
continuity + the existence of a pressure potential + boundary data. -/
theorem gStarS_from_potential
    {P s : ℝ → ℝ} {T a Tν : ℝ → ℝ} {sT T' a' Tν' ρν' : ℝ → ℝ}
    {t₁ t₂ : ℝ} {T₁ Tγ αν : ℝ}
    (hTt : ∀ t, T t ≠ 0) (hαν : αν ≠ 0) (hTνt : ∀ t, Tν t ≠ 0)
    (ha₂ : a t₂ ≠ 0) (hTγ : Tγ ≠ 0)
    (hP : ∀ t, HasDerivAt P (s (T t)) (T t))
    (hs : ∀ t, HasDerivAt s (sT t) (T t))
    (hT : ∀ t, HasDerivAt T (T' t) t) (ha : ∀ t, HasDerivAt a (a' t) t)
    (hTν : ∀ t, HasDerivAt Tν (Tν' t) t)
    (hρν : ∀ t, HasDerivAt (fun u => αν * Tν u ^ 4) (ρν' t) t)
    (hcont : ∀ t, a t * (T t * sT t * T' t)
        = -3 * a' t * (energyOf P s (T t) + P (T t)))
    (hcontν : ∀ t, a t * ρν' t
        = -3 * a' t * (αν * Tν t ^ 4 + αν * Tν t ^ 4 / 3))
    (hbefore : s (T t₁) = NeutrinoDilution.radiationEntropy 2 4 T₁)
    (hafter : s (T t₂) = NeutrinoDilution.radiationEntropy 2 0 Tγ)
    (hshare : Tν t₁ = T₁) :
    NeutrinoDilution.radiationEntropy 2 0 Tγ
        + NeutrinoDilution.radiationEntropy 0 6 (Tν t₂)
      = 2 * π ^ 2 / 45 * ((EntropyPerPhoton.gStarS : ℚ) : ℝ) * Tγ ^ 3 :=
  NeutrinoDilution.total_entropy_eq_gStarS hTγ
    (dilution_from_potential hTt hαν hTνt ha₂ hTγ hP hs hT ha hTν hρν
      hcont hcontν hbefore hafter hshare)

/-- **Euler relation (derived).** `T·s = ρ + P` is an algebraic identity of
the Legendre-transform structure — not an independent equilibrium postulate. -/
theorem potential_euler (P s : ℝ → ℝ) (x : ℝ) :
    x * s x = energyOf P s x + P x := by
  simp only [energyOf]
  ring

/-- **Gibbs–Duhem relation (derived).** Along any temperature trajectory
`T(t)`, the pressure obeys `p′ = s(T)·T′`: this is the chain rule applied to
`s = dP/dT`, not an independent postulate. -/
theorem potential_gibbs_duhem
    {P s T : ℝ → ℝ} {T' t : ℝ}
    (hP : HasDerivAt P (s (T t)) (T t))
    (hT : HasDerivAt T T' t) :
    HasDerivAt (fun u => P (T u)) (s (T t) * T') t := by
  simpa [Function.comp] using hP.comp t hT

/-- **THEOREM (plasma pressure closed form).** The log-kernel integrals
(`π⁴/45`, `7π⁴/360`, both derived via Mellin transforms) collapse the
pressure to `P = (π²/90)·(g_B + (7/8)·g_F)·T⁴`.  The `7/8` is the same
fermionic weight that appears in entropy and energy — here it comes out of
the pressure channel independently. -/
theorem plasmaPressure_eq (gB gF T : ℝ) :
    plasmaPressure gB gF T = π ^ 2 / 90 * (gB + 7 / 8 * gF) * T ^ 4 := by
  unfold plasmaPressure
  rw [RadiationEntropyRelation.boseLog_integral_value,
    RadiationEntropyRelation.fermiLog_integral_value]
  have hπ : (π : ℝ) ≠ 0 := Real.pi_ne_zero
  field_simp
  ring

/-- **THEOREM (plasma energy closed form).** The energy integrals (`π⁴/15`,
`7π⁴/120`) give `ρ = (π²/30)·(g_B + (7/8)·g_F)·T⁴`. -/
theorem plasmaEnergy_eq (gB gF T : ℝ) :
    plasmaEnergy gB gF T = π ^ 2 / 30 * (gB + 7 / 8 * gF) * T ^ 4 := by
  unfold plasmaEnergy
  rw [FermionWeightIntegral.bose_integral_value,
    FermionWeightIntegral.fermi_integral_value]
  have hπ : (π : ℝ) ≠ 0 := Real.pi_ne_zero
  field_simp
  ring

/-- **THEOREM (radiation equation of state derived).** `p = ρ/3` is forced by
the grand-canonical integrals — it is not an input anywhere in the chain. -/
theorem plasma_eos (gB gF T : ℝ) :
    plasmaPressure gB gF T = plasmaEnergy gB gF T / 3 := by
  rw [plasmaEnergy_eq, plasmaPressure_eq]
  ring

/-- **THEOREM (the plasma is a potential fluid).** The temperature derivative
of the statistical-mechanical pressure is *exactly* the `radiationEntropy` of
`NeutrinoDilution` — the object previously built from the independent entropy
integrals `∫σ_B`, `∫σ_F`.  So `s = dP/dT` holds between two independently
derived statistical-mechanical quantities: the potential structure of the
plasma is a theorem, not a definition. -/
theorem plasmaPressure_potential (gB gF x : ℝ) :
    HasDerivAt (fun T => plasmaPressure gB gF T)
      (NeutrinoDilution.radiationEntropy gB gF x) x := by
  rw [NeutrinoDilution.radiationEntropy_eq]
  have hfun : (fun T => plasmaPressure gB gF T)
      = fun T => π ^ 2 / 90 * (gB + 7 / 8 * gF) * T ^ 4 :=
    funext fun T => plasmaPressure_eq gB gF T
  rw [hfun]
  have hpow : HasDerivAt (fun T : ℝ => T ^ 4) (4 * x ^ 3) x := by
    simpa using hasDerivAt_pow 4 x
  have h := hpow.const_mul (π ^ 2 / 90 * (gB + 7 / 8 * gF))
  have hval : π ^ 2 / 90 * (gB + 7 / 8 * gF) * (4 * x ^ 3)
      = 2 * π ^ 2 / 45 * (gB + 7 / 8 * gF) * x ^ 3 := by ring
  rw [hval] at h
  exact h

/-- **THEOREM (Legendre consistency).** The abstract Legendre energy
`T·s − P` built from the pressure potential and the entropy equals the
independently derived energy integral: `∫ t³/(eᵗ∓1)` agrees with
`4·∫ t²·logkernel − ∫ t²·logkernel`.  Two separate ensemble computations
meet — this is the internal consistency of the grand-canonical structure. -/
theorem plasma_energyOf (gB gF : ℝ) :
    energyOf (plasmaPressure gB gF) (NeutrinoDilution.radiationEntropy gB gF)
      = plasmaEnergy gB gF := by
  funext T
  simp only [energyOf]
  rw [NeutrinoDilution.radiationEntropy_eq, plasmaEnergy_eq, plasmaPressure_eq]
  ring